%% Produce Fluctuation Test Statistics (Gianonne & Rossi, JAE 2010)
% This version is for the broader replication and first final data only - RTNFCI versus SPF

clear all
tic 
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 = 1 ;                                                   % Set to 1 to produce plots
m = 60 ;                                                        % Set window length in quarters 

%% 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('YSRep1973.mat','R1FF','SPF') ;                    % Data for out-turn and 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
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

load('BroadFF.mat','OosPR') ;
qtr = OosPR(:,1) ;                                      % Quarters
FCF = OosPR(:,7:11) ;                                   % RTNFCI forecasts

NF = size(SPF,1) ;                                      % Number of forecasts 
P = NF + 1 ;                                            % NF + 1
Q = round(P^0.25) ;                                     % Bandwidth for HAC estimator 

%% Calculate forecast errors & loss differentials
% Forecast 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

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

% Forecast Errors and Log Scores - 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 ;

% Forecast Errors and Log Scores - NFCI
for h = 0:4 ;
    FE(1:NF-h,h+1,2) = FCF(1:NF-h,h+1) - OutTurn(h+1:NF) ;                         
    LS(1:NF-h,h+1,2) = (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 ;

% Loss Diffferentials
LD = NaN(NF,5,3) ;                                                  % Store loss differentials here
for h = 0:4 ;                                                       % Do loop through horizon
    LD(:,h+1,1) = abs(FE(:,h+1,2)) - abs(FE(:,h+1,1)) ;             % Loss differential for MAE
    LD(:,h+1,2) = (FE(:,h+1,2)).^2 - (FE(:,h+1,1)).^2 ;             % Ditto MSE, next line is for LS - set as SPF-YSP so negative indicates YSP performs better
    LD(1:NF-h,h+1,3) = (1-OutTurn(h+1:NF)).*(log(1-SPF(1:NF-h,h+1))-log(1-FCF(1:NF-h,h+1))) + OutTurn(h+1:NF).*(log(SPF(1:NF-h,h+1))-log(FCF(1:NF-h,h+1))) ;
end ;

%% HAC variance estimator 
HACSIG = NaN(5,3) ;                                                 % Store HAC variance estimator here
for h = 0:4 ;                                                       % Do loop through horizon
    for l = 1:3 ;                                                   % Do loop through MAE / MSE
        LDhl = LD(1:NF-h,h+1,l) ;                                   % Loss differential at horizon h
        NFhl = size(LDhl,1) ;                                       % Number obs
        sigsq = 0 ;                                                 % Initialize HAC estimator
        for i = -Q+1:Q-1 ;                                          % Do loop through leads/lags
            sigsq = sigsq + (1-abs(i/Q))*(1/P)*sum(LDhl(Q:NFhl-Q+1).*LDhl(Q-i:NFhl-Q+1-i)) ;
        end ;                                                       % Previous line adds component for lead/lag i
    HACSIG(h+1,l) = sqrt(sigsq) ;                                   % Store HAC estimator
    clear LDhl NFhl i sigsq
    end ;
end ;
clear h l Q

%% Calculate Fluctuation Test Statistics

FT = NaN(NF,5,3) ;                                                      % Store fluctuation test stats here, rows = dates, cols = horizon, page = MAE, MSE, LS

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

    for t = m/2+1:NF-m/2-h+1 ;                                          % Do loop through time        

        df = t - m/2 ;                                                  % First date of window
        dl = t + m/2 - 1 ;                                              % Last date of window   
        
        if dl == NF ;                                                   % If last quarter of window is last forecast origin....                                             
             if h == 0 ;                                                % ....only calculate statistics for h = 0

                sum1ht = sum(abs(FE(df:dl,h+1,1))) ;                    % Sum MAE for SPF
                sum2ht = sum(abs(FE(df:dl,h+1,2))) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,1)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MAE    
                FT(t,h+1,1) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht       
        
                sum1ht = sum((FE(df:dl,h+1,1)).^2) ;                    % Sum MAE for SPF
                sum2ht = sum((FE(df:dl,h+1,2)).^2) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,2)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MSE
                FT(t,h+1,2) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht

                diffht = (1-OutTurn(df+h:dl+h)).*(log(1-SPF(df:dl,h+1))-log(1-FCF(df:dl,h+1))) + OutTurn(df+h:dl+h).*(log(SPF(df:dl,h+1))-log(FCF(df:dl,h+1))) ;
                Fht = HACSIG(h+1,3)^(-1)*m^(-0.5)*sum(diffht) ;         % Test stat for LS
                FT(t,h+1,3) = Fht ;                                     % Store test stat for horizon h at time t
                clear diffht Fht
            
             else ;
                FT(t,h+1,:) = NaN ;
            end ;
        
        elseif dl == NF-1 ;                                             % If last quarter of window is penultimate forecast origin....                                                                 
            if h <= 1 ;                                                 % ....only calculate statistics for h = 0 & 1
                
                sum1ht = sum(abs(FE(df:dl,h+1,1))) ;                    % Sum MAE for SPF
                sum2ht = sum(abs(FE(df:dl,h+1,2))) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,1)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MAE    
                FT(t,h+1,1) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht       
        
                sum1ht = sum((FE(df:dl,h+1,1)).^2) ;                    % Sum MAE for SPF
                sum2ht = sum((FE(df:dl,h+1,2)).^2) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,2)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MSE
                FT(t,h+1,2) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht

                diffht = (1-OutTurn(df+h:dl+h)).*(log(1-SPF(df:dl,h+1))-log(1-FCF(df:dl,h+1))) + OutTurn(df+h:dl+h).*(log(SPF(df:dl,h+1))-log(FCF(df:dl,h+1))) ;
                Fht = HACSIG(h+1,3)^(-1)*m^(-0.5)*sum(diffht) ;         % Test stat for LS
                FT(t,h+1,3) = Fht ;                                     % Store test stat for horizon h at time t
                clear diffht Fht df dl
            
            else ;
                FT(t,h+1,:) = NaN ;
            end ;

        elseif dl == NF-2 ;                                             % If last quarter of window is two from last forecast origin....
            if h <= 2 ;                                                 % ....only calculate statistics for h = 0, 1 & 2
                
                sum1ht = sum(abs(FE(df:dl,h+1,1))) ;                    % Sum MAE for SPF
                sum2ht = sum(abs(FE(df:dl,h+1,2))) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,1)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MAE    
                FT(t,h+1,1) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht       
        
                sum1ht = sum((FE(df:dl,h+1,1)).^2) ;                    % Sum MAE for SPF
                sum2ht = sum((FE(df:dl,h+1,2)).^2) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,2)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MSE
                FT(t,h+1,2) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht

                diffht = (1-OutTurn(df+h:dl+h)).*(log(1-SPF(df:dl,h+1))-log(1-FCF(df:dl,h+1))) + OutTurn(df+h:dl+h).*(log(SPF(df:dl,h+1))-log(FCF(df:dl,h+1))) ;
                Fht = HACSIG(h+1,3)^(-1)*m^(-0.5)*sum(diffht) ;         % Test stat for LS
                FT(t,h+1,3) = Fht ;                                     % Store test stat for horizon h at time t
                clear diffht Fht df dl
            
            else ;
                FT(t,h+1,:) = NaN ;
            end ;
        
        elseif dl == NF-3 ;                                             % If last quarter of window is three from last forecast origin....
            if h <= 3 ;                                                 % ....only calculate statistics for h = 0, 1, 2 & 3
                
                sum1ht = sum(abs(FE(df:dl,h+1,1))) ;                    % Sum MAE for SPF
                sum2ht = sum(abs(FE(df:dl,h+1,2))) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,1)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MAE    
                FT(t,h+1,1) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht       
        
                sum1ht = sum((FE(df:dl,h+1,1)).^2) ;                    % Sum MAE for SPF
                sum2ht = sum((FE(df:dl,h+1,2)).^2) ;                    % Ditto for YSP
                Fht = HACSIG(h+1,2)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;     % Test stat for MSE
                FT(t,h+1,2) = Fht ;                                     % Store test stat for horizon h at time t
                clear sum1ht sum2ht Fht

                diffht = (1-OutTurn(df+h:dl+h)).*(log(1-SPF(df:dl,h+1))-log(1-FCF(df:dl,h+1))) + OutTurn(df+h:dl+h).*(log(SPF(df:dl,h+1))-log(FCF(df:dl,h+1))) ;
                Fht = HACSIG(h+1,3)^(-1)*m^(-0.5)*sum(diffht) ;         % Test stat for LS
                FT(t,h+1,3) = Fht ;                                     % Store test stat for horizon h at time t
                clear diffht Fht df dl
            
            else ;
                FT(t,h+1,:) = NaN ;
            end ;
        
        elseif dl <= NF-4 ;                                             % If last quarter of window is four or more from last forecast origin calculate statistics for all h
                
            sum1ht = sum(abs(FE(df:dl,h+1,1))) ;                        % Sum MAE for SPF
            sum2ht = sum(abs(FE(df:dl,h+1,2))) ;                        % Ditto for YSP
            Fht = HACSIG(h+1,1)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;         % Test stat for MAE    
            FT(t,h+1,1) = Fht ;                                         % Store test stat for horizon h at time t
            clear sum1ht sum2ht Fht       
        
            sum1ht = sum((FE(df:dl,h+1,1)).^2) ;                        % Sum MAE for SPF
            sum2ht = sum((FE(df:dl,h+1,2)).^2) ;                        % Ditto for YSP
            Fht = HACSIG(h+1,2)^(-1)*m^(-0.5)*(sum2ht-sum1ht) ;         % Test stat for MSE
            FT(t,h+1,2) = Fht ;                                         % Store test stat for horizon h at time t
            clear sum1ht sum2ht Fht

            diffht = (1-OutTurn(df+h:dl+h)).*(log(1-SPF(df:dl,h+1))-log(1-FCF(df:dl,h+1))) + OutTurn(df+h:dl+h).*(log(SPF(df:dl,h+1))-log(FCF(df:dl,h+1))) ;
            Fht = HACSIG(h+1,3)^(-1)*m^(-0.5)*sum(diffht) ;             % Test stat for LS
            FT(t,h+1,3) = Fht ;                                         % Store test stat for horizon h at time t
            clear diffht Fht

        end ;

    end ;

    clear df dl t

