%% Matlab code to produce figures
% For paper
% 1 = plot of yield spread, NFCI and recession indictor (first final data)
% For appendix
% 2.1 = plot of yield spread, NFCI and recession indictor (advance data)
% 2.2 = rolling MAE/RMSE/LPS & FT Tests YSP versus SPF
% 2.3 = rolling MAE/RMSE/LPS & FT Tests NFCI versus SPF
% 2.4 = density plots

clear all
set(0,'defaultAxesFontName', 'Times New Roman')                 % Set font for axes
set(0,'defaultTextFontName', 'Times New Roman')                 % Set font for text
set(0,'defaultAxesFontSize', 16)                                % Set font size for axes
set(0,'defaultTextFontSize', 16)                                % Set font size for text

figflag = 2.4 ;

%% Plot Yield Spread, NFCI & R1 indicator -- First Final Data

if figflag == 1 ;

    load('YSRep.mat','qtr','R1FF','Spread') ;                   % Load data
    T1 = 63 ;                                                   % First period from which to forecast is 1968Q4
    TF = 283 ;                                                  % Final period from which to project is 2023Q4  
    qtrfc = qtr(T1:TF)' ;                                       % Dates for evaluation period
    
    % Create real time recession indicator
    R1indicator = [] ;
    for t = T1:TF ;
        R1indicator = [ R1indicator ; R1FF(t,t+1) ] ;           
    end ;
    clear R1FF t qtr
    
    load('YSRep1971.mat','qtr','NFCI')                          % Add NFCI data
    
    figure('position', [100, 100, 1000, 600])  
    tiledlayout(1,1,'TileSpacing', 'compact','Padding','compact')
    nexttile  
    plot(qtrfc,Spread(T1:TF),'b-','LineWidth',2) ;              % Plot yield spread
    hold
    plot(qtr,NFCI,'r-','LineWidth',2) ;                         % Add NFCI
    axis([1968 2024 -2 5])
    for t = 1:length(R1indicator) ;                             % Add bars for R1 indicator = 1
        if R1indicator(t) == 1 ;
            X = [ qtrfc(t)-0.125 qtrfc(t)+0.125 qtrfc(t)+0.125 qtrfc(t)-0.125] ;
            Y = [ -2 -2 5 5] ;
            patch(X,Y,[0.33 0.33 0.33],'Edgecolor','none','FaceAlpha',0.25) 
        end ;
    end ;
    set(gca,'TickLabelInterpreter','latex') 
    line([1968 2024],[0 0],'Color','k','LineStyle',':','LineWidth',1)
    legend('Yield Spread','NFCI','R1 Recession','location','southoutside','Orientation','horizontal','Interpreter','latex')

end ;

%% Plot Yield Spread, NFCI & R1 indicator -- Advance Data

if figflag == 2.1 ;
    
    load('YSRep.mat','qtr','R1AD','Spread') ;                   % Load data
    T1 = 63 ;                                                   % First period from which to forecast is 1968Q4
    TF = 283 ;                                                  % Final period from which to project is 2023Q4  
    qtrfc = qtr(T1:TF)' ;                                       % dates for evaluation period
    
    % Create real time recession indicator
    R1indicator = [] ;
    for t = T1:TF ;
        R1indicator = [ R1indicator ; R1AD(t,t+1) ] ;           
    end ;
    clear R1AD t qtr
    
    load('YSRep1971.mat','qtr','NFCI')

    figure('position', [100, 100, 1000, 600])  
    tiledlayout(1,1,'TileSpacing', 'compact','Padding','compact')
    nexttile  
    plot(qtrfc,Spread(T1:TF),'b-','LineWidth',2) ;              % Plot yield spread
    hold
    plot(qtr,NFCI,'r-','LineWidth',2) ;                         % Add NFCI
    axis([1968 2024 -2 5])
    for t = 1:length(R1indicator) ;                             % Add bars for R1 indicator = 1
        if R1indicator(t) == 1 ;
            X = [ qtrfc(t)-0.125 qtrfc(t)+0.125 qtrfc(t)+0.125 qtrfc(t)-0.125] ;
            Y = [ -2 -2 5 5] ;
            patch(X,Y,[0.33 0.33 0.33],'Edgecolor','none','FaceAlpha',0.25) 
        end ;
    end ;
    set(gca,'TickLabelInterpreter','latex') 
    line([1968 2024],[0 0],'Color','k','LineStyle',':','LineWidth',1)
    legend('Yield Spread','NFCI','R1 Recession','location','southoutside','Orientation','horizontal','Interpreter','latex')

