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

clear all
tic 
rel = 3 ;                                               % Set to 1 to use first final data as outturn, 2 to use advance data, 3 to use final vintage

%% Set output file
if rel == 1 ;
    delete('NarrowFF.out')
    diary('NarrowFF.out ')
elseif rel == 2 ;
    delete('NarrowAD.out')
    diary('NarrowAD.out ')
elseif rel == 3 ;
    delete('NarrowFV.out')
    diary('NarrowFV.out ')
end ;

%% Load Data for Out-turns (1953Q2 - 2023Q4) & Forecasts
% Forecasts are made at horizons of 0 to 4 quarters, from 1968Q4 to 2023Q4
% Each matrix has 283 rows: recession indicators for 1953Q2 - 2023Q4....
% ....and 284 columns for vintages 1953Q2 - 2024Q1 (missing data for 1953Q2-1965Q3)

load('YSRep.mat','R1FF','R1AD','R1FV','SPF') ;          % Data for out-turn and SPF

T1 = 63 ;                                               % First forecast origin is 1968Q4
TF = 283 ;                                              % Final forecast origin is 2023Q4  

% Create Out-turns

OutTurn = [] ;                                          % Store out-turn here
if rel == 1 ;                                           % Use first final
    for t = T1:TF ;
        OutTurn = [ OutTurn ; R1FF(t,t+1) ] ;           % Out-turn is first-final release for quarter t from 3rd month of t+1 quarter
    end ;
    clear R1FF t
elseif rel == 2 ;                                       % Use advance data
    for t = T1:TF ;
        OutTurn = [ OutTurn ; R1AD(t,t+1) ] ;           % Out-turn is advance-data release for quarter t from 2nd month of t+1 quarter
    end ;
    clear R1AD t
elseif rel == 3 ;                                       % Use first final
    OutTurn = R1FV(T1:TF) ;                             % Out-turn is final vintage release for quarter t
    clear R1FV t
end ;

% Load Probit Forecasts

if rel == 1 ;                                           % Use first final
    load('NarrowFF.mat','OosPR') ;                      % Load probit forecasts
elseif rel == 2 ;                                       % Use advance data
    load('NarrowAD.mat','OosPR') ;                      % Load probit forecasts
elseif rel == 3 ;                                       % Use first final
    load('NarrowFF.mat','OosPR') ;                      % Load probit forecasts (note - first final forecasts evaluated against final vintage outturns)
end ;
qtr = OosPR(:,1) ;                                      % Quarters
YSF = OosPR(:,2:6) ;                                    % Yield spread forecasts
clear OosPR

%% Calculate forecast errors & log-scores
% Forecasts errors and log-scores are for:
% h = 0: forecasts made in 1968Q4-2023Q4 for 1968Q4-2023Q4
% h = 1: forecasts made in 1968Q4-2023Q3 for 1969Q1-2023Q4
% h = 2: forecasts made in 1968Q4-2023Q2 for 1969Q2-2023Q4
% h = 3: forecasts made in 1968Q4-2023Q1 for 1969Q3-2023Q4
% h = 4: forecasts made in 1968Q4-2022Q4 for 1969Q4-2023Q4

NF = size(SPF,1) ;              % Number forecasts
FE = NaN(NF,5,2) ;              % Store forecast errors here, rows = qtrs, cols = horizon, pages = SPF/YSF
LS = NaN(NF,5,2) ;              % Store log scores here, rows = qtrs, cols = horizon, pages = SPF/YSF

% SPF
for h = 0:4 ;
    FE(1:NF-h,h+1,1) = SPF(1:NF-h,h+1) - OutTurn(h+1:NF) ;                         
    LS(1:NF-h,h+1,1) = (1-OutTurn(h+1:NF)).*log(1-SPF(1:NF-h,h+1)) + OutTurn(h+1:NF).*log(SPF(1:NF-h,h+1)) ;
end ;

% Yield Spread 
for h = 0:4 ;
    FE(1:NF-h,h+1,2) = YSF(1:NF-h,h+1) - OutTurn(h+1:NF) ;                         
    LS(1:NF-h,h+1,2) = (1-OutTurn(h+1:NF)).*log(1-YSF(1:NF-h,h+1)) + OutTurn(h+1:NF).*log(YSF(1:NF-h,h+1)) ;
end ;

%% RW's Sample

% Calculate Forecast Statistics 

diary off 
NF = 154 ;                                              % Evaluate forecasts upto 2007Q1  
NM = size(FE,3) ;                                       % # models is 3rd dimension of FE matrix