end ;

clear FE LS LD h

%% Critical Values

mu = round(m/P,1) ;

    if mu == 0.1 ;      cv05 = 3.393 ;  cv10 = 3.170 ;
elseif mu == 0.2 ;      cv05 = 3.179 ;  cv10 = 2.948 ;
elseif mu == 0.3 ;      cv05 = 3.012 ;  cv10 = 2.766 ;
elseif mu == 0.4 ;      cv05 = 2.890 ;  cv10 = 2.626 ;
elseif mu == 0.5 ;      cv05 = 2.779 ;  cv10 = 2.500 ;
elseif mu == 0.6 ;      cv05 = 2.634 ;  cv10 = 2.356 ;
elseif mu == 0.7 ;      cv05 = 2.560 ;  cv10 = 2.248 ;
elseif mu == 0.8 ;      cv05 = 2.433 ;  cv10 = 2.080 ;
elseif mu == 0.9 ;      cv05 = 2.248 ;  cv10 = 1.975 ;
end ;

clear mu P 

%% Save FT statistics and critical values

    if m == 40 ;
    save('BroadFT40.mat', 'FT', 'cv05', 'cv10') ;
elseif m == 60 ;
    save('BroadFT60.mat', 'FT', 'cv05', 'cv10') ;
    end ;

