clear
close all
clc
% =========================================================================
% DESCRIPTION 
% This script loads in the transformed y variables to be predicted - CPI 
% and IP - and compares the performance of 1) factor estimated 
% with Hamilton & Xi's approach, 2) factor estimated with 	
% McCraken & Ng's approach, and 3) CFNAI. Comparisons are based on the 
% ratio of MSE from their out-of-sample forecasts relative to MSE from 
% forecasts with only autoregressive lags. 

% -------------------------------------------------------------------------
% INPUT FILES 
%       mnfactor_'year'.mat - MAT-file containing FRED factors over using  
%       McCracken & Ng's method. Sample starts from 1960:03 and ends by the 
%       year indicated by the file name. This data is generated by code 
%       "main_replicate.m".
%
%       hxfactor_'year'.mat - MAT-file containing cyclical factor using
%       Hamilton & Xi's approach. By default sample starts from 1963:01 and  
%       ends by the year indicated by the file name. This data is generated
%       by code "application_FREDMD.m".
%
%       yvar_'year'.mat - MAT-file containing y variables to be predicted. 
%       This data is generated by code "generate_yvar.m".        
%
%       cfnai_2023.mat - CFNAI data 
%
%       hxfactor_real_2023.mat - All variables detrended using regression 
%       coefficients estimated from data only through end of 2014

% OUTPUT FILES
%       MSE_MN_vs_CFNAI_vs_HX.xlsx - Excel file 
%       containing mean squared errors of forecasts with autoregressive
%       lags, and relative ratio of MSE from froecasts with McCraken & 
%       Ng's factor, CFNAI, and Hamilton & Xi's factor compared to  
%       autoregression.
%
%       MSE_MN.xlsx - Excel file containing mean squared errors of  
%       autoregression and relative ratio  of MSE from forecasts with 
%       McCraken & Ng's factor
% -------------------------------------------------------------------------
% BREAKDOWN OF THE SCRIPT
% 
% Part 1: Choose factors used for forecasting exercise
%
% Part 2: Set up
%
% Part 3: Load in data
%
% Part 4- Aggregate y variables to be forecasted 
%
% Part 5: Extract the sample period used for forecasting exercise.
%
% Part 6: Perform out-of-sample forcasts and calculate mean squared errors.
%          
%
% -------------------------------------------------------------------------
% AUXILIARY FUNCTIONS
%
%   mse() - performs a recursive out-of-sample forecasting exercise. 

% =========================================================================
% Part 1: Choose factors used for forecasting exercise

% =========================================================================
% Part 2: Set up

% Key dates and load in data
    
    % Date of the first factor used for forecasting 
    date_start_factor = 1967+2/12;  % March 1967
    
    % Out-of-sample prediction periods - each row contains a starting date
    % and an ending date of forecast. 
    dates_forecast = [1970, 1996+11/12; % Jan 1970 - Dec 1996
                  1997, 2014+11/12;     % Jan 1997 - Dec 2014
                  2015, 2024+8/12];    % Jan 2015 - Aug 2024
              
    % Number of lags used to forecast IP, Emp, and CPI
    
    numlag_y = [6,6,6; % number of autoregressive lags, CPI, h=1, 6, 12 
                6,6,6];  % number of autoregressive lags, IP, h=1, 6, 12 
    
    numlag_f = [6,6,6; % number of factor lags, CPI, h=1, 6, 12 
                6,6,6]; % number of factor lags, IP, h=1, 6, 12 
                 
    % Models 
models_cell = {'AR' , 'CF', 'MN1','HX1', 'MN2', 'HX2', 'HX1*' , 'HX2*'};

% =========================================================================
% Part 3: load in data 

R = size(dates_forecast,1); % number of forecasting exercise to perform

load hxfactor_real_2024;

for r = 1:R
    
    version = dates_forecast(r,2);
    
    if version == 1996+11/12
        input_yvar = 'yvar_2015.mat';
        input_mn = 'mnfactor_all_1997.mat';
        input_hx = 'hxfactor_1997.mat';
    elseif version == 2014+11/12
        input_yvar = 'yvar_2015.mat';
        input_mn = 'mnfactor_all_2015.mat';
        input_hx = 'hxfactor_2015.mat';
    elseif version == 2023+5/12
        input_yvar = 'yvar_2023.mat';
        input_mn = 'mnfactor_all_2023.mat';
        input_hx = 'hxfactor_2023.mat';
    elseif version == 2024+8/12
        input_yvar = 'yvar_2024.mat';
        input_mn = 'mnfactor_all_2024.mat';
        input_hx = 'hxfactor_2024.mat';

    end
    input_cfnai = 'cfnai_2024.mat';
    
    load(input_yvar) % y variable to be forecasted 

    load(input_mn) % FRED factors using McKracken & Ng's method
        
        load(input_cfnai) % CFNAI over reduced sample 

        load(input_hx) % HX factor using Hamilton & Xi's method

