%% Generate Recession Forecasts
% This version is for the broader replication
clear all
tic 
rel = 2 ;                                           % Set to 1 to use first final data & 2 to use advance data

%% Load Data (1973Q1 - 2023Q4)
% First obs on R1 used in estimation is 1974Q2 
% First period from which to generate forecasts is 1988Q1 -- 61st obs
% Last period from which to generate forecasts is 2021Q2 -- 194th obs
load('YSRep1973.mat','R1FF','R1AD','Spread','NFCIFF','NFCIAD','qtr')

%% Generate OOS Forecasts -- expanding window for estimation

T1 = 61 ;                                           % First period from which to forecast is 1988Q1
TF = 194 ;                                          % Final period from which to forecast is 2021Q2
OosPR = [] ;                                        % Store predictions here

for t = T1:TF ;
    
    % Set up x matrices for different forecast horizons
    % Estimate R1(t-1) = a + b*X(t-2-h) 
    if rel == 1 ;
        y = R1FF(6:t-1,t) ;         % First final R1 1974Q2 to t-1 from vintage t
    elseif rel == 2 ;
        y = R1AD(6:t-1,t) ;         % Advance data R1 1974Q2 to t-1 from vintage t
    end ;

    x0s = Spread(5:t-2) ;           % Term spread for h = 0
    x1s = Spread(4:t-3) ;           % ditto, h = 1
    x2s = Spread(3:t-4) ;           % ditto, h = 2
    x3s = Spread(2:t-5) ;           % ditto, h = 3
    x4s = Spread(1:t-6) ;           % ditto, h = 4
    
    if rel == 1 ;                   % First final          
        x0f = NFCIFF(5:t-2,t) ;     % RTNFCI h = 0 from vintage t
        x1f = NFCIFF(4:t-3,t) ;     % ditto, h = 1
        x2f = NFCIFF(3:t-4,t) ;     % ditto, h = 2
        x3f = NFCIFF(2:t-5,t) ;     % ditto, h = 3
        x4f = NFCIFF(1:t-6,t) ;     % ditto, h = 4
    elseif rel == 2 ;               % Advance data
        x0f = NFCIAD(5:t-2,t) ;     % RTNFCI h = 0 from vintage t
        x1f = NFCIAD(4:t-3,t) ;     % ditto, h = 1
        x2f = NFCIAD(3:t-4,t) ;     % ditto, h = 2
        x3f = NFCIAD(2:t-5,t) ;     % ditto, h = 3
        x4f = NFCIAD(1:t-6,t) ;     % ditto, h = 4
    end ;

    % Estimate probits - yield spread 
    b0s = glmfit(x0s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % Estimate probit for h = 0
    b1s = glmfit(x1s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 1
    b2s = glmfit(x2s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 2
    b3s = glmfit(x3s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 3
    b4s = glmfit(x4s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 4

    % Estimate probits - real time NFCI 
    b0f = glmfit(x0f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % Estimate probit for h = 0
    b1f = glmfit(x1f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 1
    b2f = glmfit(x2f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 2
    b3f = glmfit(x3f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 3
    b4f = glmfit(x4f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto h = 4

    % Out of sample probit fitted values - yield spread 
    f0s = b0s(1) + b0s(2)*Spread(t-1) ;                 % Fitted value h = 0
    f1s = b1s(1) + b1s(2)*Spread(t-1) ;                 % ditto, h = 1 
    f2s = b2s(1) + b2s(2)*Spread(t-1) ;                 % ditto, h = 2
    f3s = b3s(1) + b3s(2)*Spread(t-1) ;                 % ditto, h = 3
    f4s = b4s(1) + b4s(2)*Spread(t-1) ;                 % ditto, h = 4

    % Out of sample probit fitted values - NFCI 

    if rel == 1 ;                                       % First final          
        f0f = b0f(1) + b0f(2)*NFCIFF(t-1,t) ;           % Fitted value h = 0
        f1f = b1f(1) + b1f(2)*NFCIFF(t-1,t) ;           % ditto, h = 1 
        f2f = b2f(1) + b2f(2)*NFCIFF(t-1,t) ;           % ditto, h = 2
        f3f = b3f(1) + b3f(2)*NFCIFF(t-1,t) ;           % ditto, h = 3
        f4f = b4f(1) + b4f(2)*NFCIFF(t-1,t) ;           % ditto, h = 4
    elseif rel == 2 ;                                   % Advance data
        f0f = b0f(1) + b0f(2)*NFCIAD(t-1,t) ;           % Fitted value h = 0
        f1f = b1f(1) + b1f(2)*NFCIAD(t-1,t) ;           % ditto, h = 1 
        f2f = b2f(1) + b2f(2)*NFCIAD(t-1,t) ;           % ditto, h = 2
        f3f = b3f(1) + b3f(2)*NFCIAD(t-1,t) ;           % ditto, h = 3
        f4f = b4f(1) + b4f(2)*NFCIAD(t-1,t) ;           % ditto, h = 4
    end ;

    % Predicted recession probabilities - yield spread 
    p0s = normcdf(f0s) ;                                % Prob R(t) = 1                         
    p1s = normcdf(f1s) ;                                % Prob R(t+1) = 1                           
    p2s = normcdf(f2s) ;                                % Prob R(t+2) = 1                                                          
    p3s = normcdf(f3s) ;                                % Prob R(t+3) = 1                     
    p4s = normcdf(f4s) ;                                % Prob R(t+4) = 1                        
   
    % Predicted recession probabilities - NFCI 
    p0f = normcdf(f0f) ;                                % Prob R(t) = 1                         
    p1f = normcdf(f1f) ;                                % Prob R(t+1) = 1                           
    p2f = normcdf(f2f) ;                                % Prob R(t+2) = 1                                                          
    p3f = normcdf(f3f) ;                                % Prob R(t+3) = 1                     
    p4f = normcdf(f4f) ;                                % Prob R(t+4) = 1                        

    % Store forecasts
    pt = [ qtr(t) p0s p1s p2s p3s p4s ] ;               % Time t forecasts - yield spread 
    pt = [ pt p0f p1f p2f p3f p4f ] ;                   % add real time NFCI forecasts
    OosPR = [ OosPR ; pt ] ;                            % Store time t date and forecasts

    clear y x0s x1s x2s x3s x4s b0s b1s b2s b3s b4s f0s f1s f2s f3s f4s p0s p1s p2s p3s p4s
    clear x0f x1f x2f x3f x4f b0f b1f b2f b3f b4f f0f f1f f2f f3f f4f p0f p1f p2f p3f p4f pt

end ;

% Write objects to file

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

toc