end ;

%% Plot Metrics & Fluctuation Test for Moving Windows, h = 4 -- Narrow Replication

if figflag == 2.2 ;                                          
    
    % Load Forecast Metrics
    load('NarrowFS.mat','SPF_FE','YSF_FE','SPF_LS','YSF_LS') ;          % Load forecast errors and log-scores
    qtr = linspace(1968.75,2023.75,size(SPF_FE,1))' ;                   % Dates 1968Q4 - 2023Q4
    NF = size(SPF_FE,1) ;

    % 40 quarter windows

    W = 40 ;
 
    SPF_MAE = NaN(NF,1) ;                                               % Store recursive MAE for SPF here
    YSF_MAE = NaN(NF,1) ;                                               % Ditto yield spread
    SPF_RMSE = NaN(NF,1) ;                                              % Store recursive RMSE for SPF here
    YSF_RMSE = NaN(NF,1) ;                                              % Ditto yield spread
    SPF_LPS = NaN(NF,1) ;                                               % Store recursive LPS for SPF here
    YSF_LPS = NaN(NF,1) ;                                               % Ditto yield spread

    for t = W/2+1:NF-W/2+1 ;                                            % Do loop through time
        df = t - W/2 ;                                                  % First date of window
        dl = t + W/2 - 1 ;                                              % Last date of window
        SPF_MAE(t) = mean(abs(SPF_FE(df:dl,5))) ;                       % MAE SPF
        YSF_MAE(t) = mean(abs(YSF_FE(df:dl,5))) ;                       % MAE yield spread
        SPF_RMSE(t) = sqrt(mean(SPF_FE(df:dl,5).^2)) ;                  % RMSE SPF
        YSF_RMSE(t) = sqrt(mean(YSF_FE(df:dl,5).^2)) ;                  % RMSE yield spread
        SPF_LPS(t) = -mean(SPF_LS(df:dl,5)) ;                           % LPS SPF
        YSF_LPS(t) = -mean(YSF_LS(df:dl,5)) ;                           % LPS yield spread
    end ;
    
    relMAE40 = YSF_MAE./SPF_MAE ;                                       % Relative MAE
    relRMSE40 = YSF_RMSE./SPF_RMSE ;                                    % Relative RMSE
    relLPS40 = YSF_LPS./SPF_LPS ;                                       % Relative LPS
    
    % Fluctuation Test Statistics
    
    load('NarrowFT40.mat','FT','cv05','cv10')
    FT40_MAE = FT(:,5,1) ;                                              % FT stats - MAE
    FT40_MSE = FT(:,5,2) ;                                              % FT stats - MSE   
    FT40_LPS = FT(:,5,3) ;                                              % FT stats - LPS
    cv10m40 = cv10 ;

    clear W SPF_MAE YSF_MAE SPF_RMSE YSF_RMSE SPF_LPS YSF_LPS FT cv05 cv10
    
    % 60 quarter windows

    W = 60 ;
 
    SPF_MAE = NaN(NF,1) ;                                               % Store recursive MAE for SPF here
    YSF_MAE = NaN(NF,1) ;                                               % Ditto yield spread
    SPF_RMSE = NaN(NF,1) ;                                              % Store recursive RMSE for SPF here
    YSF_RMSE = NaN(NF,1) ;                                              % Ditto yield spread
    SPF_LPS = NaN(NF,1) ;                                               % Store recursive LPS for SPF here
    YSF_LPS = NaN(NF,1) ;                                               % Ditto yield spread

    for t = W/2+1:NF-W/2+1 ;                                            % Do loop through time
        df = t - W/2 ;                                                  % First date of window
        dl = t + W/2 - 1 ;                                              % Last date of window
        SPF_MAE(t) = mean(abs(SPF_FE(df:dl,5))) ;                       % MAE SPF
        YSF_MAE(t) = mean(abs(YSF_FE(df:dl,5))) ;                       % MAE yield spread
        SPF_RMSE(t) = sqrt(mean(SPF_FE(df:dl,5).^2)) ;                  % RMSE SPF
        YSF_RMSE(t) = sqrt(mean(YSF_FE(df:dl,5).^2)) ;                  % RMSE yield spread
        SPF_LPS(t) = -mean(SPF_LS(df:dl,5)) ;                           % LPS SPF
        YSF_LPS(t) = -mean(YSF_LS(df:dl,5)) ;                           % LPS yield spread
    end ;
    
    relMAE60 = YSF_MAE./SPF_MAE ;                                       % Relative MAE
    relRMSE60 = YSF_RMSE./SPF_RMSE ;                                    % Relative RMSE
    relLPS60 = YSF_LPS./SPF_LPS ;                                       % Relative LPS
    
    % Fluctuation Test Statistics
    
    load('NarrowFT60.mat','FT','cv05','cv10')
    FT60_MAE = FT(:,5,1) ;                                              % FT stats - MAE
    FT60_MSE = FT(:,5,2) ;                                              % FT stats - MSE   
    FT60_LPS = FT(:,5,3) ;                                              % FT stats - LPS
    cv10m60 = cv10 ;

    clear W SPF_MAE YSF_MAE SPF_RMSE YSF_RMSE SPF_LPS YSF_LPS FT cv05 cv10
   
    % Plot Figure

    figure('position', [100, 100, 1000, 1000])  
    tiledlayout(2,2,'Padding','compact')
    set(gca,'TickLabelInterpreter','latex') 

    nexttile
    h1 = plot(qtr,relMAE40,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,relRMSE40,'b','LineWidth',2) ;
    h3 = plot(qtr,relLPS40,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 1 1],'Color','black','LineStyle',':','LineWidth',2) ;
    xlim([1973 2019])
    ylim([0.5 1.5])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 40-quarter Rolling Window','Interpreter','latex')
    ylabel('Relative Forecast Metric','Interpreter','latex')
    legend([ h1 h2 h3],'MAE','RMSE','LPS','Interpreter','latex')
    title('(a) Rolling Window Metrics, $m = 40$','Interpreter','latex')
    
    nexttile
    h1 = plot(qtr,FT40_MAE,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,FT40_MSE,'b','LineWidth',2) ;
    h3 = plot(qtr,FT40_LPS,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 0 0],'Color','black','LineStyle',':','LineWidth',2) ;
    h5 = line([ qtr(1) qtr(NF)],[ cv10m40 cv10m40],'Color','black','LineStyle','--','LineWidth',2) ;
    h6 = line([ qtr(1) qtr(NF)],[ -cv10m40 -cv10m40],'Color','black','LineStyle','--','LineWidth',2) ;
    xlim([1973 2019])
    ylim([-6 6])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 40-quarter Rolling Window','Interpreter','latex')
    ylabel('Test Statistic','Interpreter','latex')
    legend('MAE','MSE','LPS','Interpreter','latex')
    title('(b) Rolling Window Fluctuation Test, $m = 40$','Interpreter','latex')

    nexttile
    h1 = plot(qtr,relMAE60,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,relRMSE60,'b','LineWidth',2) ;
    h3 = plot(qtr,relLPS60,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 1 1],'Color','black','LineStyle',':','LineWidth',2) ;
    xlim([1973 2019])
    ylim([0.5 1.5])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 60-quarter Rolling Window','Interpreter','latex')
    ylabel('Relative Forecast Metric','Interpreter','latex')
    legend([ h1 h2 h3],'MAE','RMSE','LPS','Interpreter','latex')
    title('(c) Rolling Window Metrics, $m = 60$','Interpreter','latex')

    nexttile
    h1 = plot(qtr,FT60_MAE,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,FT60_MSE,'b','LineWidth',2) ;
    h3 = plot(qtr,FT60_LPS,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 0 0],'Color','black','LineStyle',':','LineWidth',2) ;
    h5 = line([ qtr(1) qtr(NF)],[ cv10m60 cv10m60],'Color','black','LineStyle','--','LineWidth',2) ;
    h6 = line([ qtr(1) qtr(NF)],[ -cv10m60 -cv10m60],'Color','black','LineStyle','--','LineWidth',2) ;
    xlim([1973 2019])
    ylim([-6 6])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 60-quarter Rolling Window','Interpreter','latex')
    ylabel('Test Statistic','Interpreter','latex')
    legend('MAE','MSE','LPS','Interpreter','latex')
    title('(d) Rolling Window Fluctuation Test, $m = 60$','Interpreter','latex')
    