% ========================================================================
% Part 4: Aggregate y variables to be forecasted - CPI and IP

    % Variable names

    yname=['CPI       ';'IP        '];

    yname_cell=cellstr(yname); % cell

    % Horizons of 1, 6, and 12 months
    horizon=[1,6,12]; % array
    h_string=['1 ';'6 ';'12']; % string
    h_cell=cellstr(h_string); % cell


% -------------------------------------------------------------------------
% Combine y variables

    % Aggregate by horizon
    y_1 = [CPI_1,IP_1]; % CPI, IP h=1
    y_6 = [CPI_6,IP_6]; % CPI, IP h=6
    y_12 = [CPI_12,IP_12]; % CPI, IP h=12

    % Combine into one matrix
    Y_All(:,:,1)=y_1;
    Y_All(:,:,2)=y_6;
    Y_All(:,:,3)=y_12;

% =========================================================================
% Part 5: Sample used for forecasting in y variable and factor

% Determine rows corresponding to starting date of factor and ending date 
% of forecast, and keep data in between

tol = 1.e-10;   % this is tolerance to allow for floating-point precision

    % starting date and ending
    date_start_forecast = dates_forecast(r,1);
    date_end_forecast = dates_forecast(r,2);

    % row of date_start_factor in y variables
    [~,start_y] = ismembertol(date_start_factor,dates_num_yvar,tol); 

    % row of date_forecast_end in yvar
    [~,end_y] = ismembertol(date_end_forecast,dates_num_yvar,tol); 

    % keep y variable in between 
    Y = Y_All(start_y:end_y,:,:);
    
    % Uniform date vector 
    dates_num = dates_num_yvar(start_y: end_y);

% -------------------------------------------------------------------------
% Aggregate factors

 % MNfactors
    % row of date_start_factor in FRED factors
    [~,start_mn] = ismembertol(date_start_factor,dates_num_mn,tol);

    % row of date_end_forecast in FRED factors
    [~,end_mn] = ismembertol(date_end_forecast,dates_num_mn,tol); 

% HX factors
     % row of date_fstart in OLS filtered factors
     [~,start_hx] = ismembertol(date_start_factor,dates_num_hx,tol); 

     % row of date_pend in OLS filtered factors
     [~,end_hx] = ismembertol(date_end_forecast,dates_num_hx,tol);

 % CFNAI
     % row of date_start_factor in CFNAI
     [~,start_cfnai] = ismembertol(date_start_factor,dates_num_cfnai,tol); 

     % row of date_end_forecast in CFNAI 
     [~,end_cfnai] = ismembertol(date_end_forecast,dates_num_cfnai,tol); 

     cfnai = cfnai(start_cfnai:end_cfnai);

% add real-time HX factors if r == R


% Summarize all factors in a matrix
    factors_All = cfnai;   
    factors_All = [factors_All fhat_mn(start_mn: end_mn,1)];
    factors_All = [factors_All fhat_hx(start_hx: end_hx,1)];
    factors_All = [factors_All fhat_mn(start_mn: end_mn,2)];
    factors_All = [factors_All fhat_hx(start_hx: end_hx,2)];
    if r < R
        factors_All = [factors_All fhat_hx(start_hx: end_hx,1)];
        factors_All = [factors_All fhat_hx(start_hx: end_hx,2)];  
    elseif r == R 
        factors_All = [factors_All fhat_real_hx(start_hx: end_hx,1:2)];
    end  
    
    [T,f_n] = size(factors_All); % dimensions of factor matrix  

% =========================================================================
% Part 6: Out-of-sample prediction and MSE Calculation

% Determine the row corresponding the date_start_forecast (determine the 
% starting point of out-of-sample prediction in Y)
    [~,pstart] = ismembertol(date_start_forecast,dates_num,tol); 