MAE = [] ;                                              % Store mean absolute error here
RMSE = [] ;                                             % Store root mean squared error here
LPS = [] ;                                              % Store log probability score here

for h = 0:4 ;                                           % Do loop through horizons

    MAEh = NaN(NM,1) ;                                  % Store horizon h MAE statistics here
    RMSEh = NaN(NM,1) ;                                 % ditto RMSE
    LPSh = NaN(NM,1) ;                                  % ditto LPS

    for m = 1:NM ;                                      % Do loop through models
        MAEh(m) = mean(abs(FE(1:NF-h,h+1,m))) ;         % MAE for model m at horizon h
        RMSEh(m) = sqrt(mean(FE(1:NF-h,h+1,m).^2)) ;    % ditto RMSE
        LPSh(m) = -mean(LS(1:NF-h,h+1,m)) ;             % ditto LPS
    end ;
    MAE = [ MAE MAEh ] ;                                % Store horizon h MAEs
    RMSE = [ RMSE RMSEh ] ;                             % ditto RMSE
    LPS = [ LPS LPSh ] ;                                % ditto LPS
    clear MAEh RMSEh LPSh m

end ;

% Diebold-Mariano Statistics Versus SPF

DM_MAE = NaN(5,1) ;                                             % Store DM for MAE here
DM_MSE = NaN(5,1) ;                                             % ditto MSE
DM_LPS = NaN(5,1) ;                                             % ditto LPS

for h = 0:4 ;                                                   % Do loop through horizons

    diffMAE = abs(FE(1:NF-h,h+1,1)) - abs(FE(1:NF-h,h+1,2)) ;   % Difference MAE    
    diffMSE = FE(1:NF-h,h+1,1).^2 - FE(1:NF-h,h+1,2).^2 ;       % ditto MSE (next line is same for LPS)              
    diffLPS = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-YSF(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(YSF(1:NF-h,h+1))) ;

    x = ones(length(diffMAE),1) ;                               % Intercept

    [ ~, seMAE, bhatMAE ] = hac(x,diffMAE,Intercept=false) ;    % OLS with HAC standard errors - MAE
    [ ~, seMSE, bhatMSE ] = hac(x,diffMSE,Intercept=false) ;    % ditto MSE
    [ ~, seLPS, bhatLPS ] = hac(x,diffLPS,Intercept=false) ;    % ditto LPS               

    DM_MAE(h+1) =  bhatMAE/seMAE ;                              % DM stat - MAE
    DM_MSE(h+1) =  bhatMSE/seMSE ;                              % ditto MSE
    DM_LPS(h+1) =  -bhatLPS/seLPS ;                             % ditto LPS (note -bhat, so postive stat implies improvement on SPF)

    clear diffMAE diffMSE diffLPS x seMAE seMSE seLPS bhatMAE bhatMSE bhatLPS

end ;

% Print Results

diary on
disp('RW Sample 1968Q4 - 2007Q1') ;
disp('*************************') ;
disp(blanks(1))
disp([ 'Last quarter from which forecasts are generated = ' num2str(qtr(NF)) ])
disp(blanks(1))
disp('           MAE        RMSE        LPS  ')
disp('---------------------------------------')
for h = 0:4 ;
    FCSTAT = [ MAE(:,h+1)  RMSE(:,h+1)  LPS(:,h+1) ] ;
    disp([ 'h = ' num2str(h) ])
    disp(['      SPF ' num2str(round(FCSTAT(1,:),5)) ])
    disp(['      YSF ' num2str(round(FCSTAT(2,:),5)) ])
    clear FCSTAT
end ;
disp(blanks(1))

disp('Diebold-Mariano Statistics ')
disp('       h       MAE       MSE       LPS  ')
disp('----------------------------------------')
H = [ 0 ; 1 ; 2 ; 3 ; 4 ] ;
DMSTAT = [ H DM_MAE  DM_MSE DM_LPS ] ;
disp(round(DMSTAT,5))
disp(blanks(2))
clear FCSTAT MAE RMSE LPS DMSTAT DM_MAE DM_MSE DM_LPS

%% CM's Sample, first final only

