%% Generate Recession Forecasts
% This version uses the yield curve inversion indicator - first final R1 data only 
% (forecasts begin in 1975Q1, before then insufficient inversions for estimation)
clear all
tic 

%% Load Data (1953Q2 - 2023Q4)
% First obs on R1 used in estimation is 1953Q3 (for h=0) to 1954Q3 (for h=4)
% First period from which to generate forecasts is 1975Q1 -- 88th obs
% Last period from which to generate forecasts is 2023Q4 -- 283rd obs
load('YSRep.mat','R1FF','Spread','SpreadNeg','qtr')

%% Generate OOS Forecasts -- expanding window for estimation

T1 = 88 ;                                           % First period from which to forecast is 1975Q1
TF = 283 ;                                          % Final period from which to forecast is 2023Q4
OosPR = [] ;                                        % Store predictions here

for t = T1:TF ;

    % Set up x matrices for different forecast horizons
    % Estimate R(t-1) = a + b*X(t-2-h) 

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

    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

    x0i = SpreadNeg(5:t-2) ;                            % Inversion dummy for h = 0
    x1i = SpreadNeg(4:t-3) ;                            % ditto, h = 1
    x2i = SpreadNeg(3:t-4) ;                            % ditto, h = 2
    x3i = SpreadNeg(2:t-5) ;                            % ditto, h = 3
    x4i = SpreadNeg(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 - inversion 
    b0i = glmfit(x0i,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % Estimate probit for h = 0
    b1i = glmfit(x1i,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 1
    b2i = glmfit(x2i,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 2
    b3i = glmfit(x3i,y,'binomial','link','probit','Options',statset('MaxIter',1000)) ;      % ditto for h = 3
    b4i = glmfit(x4i,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 - inversion 
    f0i = b0i(1) + b0i(2)*SpreadNeg(t-1) ;              % Fitted value h = 0
    f1i = b1i(1) + b1i(2)*SpreadNeg(t-1) ;              % ditto, h = 1 
    f2i = b2i(1) + b2i(2)*SpreadNeg(t-1) ;              % ditto, h = 2
    f3i = b3i(1) + b3i(2)*SpreadNeg(t-1) ;              % ditto, h = 3
    f4i = b4i(1) + b4i(2)*SpreadNeg(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 - inversion 
    p0i = normcdf(f0i) ;                                % Prob R(t) = 1                         
    p1i = normcdf(f1i) ;                                % Prob R(t+1) = 1                           
    p2i = normcdf(f2i) ;                                % Prob R(t+2) = 1                                                          
    p3i = normcdf(f3i) ;                                % Prob R(t+3) = 1                     
    p4i = normcdf(f4i) ;                                % Prob R(t+4) = 1                        

    % Store forecasts
    pt = [ qtr(t) p0s p1s p2s p3s p4s ] ;               % Time t forecasts - yield spread 
    pt = [ pt p0i p1i p2i p3i p4i ] ;                   % add inversion dummy forecast

    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 x0i x1i x2i x3i x4i b0i b1i b2i b3i b4i f0i f1i f2i f3i f4i p0i p1i p2i p3i p4i pt

end ;

% Write objects to file

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

toc