% -------------------------------------------------------------------------
% Estimate MSE


    for i = 1:2 
        for j = 1:3
            h = horizon(j); % prediction horizon
            nlag_y = numlag_y(i,j); % number of autoregressive lags
            nlag_f = numlag_f(i,j); % number of factor lags
            
            % initiate regressor matrices
            X_base = []; % regressor matrix with only autoregressive lags
            X_f = []; % regressor matrix with autoregressive and factor lags
            
            % form autoregressive lags
            if nlag_y > 0
                for k = 1:nlag_y
                    ylag(:,k) = Y_All(start_y+nlag_f-k: end_y-h-k+1,i,1);
                end
                X_base = [X_base ylag];
                for k = 1:f_n
                    X_f(:,1:nlag_y,k) = ylag;
                end
            else
            end
            
            % form factor lags
            if nlag_f > 0
                for k = 1:f_n
                    for s = 1:nlag_f
                        flag(:,s) = factors_All(s: T-h-nlag_f+s,k);  
                    end
                    X_f(1:T-h-nlag_f+1,nlag_y+1:nlag_y+nlag_f,k) = flag;
                    clear flag
                end
            else
            end
            
            % add an intercept
            cons = ones(size(X_f,1),1);
            X_base = [X_base cons];
            for k = 1:f_n
                X_f(1:T-h-nlag_f+1,nlag_y+nlag_f+1,k) = cons;
            end
            
            % sample of y variable used for estimation 
            Y1 = Y(:,i,j);

            % estimate mse 
            mse_base = mse(Y1,X_base,h,pstart,nlag_f);
            for k = 1:f_n
                X = X_f(:,:,k);
                mse_factor(1,k) = mse(Y1,X,h,pstart,nlag_f);
                ratio{1,k} = round(mse_factor(1,k)/mse_base,2);
            end
            
            % save results to cell_MSE
            cell_MSE{3,1} = strcat(string(floor(date_start_forecast)),'-',string(ceil(date_end_forecast)));
            cell_MSE{2+j,2}= strcat('h=',h_cell{j});
            cell_MSE(2,(f_n+2)*(i-1)+3:(f_n+2)*i+1) = models_cell;
            cell_MSE{1,(f_n+2)*(i-1)+4} = yname_cell{i};
            cell_MSE{2+j,(f_n+2)*(i-1)+3} = round(mse_base,2);
            cell_MSE(2+j,(f_n+2)*(i-1)+4:(f_n+2)*i+1) = ratio;
            
            
            clear ylag X_base X_f flag_fred flag_cfnai flag_ols Y1 Y2 
            
        end   
    end
    
    % print results to Excel sheet  
        output_excel = 'MSE_MN_vs_CFNAI_vs_HX_all.xlsx'; % excel name
     

    if r == 1
        writecell(cell_MSE, output_excel, 'Sheet',1,'Range','A1')
    elseif r==2
        MSE = cell_MSE(3:5,:);
        writecell(MSE, output_excel, 'Sheet',1,'Range','A7')
    elseif r==3
        MSE = cell_MSE(3:5,:);
        writecell(MSE, output_excel, 'Sheet',1,'Range','A11')
    elseif r==4
        MSE = cell_MSE(3:5,:);
        writecell(MSE, output_excel, 'Sheet',1,'Range','A15')
    end 

    clear Y_All factor_All cell_MSE
end
    
%% Subfunction for Out-of-Sample Inference
% This function performs a recursive out-of-sample forecasting exercise. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
    % Y:            regressand
    % X:             
    % h:            forecasting horizon
    % pstart:       index of the starting point in Y to perform out-of-sample
    %               prediction. For example, the first out-of-sample 
    %               prediction uses all data in Y before pstart. 
    % n_factor:     number of factor lags             
% OUTPUT
    % MSE:          the value of Mean Squared Error from forecasting  
    
function predict_mse = mse(Y,X,h,pstart,n_factor)

T = size(Y,1);

% Initiate a matrix to store forecast results
Y_predict = NaN*ones(T-h-pstart+2,1);

% Recursively forecast y 
for t = pstart-1:T-h
    Y_reg = Y(h+n_factor:t);
    X_reg = X(1:t-h-n_factor+1,:);
    b = ((X_reg'*X_reg)^(-1))*(X_reg'*Y_reg);
    
    % Form prediction
    X_forecast = X(t-n_factor+1,:);
    Yhat = X_forecast * b;
    Y_predict(t-pstart+2) = Yhat;
end

% Calculate mean squared error
predict_mse = mean((Y(pstart+h-1:end,:)-Y_predict).^2);
end