if rel == 1 ;

    % Calculate Forecast Statistics 

    diary off 
    NF = 178 ;                                              % Evaluate forecasts upto those made in 2013Q1  
    NM = size(FE,3) ;                                       % Number of models is 3rd dimension of FE matrix

    MAE = [] ;                                              % Store mean absolute error here
    RMSE = [] ;                                             % Store root mean squared error here
    LPS = [] ;                                              % Store log probability score here

    for h = 0:4 ;                                           % Do loop through horizons

        MAEh = NaN(NM,1) ;                                  % Store horizon h MAE statistics here
        RMSEh = NaN(NM,1) ;                                 % ditto RMSE
        LPSh = NaN(NM,1) ;                                  % ditto LPS

        for m = 1:NM ;                                      % Do loop through models
            MAEh(m) = mean(abs(FE(1:NF-h,h+1,m))) ;         % MAE for model m at horizon h
            RMSEh(m) = sqrt(mean(FE(1:NF-h,h+1,m).^2)) ;    % ditto RMSE
            LPSh(m) = -mean(LS(1:NF-h,h+1,m)) ;             % ditto LPS
        end ;
        MAE = [ MAE MAEh ] ;                                % Store horizon h MAEs
        RMSE = [ RMSE RMSEh ] ;                             % ditto RMSE
        LPS = [ LPS LPSh ] ;                                % ditto LPS
        clear MAEh RMSEh LPSh m

    end ;

    % Diebold-Mariano Statistics Versus SPF

    DM_MAE = NaN(5,1) ;                                             % Store DM for MAE here
    DM_MSE = NaN(5,1) ;                                             % ditto MSE
    DM_LPS = NaN(5,1) ;                                             % ditto LPS

    for h = 0:4 ;                                                   % Do loop through horizons

        diffMAE = abs(FE(1:NF-h,h+1,1)) - abs(FE(1:NF-h,h+1,2)) ;   % Difference MAE    
        diffMSE = FE(1:NF-h,h+1,1).^2 - FE(1:NF-h,h+1,2).^2 ;       % ditto MSE (next line is same for LPS)              
        diffLPS = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-YSF(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(YSF(1:NF-h,h+1))) ;

        x = ones(length(diffMAE),1) ;                               % Intercept

        [ ~, seMAE, bhatMAE ] = hac(x,diffMAE,Intercept=false) ;    % OLS with HAC standard errors - MAE
        [ ~, seMSE, bhatMSE ] = hac(x,diffMSE,Intercept=false) ;    % ditto MSE
        [ ~, seLPS, bhatLPS ] = hac(x,diffLPS,Intercept=false) ;    % ditto LPS               

        DM_MAE(h+1) =  bhatMAE/seMAE ;                              % DM stat - MAE
        DM_MSE(h+1) =  bhatMSE/seMSE ;                              % ditto MSE
        DM_LPS(h+1) =  -bhatLPS/seLPS ;                             % ditto LPS (note -bhat, so postive stat implies improvement on SPF)

        clear diffMAE diffMSE diffLPS x seMAE seMSE seLPS bhatMAE bhatMSE bhatLPS

    end ;

    % Print Results

    diary on
    disp('CM Sample 1968Q4 - 2013Q1') ;
    disp('*************************') ;
    disp(blanks(1))
    disp([ 'Last quarter from which forecasts are generated = ' num2str(qtr(NF)) ])
    disp(blanks(1))
    for h = 0:4 ;
        FCSTAT = [ MAE(:,h+1)  RMSE(:,h+1)  LPS(:,h+1) ] ;
        disp([ 'h = ' num2str(h) ])
        disp(['      SPF ' num2str(round(FCSTAT(1,:),5)) ])
        disp(['      YSF ' num2str(round(FCSTAT(2,:),5)) ])
        clear FCSTAT
    end ;
    disp(blanks(1))

    disp('Diebold-Mariano Statistics ')
    disp('       h       MAE       MSE       LPS  ')
    disp('----------------------------------------')
    H = [ 0 ; 1 ; 2 ; 3 ; 4 ] ;
    DMSTAT = [ H DM_MAE  DM_MSE DM_LPS ] ;
    disp(round(DMSTAT,5))
    disp(blanks(2))
    clear FCSTAT MAE RMSE LPS DMSTAT DM_MAE DM_MSE DM_LPS

end ;

%% Full Sample 

% Calculate Forecast Statistics 

diary off 
NF = size(SPF,1) ;                                      % Evaluate forecasts upto those made in 2022Q4  
NM = size(FE,3) ;                                       % Number of models is 3rd dimension of FE matrix

MAE = [] ;                                              % Store mean absolute error here
RMSE = [] ;                                             % Store root mean squared error here
LPS = [] ;                                              % Store log probability score here

