%% Evaluate 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 
delete('Inversion.out')
diary('Inversion.out')

%% Load Data (1953Q2 - 2023Q4)for Out-turns & SPF Forecasts
% Forecasts are made at horizons of 0 to 4 quarters, from 1974Q1 to 2023Q4
% R1FF is vintage from last month of quarter - recession indicators for 1953Q2 - 2023Q4 (279 rows)....
% ....and 284 columns for vintages 1953Q2 - 2024Q1 (missing data in columns for 1953Q2-1965Q3)

load('YSRep.mat','qtr','R1FF','SPF') ;

SPF = SPF(26:end,:) ;                                   % Remove SPF forecasts prior to 1975Q1
T1 = 88 ;                                               % First period from which to forecast is 1975Q1
TF = 283 ;                                              % Final period from which to project is 2023Q4  

OutTurn = [] ;
for t = T1:TF ;
    OutTurn = [ OutTurn ; R1FF(t,t+1) ] ;               % Out-turn is first final release for quarter t from 3rd month t+1 quarter
end ;
clear R1FF t

load('Inversion.mat','OosPR') ;
qtr = OosPR(:,1) ;                                      % Quarters
YSF = OosPR(:,2:6) ;                                    % Yield spread forecasts
YIF = OosPR(:,7:11) ;                                   % Yield inversion forecasts

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

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 ;

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


%% RW's sample

% Calculate Forecast Statistics 

NF = 129 ;                                              % Number of forecasts to 2007Q1
NM = size(FE,3) ;                                       % Number of models
diary off 

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 ;                         % Yield spread
    elseif m == 3 ;         FCM = YIF ;                         % Yield inversion
    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 Output

diary on
disp('1975Q1 - 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(['      YIF ' 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 ;
disp(blanks(1))


%% Full Sample

% Calculate Forecast Statistics 

NF = length(SPF(:,1)) ;                                 % Number of forecasts
NM = size(FE,3) ;                                       % Number of models
diary off 

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 ;                         % Yield spread
    elseif m == 3 ;         FCM = YIF ;                         % Yield inversion
    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 Output

diary on
disp('1975Q1 - 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)) ])
    disp(['      YIF ' 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 ;
disp(blanks(1))

toc 
