%% Generate Recession Forecasts
% This version uses final vintage NFCI (forecasts begin in 1982Q2) - first final R1 data only
clear all
tic 

%% Load Data (1971Q1 - 2023Q4)
% Data begins in 1971Q1
% First obs on R1 used in estimation is 1972Q2 
% First period from which to generate forecasts is 1982Q2 -- 46th obs
% Last period from which to generate forecasts is 2023Q4 -- 212th obs
load('YSRep1971.mat','R1FF','Spread','NFCI','qtr')

%% Generate OOS Forecasts -- expanding window for estimation

T1 = 46 ;                                           % First period from which to forecast is 1982Q2
TF = 212 ;                                          % Final period from which to forecast is 2023Q4
OosPR = [] ;                                        % Store predictions here

for t = T1:TF ;

    % Set up x matrices for different forecast horizons
    % Using first final data and so in mid quarter the current vintage contains data on R1 from t-1
    % Estimate R(t-1) = a + b*Spread(t-2-h)

    y = R1FF(6:t-1,t) ;             % R1 indicator 1972Q2 to t-1 from vintage t

    x0s = Spread(5:t-2) ;           % Yield spread 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

    x0f = NFCI(5:t-2) ;             % NFCI h = 0
    x1f = NFCI(4:t-3) ;             % ditto, h = 1
    x2f = NFCI(3:t-4) ;             % ditto, h = 2
    x3f = NFCI(2:t-5) ;             % ditto, h = 3
    x4f = NFCI(1:t-6) ;             % ditto, h = 4

    % 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 for h = 1
    b2s = glmfit(x2s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 2
    b3s = glmfit(x3s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 3
    b4s = glmfit(x4s,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 4

   % Estimate probits - 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 for h = 1
    b2f = glmfit(x2f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 2
    b3f = glmfit(x3f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 3
    b4f = glmfit(x4f,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for 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 
    f0f = b0f(1) + b0f(2)*NFCI(t-1) ;                   % Fitted value h = 0
    f1f = b1f(1) + b1f(2)*NFCI(t-1) ;                   % ditto, h = 1 
    f2f = b2f(1) + b2f(2)*NFCI(t-1) ;                   % ditto, h = 2
    f3f = b3f(1) + b3f(2)*NFCI(t-1) ;                   % ditto, h = 3
    f4f = b4f(1) + b4f(2)*NFCI(t-1) ;                   % ditto, h = 4

    % 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 NFCI 
    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

save('FinalNFCI.mat', 'OosPR')

toc