for h = 0:4 ;                                           % Do loop through horizons

    MAEh = NaN(NM,1) ;                                  % Store horizon h MAE statistics here
    RMSEh = NaN(NM,1) ;                                 % ditto RMSE
    LPSh = NaN(NM,1) ;                                  % ditto LPS

    for m = 1:NM ;                                      % Do loop through models
        MAEh(m) = mean(abs(FE(1:NF-h,h+1,m))) ;         % MAE for model m at horizon h
        RMSEh(m) = sqrt(mean(FE(1:NF-h,h+1,m).^2)) ;    % ditto RMSE
        LPSh(m) = -mean(LS(1:NF-h,h+1,m)) ;             % ditto LPS
    end ;
    MAE = [ MAE MAEh ] ;                                % Store horizon h MAEs
    RMSE = [ RMSE RMSEh ] ;                             % ditto RMSE
    LPS = [ LPS LPSh ] ;                                % ditto LPS
    clear MAEh RMSEh LPSh m

end ;

% Diebold-Mariano Statistics Versus SPF

DM_MAE = NaN(5,1) ;                                             % Store DM for MAE here
DM_MSE = NaN(5,1) ;                                             % ditto MSE
DM_LPS = NaN(5,1) ;                                             % ditto LPS

for h = 0:4 ;                                                   % Do loop through horizons

    diffMAE = abs(FE(1:NF-h,h+1,1)) - abs(FE(1:NF-h,h+1,2)) ;   % Difference MAE    
    diffMSE = FE(1:NF-h,h+1,1).^2 - FE(1:NF-h,h+1,2).^2 ;       % ditto MSE (next line is same for LPS)              
    diffLPS = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-YSF(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(YSF(1:NF-h,h+1))) ;

    x = ones(length(diffMAE),1) ;                               % Intercept

    [ ~, seMAE, bhatMAE ] = hac(x,diffMAE,Intercept=false) ;    % OLS with HAC standard errors - MAE
    [ ~, seMSE, bhatMSE ] = hac(x,diffMSE,Intercept=false) ;    % ditto MSE
    [ ~, seLPS, bhatLPS ] = hac(x,diffLPS,Intercept=false) ;    % ditto LPS               

    DM_MAE(h+1) =  bhatMAE/seMAE ;                              % DM stat - MAE
    DM_MSE(h+1) =  bhatMSE/seMSE ;                              % ditto MSE
    DM_LPS(h+1) =  -bhatLPS/seLPS ;                             % ditto LPS (note -bhat, so postive stat implies improvement on SPF)

    clear diffMAE diffMSE diffLPS x seMAE seMSE seLPS bhatMAE bhatMSE bhatLPS

end ;

% Print Results for Full sample
diary on
disp('Full Sample 1968Q4 - 2023Q4') ;
disp('***************************') ;
disp(blanks(1))
disp([ 'Last quarter from which forecasts are generated = ' num2str(qtr(NF)) ])
disp(blanks(1))
for h = 0:4 ;
    FCSTAT = [ MAE(:,h+1)  RMSE(:,h+1)  LPS(:,h+1) ] ;
    disp([ 'h = ' num2str(h) ])
    disp(['      SPF ' num2str(round(FCSTAT(1,:),5)) ])
    disp(['      YSF ' num2str(round(FCSTAT(2,:),5)) ])
    clear FCSTAT
end ;
disp(blanks(1))

disp('Diebold-Mariano Statistics ')
disp('       h       MAE       MSE       LPS  ')
disp('----------------------------------------')
H = [ 0 ; 1 ; 2 ; 3 ; 4 ] ;
DMSTAT = [ H DM_MAE  DM_MSE DM_LPS ] ;
disp(round(DMSTAT,5))
disp(blanks(2))
clear FCSTAT MAE RMSE LPS DMSTAT DM_MAE DM_MSE DM_LPS

%% Write forecast errors and log-scores to file for first-final version only (for h = 1, 2023Q4 = NaN, h=2 2023Q3 & 2023Q4 = NaN etc)

SPF_FE = FE(:,:,1) ;                                    % SPF forecast errors
YSF_FE = FE(:,:,2) ;                                    % YSF forecast errors 
SPF_LS = LS(:,:,1) ;                                    % SPF log scores
YSF_LS = LS(:,:,2) ;                                    % YSF log scores

if rel == 1 ;
    save('NarrowFS.mat', 'SPF_FE', 'YSF_FE', 'SPF_LS', 'YSF_LS')
end ;

toc