end ;

%% Plot Metrics & Fluctuation Test for Moving Windows, h = 4 --- Broad Replication

if figflag == 2.3 ;                                          
    
    % Load Forecast Metrics
    load('BroadFS.mat','SPF_FE','FCI_FE','SPF_LS','FCI_LS') ;           % Load forecast errors and log-scores
    qtr = linspace(1988,2021.25,size(SPF_FE,1))' ;                      % Dates 1988Q1 - 2021Q2
    NF = size(SPF_FE,1) ;

    % 40 quarter windows

    W = 40 ;
 
    SPF_MAE = NaN(NF,1) ;                                               % Store recursive MAE for SPF here
    FCF_MAE = NaN(NF,1) ;                                               % Ditto NFCI
    SPF_RMSE = NaN(NF,1) ;                                              % Store recursive RMSE for SPF here
    FCF_RMSE = NaN(NF,1) ;                                              % Ditto NFCI
    SPF_LPS = NaN(NF,1) ;                                               % Store recursive LPS for SPF here
    FCF_LPS = NaN(NF,1) ;                                               % Ditto NFCI

    for t = W/2+1:NF-W/2+1 ;                                            % Do loop through time
        df = t - W/2 ;                                                  % First date of window
        dl = t + W/2 - 1 ;                                              % Last date of window
        SPF_MAE(t) = mean(abs(SPF_FE(df:dl,5))) ;                       % RMSE SPF
        FCF_MAE(t) = mean(abs(FCI_FE(df:dl,5))) ;                       % RMSE NFCI
        SPF_RMSE(t) = sqrt(mean(SPF_FE(df:dl,5).^2)) ;                  % RMSE SPF
        FCF_RMSE(t) = sqrt(mean(FCI_FE(df:dl,5).^2)) ;                  % RMSE NFCI
        SPF_LPS(t) = -mean(SPF_LS(df:dl,5)) ;                           % RMSE SPF
        FCF_LPS(t) = -mean(FCI_LS(df:dl,5)) ;                           % RMSE NFCI
    end ;
    
    relMAE40 = FCF_MAE./SPF_MAE ;                                       % Relative MAE
    relRMSE40 = FCF_RMSE./SPF_RMSE ;                                    % Relative RMSE
    relLPS40 = FCF_LPS./SPF_LPS ;                                       % Relative LPS
    
    % Fluctuation Test
    
    load('BroadFT40.mat','FT','cv05','cv10')
    FT40_MAE = FT(:,5,1) ;                                              % FT stats - MAE
    FT40_MSE = FT(:,5,2) ;                                              % FT stats - MSE   
    FT40_LPS = FT(:,5,3) ;                                              % FT stats - LPS
    cv10m40 = cv10 ;

    clear W SPF_MAE FCF_MAE SPF_RMSE FCF_RMSE SPF_LPS FCF_LPS FT cv05 cv10
    
     % 60 quarter windows

    W = 60 ;
 
    SPF_MAE = NaN(NF,1) ;                                               % Store recursive MAE for SPF here
    FCF_MAE = NaN(NF,1) ;                                               % Ditto NFCI
    SPF_RMSE = NaN(NF,1) ;                                              % Store recursive RMSE for SPF here
    FCF_RMSE = NaN(NF,1) ;                                              % Ditto NFCI
    SPF_LPS = NaN(NF,1) ;                                               % Store recursive LPS for SPF here
    FCF_LPS = NaN(NF,1) ;                                               % Ditto NFCI

    for t = W/2+1:NF-W/2 ;                                              % Do loop through time
        df = t - W/2 ;                                                  % First date of window
        dl = t + W/2 - 1 ;                                              % Last date of window
        SPF_MAE(t) = mean(abs(SPF_FE(df:dl,5))) ;                       % RMSE SPF
        FCF_MAE(t) = mean(abs(FCI_FE(df:dl,5))) ;                       % RMSE NFCI
        SPF_RMSE(t) = sqrt(mean(SPF_FE(df:dl,5).^2)) ;                  % RMSE SPF
        FCF_RMSE(t) = sqrt(mean(FCI_FE(df:dl,5).^2)) ;                  % RMSE NFCI
        SPF_LPS(t) = -mean(SPF_LS(df:dl,5)) ;                           % RMSE SPF
        FCF_LPS(t) = -mean(FCI_LS(df:dl,5)) ;                           % RMSE NFCI
    end ;
    
    relMAE60 = FCF_MAE./SPF_MAE ;                                       % Relative MAE
    relRMSE60 = FCF_RMSE./SPF_RMSE ;                                    % Relative MSE
    relLPS60 = FCF_LPS./SPF_LPS ;                                       % Relative LPS
    
    load('BroadFT60.mat','FT','cv05','cv10')
    FT60_MAE = FT(:,5,1) ;                                              % FT stats - MAE
    FT60_MSE = FT(:,5,2) ;                                              % FT stats - MSE  
    FT60_LPS = FT(:,5,3) ;                                              % FT stats - LPS
    cv10m60 = cv10 ;

    clear W SPF_MAE FCF_MAE SPF_RMSE FCF_RMSE SPF_LPS FCF_LPS FT cv05 cv10
   
    % Plot Figure

    figure('position', [100, 100, 1000, 1000])  
    tiledlayout(2,2,'Padding','compact')
    set(gca,'TickLabelInterpreter','latex') 

    nexttile
    h1 = plot(qtr,relMAE40,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,relRMSE40,'b','LineWidth',2) ;
    h3 = plot(qtr,relLPS40,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 1 1],'Color','black','LineStyle',':','LineWidth',2) ;
    xlim([1992 2016])
    ylim([0.5 1.5])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 40-quarter Rolling Window','Interpreter','latex')
    ylabel('Relative Forecast Metric','Interpreter','latex')
    legend([ h1 h2 h3],'MAE','RMSE','LPS','Interpreter','latex')
    title('(a) Rolling Window Metrics, $m = 40$','Interpreter','latex')
    
    nexttile
    h1 = plot(qtr,FT40_MAE,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,FT40_MSE,'b','LineWidth',2) ;
    h3 = plot(qtr,FT40_LPS,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 0 0],'Color','black','LineStyle',':','LineWidth',2) ;
    h5 = line([ qtr(1) qtr(NF)],[ cv10m40 cv10m40],'Color','black','LineStyle','--','LineWidth',2) ;
    h6 = line([ qtr(1) qtr(NF)],[ -cv10m40 -cv10m40],'Color','black','LineStyle','--','LineWidth',2) ;
    xlim([1992 2016])
    ylim([-6 6])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 40-quarter Rolling Window','Interpreter','latex')
    ylabel('Test Statistic','Interpreter','latex')
    legend('MAE','MSE','LPS','Interpreter','latex')
    title('(b) Rolling Window Fluctuation Test, $m = 40$','Interpreter','latex')

    nexttile
    h1 = plot(qtr,relMAE60,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,relRMSE60,'b','LineWidth',2) ;
    h3 = plot(qtr,relLPS60,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 1 1],'Color','black','LineStyle',':','LineWidth',2) ;
    xlim([1992 2016])
    ylim([0.5 1.5])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 60-quarter Rolling Window','Interpreter','latex')
    ylabel('Relative Forecast Metric','Interpreter','latex')
    legend([ h1 h2 h3],'MAE','RMSE','LPS','Interpreter','latex')
    title('(c) Rolling Window Metrics, $m = 60$','Interpreter','latex')

    nexttile
    h1 = plot(qtr,FT60_MAE,'r','LineWidth',2) ;
    hold 
    h2 = plot(qtr,FT60_MSE,'b','LineWidth',2) ;
    h3 = plot(qtr,FT60_LPS,'g','LineWidth',2) ;
    h4 = line([ qtr(1) qtr(NF)],[ 0 0],'Color','black','LineStyle',':','LineWidth',2) ;
    h5 = line([ qtr(1) qtr(NF)],[ cv10m60 cv10m60],'Color','black','LineStyle','--','LineWidth',2) ;
    h6 = line([ qtr(1) qtr(NF)],[ -cv10m60 -cv10m60],'Color','black','LineStyle','--','LineWidth',2) ;
    xlim([1992 2016])
    ylim([-6 6])
    set(gca,'TickLabelInterpreter','latex') 
    xlabel('Midpoint of 60-quarter Rolling Window','Interpreter','latex')
    ylabel('Test Statistic','Interpreter','latex')
    legend('MAE','MSE','LPS','Interpreter','latex')
    title('(d) Rolling Window Fluctuation Test, $m = 60$','Interpreter','latex')    

