%% Generate Recession Forecasts
% This version is for the narrow replication 

clear all
tic 
rel = 2 ;                                           % Set to 1 to use first final data & 2 to use advance data

%% Load Data (1953Q2 - 2023Q4)
% First obs on R1 used as dependent variable in estimation is 1955Q1 -- 8th obs
% First period from which to generate forecasts is 1968Q4 -- 63rd obs
% Last period from which to generate forecasts is 2023Q4 -- 283rd obs
load('YSRep.mat','R1FF','R1AD','Spread','qtr')

%% Generate OOS Forecasts -- expanding window for estimation

T1 = 63 ;                                           % First forecast origin is 1968Q4
TF = 283 ;                                          % Final forecast origin is 2023Q4
OosPR = [] ;                                        % Store predictions here

for t = T1:TF ;
    
    % Set up x matrices for different forecast horizons
    % and estimate R(t-1) = a + b*Spread(t-2-h)
    
    if rel == 1 ;
        y = R1FF(8:t-1,t) ;                                         % First Final R1 1955Q1 to t-1 from data vintage t
    elseif rel == 2 ;
        y = R1AD(8:t-1,t) ;                                         % Advance Data R1 1955Q1 to t-1 from data vintage t
    end ;

    x0 = Spread(7:t-2) ;                                            % Term spread for h = 0 probit
    x1 = Spread(6:t-3) ;                                            % ditto, h = 1
    x2 = Spread(5:t-4) ;                                            % ditto h = 2
    x3 = Spread(4:t-5) ;                                            % ditto, h = 3
    x4 = Spread(3:t-6) ;                                            % ditto, h = 4
    
    % Estimate probits
    b0 = glmfit(x0,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;   % Estimate probit for h = 0
    b1 = glmfit(x1,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;   % ditto, h = 1
    b2 = glmfit(x2,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;   % ditto, h = 2
    b3 = glmfit(x3,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;   % ditto, h = 3
    b4 = glmfit(x4,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;   % ditto, h = 4

    % Out of sample forecast for probit fitted values at time t
    f0 = b0(1) + b0(2)*Spread(t-1) ;                                % Fitted value for h = 0
    f1 = b1(1) + b1(2)*Spread(t-1) ;                                % ditto, h = 1 
    f2 = b2(1) + b2(2)*Spread(t-1) ;                                % ditto, h = 2
    f3 = b3(1) + b3(2)*Spread(t-1) ;                                % ditto, h = 3
    f4 = b4(1) + b4(2)*Spread(t-1) ;                                % ditto, h = 4

    % Predicted recession probabilities at time t
    p0 = normcdf(f0) ;                                              % Prob R(t) = 1                         
    p1 = normcdf(f1) ;                                              % Prob R(t+1) = 1                           
    p2 = normcdf(f2) ;                                              % Prob R(t+2) = 1                                                          
    p3 = normcdf(f3) ;                                              % Prob R(t+3) = 1                     
    p4 = normcdf(f4) ;                                              % Prob R(t+4) = 1                        

    % Store forecasts
    pt = [qtr(t) p0 p1 p2 p3 p4 ] ;                                 % All time t forecasts 
    OosPR = [ OosPR ; pt ] ;                                        % Store time t date and forecasts
   
    clear y x0 x1 x2 x3 x4 b0 b1 b2 b3 b4 f0 f1 f2 f3 f4 p0 p1 p2 p3 p4 pt

end ;

% Write objects to file

if rel == 1 ;
    save('NarrowFF.mat', 'OosPR')
elseif rel == 2 ;
    save('NarrowAD.mat', 'OosPR')
end ;

toc
