%% Evaluate Forecasts
% This version is for the broader 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

if rel == 1 ;
    delete('BroadFF.out')
    diary('BroadFF.out ')
elseif rel == 2 ;
    delete('BroadAD.out')
    diary('BroadAD.out ')
elseif rel == 3 ;
    delete('BroadFV.out')
    diary('BroadFV.out ')
end ;

%% Load Data for Out-turns (1973Q1 - 2023Q4) & SPF Forecasts
% Forecasts are made at horizons of 0 to 4 quarters, from 1988Q1 to 2021Q2
% Each matrix has 212 rows: recession indicators for 1973Q1 - 2023Q4....
% ....and 213 columns for vintages 1971Q1 - 2024Q1 (missing data in columns for 1953Q2-1965Q3)

load('YSRep1973.mat','R1FF','R1AD','R1FV','SPF') ;

SPF = SPF(61:194,:) ;                                   % Remove SPF forecasts prior to 1988Q1 and after 2021Q2
T1 = 61 ;                                               % First period from which to forecast is 1988Q1
TF = 194 ;                                              % Final period from which to forecast is 2021Q2

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('BroadFF.mat','OosPR') ;                       % Load probit forecasts
elseif rel == 2 ;                                       % Use advance data
    load('BroadAD.mat','OosPR') ;                       % Load probit forecasts
elseif rel == 3 ;                                       % Use first final
    load('BroadFF.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
FCF = OosPR(:,7:11) ;                                   % RTNFCI forecasts

%% Calculate forecast errors & log-scores
% Forecasts errors and log-scores are for:
% h = 0: forecasts made in 1988Q1-2021Q2 for 1988Q1-2021Q2
% h = 1: forecasts made in 1988Q1-2021Q1 for 1988Q2-2021Q2
% h = 2: forecasts made in 1988Q1-2020Q4 for 1988Q3-2021Q2
% h = 3: forecasts made in 1988Q1-2020Q3 for 1988Q4-2021Q2
% h = 4: forecasts made in 1988Q1-2020Q2 for 1989Q1-2021Q2

NF = size(SPF,1) ;
FE = NaN(NF,5,3) ;              % Store forecast errors here, rows = qtrs, cols = horizon, pages = models
LS = NaN(NF,5,3) ;              % Store log scores here, rows = qtrs, cols = horizon, pages = models

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

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

%% Sample Ends 2007Q1 
 
% Calculate Forecast Statistics

diary off 
NF = 77 ;                                               % Evaluate forecasts upto those made in 2007Q1  
NM = size(FE,3) ;                                       % Number of models

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(NM-1,5) ;                                          % Store DM for MAE here
DM_MSE = NaN(NM-1,5) ;                                          % ditto MSE
DM_LPS = NaN(NM-1,5) ;                                          % ditto LPS

for m = 2:NM ;                                                  % Do loop through models
for h = 0:4 ;                                                   % Do loop through horizons

    diffMAE = abs(FE(1:NF-h,h+1,1)) - abs(FE(1:NF-h,h+1,m)) ;   % Difference MAE    
    diffMSE = FE(1:NF-h,h+1,1).^2 - FE(1:NF-h,h+1,m).^2 ;       % ditto MSE (next line is same for LPS)              
        if m == 2 ;         FCM = YSF ;                         % Use yield spread for log scores
    elseif m == 3 ;         FCM = FCF ;                         % or real time NFCI
    end ;
    diffLPS = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-FCM(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(FCM(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(m-1,h+1) =  bhatMAE/seMAE ;                          % DM stat - MAE
    DM_MSE(m-1,h+1) =  bhatMSE/seMSE ;                          % ditto MSE
    DM_LPS(m-1,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 ;
end ;
    
% Print Results

diary on
disp('1988Q1 - 2007Q1') ;
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)) ])
    disp(['      FCI ' num2str(round(FCSTAT(3,:),5)) ])
    clear FCSTAT
end ;
disp(blanks(1))

disp('Diebold-Mariano Statistics ')
disp('             MAE       MSE       LPS  ')
disp('------------------------------')
for h = 1:5 ;
    disp(['h = ' num2str(h-1) ])
    DMSTAT = [ DM_MAE(:,h)  DM_MSE(:,h) DM_LPS(:,h) ] ;
    disp(round(DMSTAT,5))
end ;

clear FCSTAT MAE RMSE LPS DMSTAT DM_MAE DM_MSE DM_LPS

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

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(NM-1,5) ;                                          % Store DM for MAE here
DM_MSE = NaN(NM-1,5) ;                                          % ditto MSE
DM_LPS = NaN(NM-1,5) ;                                          % ditto LPS

for m = 2:NM ;                                                  % Do loop through models
for h = 0:4 ;                                                   % Do loop through horizons

    diffMAE = abs(FE(1:NF-h,h+1,1)) - abs(FE(1:NF-h,h+1,m)) ;   % Difference MAE    
    diffMSE = FE(1:NF-h,h+1,1).^2 - FE(1:NF-h,h+1,m).^2 ;       % ditto MSE (next line is same for LPS)              
        if m == 2 ;         FCM = YSF ;                         % Use yield spread for log scores
    elseif m == 3 ;         FCM = FCF ;                         % or real time NFCI
    end ;
    diffLPS = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-FCM(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(FCM(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(m-1,h+1) =  bhatMAE/seMAE ;                          % DM stat - MAE
    DM_MSE(m-1,h+1) =  bhatMSE/seMSE ;                          % ditto MSE
    DM_LPS(m-1,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 ;
end ;

% Print Results

diary on
disp('1988Q1 - 2021Q2') ;
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)) ])
    disp(['      FCI ' num2str(round(FCSTAT(3,:),5)) ])
    clear FCSTAT
end ;
disp(blanks(1))

disp('Diebold-Mariano Statistics ')
disp('       h       MAE       MSE       LPS  ')
disp('----------------------------------------')
for h = 1:5 ;
    disp(['h = ' num2str(h-1) ])
    DMSTAT = [ DM_MAE(:,h)  DM_MSE(:,h) DM_LPS(:,h) ] ;
    disp(round(DMSTAT,5))
end ;

clear FCSTAT MAE RMSE LPS DMSTAT DM_MAE DM_MSE DM_LPS

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

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

if rel == 1 ;
    save('BroadFS.mat', 'SPF_FE', 'YSF_FE', 'FCI_FE', 'SPF_LS', 'YSF_LS', 'FCI_LS')
end ;
toc