end ;

%% Density Plots 1971:1 - 2023:4

if figflag == 2.4 ;
    
    load('YSRep1971.mat','Spread','NFCI')

    % Spread 
    SpreadPts = linspace(-3,6,1000) ;
    SpreadEmp = ksdensity(Spread,SpreadPts) ;
    SpreadNorm = mean(Spread) + std(Spread)*randn(1000000,1) ;
    SpreadNor = ksdensity(SpreadNorm,SpreadPts) ;

    % NFCI
    NFCIPts = linspace(-4,6,1000) ;
    NFCIEmp = ksdensity(NFCI,NFCIPts) ;
    NFCINorm = mean(NFCI) + std(NFCI)*randn(1000000,1) ;
    NFCINor = ksdensity(NFCINorm,NFCIPts) ;

    figure('position', [100, 100, 1100, 600])  
    tiledlayout(1,2,'TileSpacing', 'compact','Padding','compact')

    nexttile
    plot(SpreadPts,SpreadEmp,'LineWidth',2,'Color','r') ;
    hold 
    plot(SpreadPts,SpreadNor,'LineWidth',2,'Color','b') ;
    set(gca,'TickLabelInterpreter','latex') 
    title('(a) Yield Spread','Interpreter','latex')
    legend('Empirical','Normal','Interpreter','latex')

    nexttile
    plot(NFCIPts,NFCIEmp,'LineWidth',2,'Color','r') ;
    hold 
    plot(NFCIPts,NFCINor,'LineWidth',2,'Color','b') ;
    set(gca,'TickLabelInterpreter','latex') 
    title('(b) NFCI','Interpreter','latex')
    legend('Empirical','Normal','Interpreter','latex')

end ;