%% Plots

if figflag == 1 ;
    
    figure('position', [100, 100, 950, 1200])  

    subplot(5,3,1)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,1,1),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MAE and h = 0')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,2)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,1,2),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MSE and h = 0')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,3)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,1,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('LS and h = 0')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,4)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,2,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MAE and h = 1')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,5)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,2,2),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MSE and h = 1')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,6)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,2,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('LS and h = 1')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,7)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,2,1),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MAE and h = 2')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,8)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,3,2),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MSE and h = 2')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,9)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,2,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('LS and h = 2')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,10)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,4,1),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MAE and h = 3')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,11)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,4,2),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MSE and h = 3')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,12)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,4,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('LS and h = 3')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,13)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,5,1),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MAE and h = 4')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,14)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,5,2),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('MSE and h = 4')
    ylim([-5 5])
    xlim([1992 2017])

    subplot(5,3,15)
    plot(qtr(m/2+1:NF-m/2),FT(m/2+1:NF-m/2,5,3),'b-','linewidth',1.5) ;
    hold 
    line([ qtr(m/2+1) qtr(NF-m/2)],[cv10 cv10],'Color','red','LineStyle','--')
    line([ qtr(m/2+1) qtr(NF-m/2)],[-cv10 -cv10],'Color','red','LineStyle','--')
    title('LS and h = 4')
    ylim([-5 5])
    xlim([1992 2017])

end ;

toc


