clear
close all
clc

% =========================================================================
% DESCRIPTION 
% This script loads in a FRED-MD dataset, detrends using OLS, collects a 
% balanced subpanel, and estimates factors
% -------------------------------------------------------------------------
% INPUT FILES 
%       FRED-MD data in CSV format (published on 
%       https://research.stlouisfed.org/econ/mccracken/fred-databases).

% OUTPUT FILES
%       hx_cyc_comp_'year'.xlsx: spreadsheet with cyclical component 
%       estimates using Hamilton & Xi's OLS detrending method.
%
%       hxfactor'year'.csv: CSV file that contains the first k factor 
%       estimates from the cyclical component estimates. 
%
%       dropped_'year'.xlsx: spreadsheet that indicates whether each series 
%       is dropped due to missing data or is kept and used for estimation.
%
%       hxfactor'year'.mat: MAT file that contains the first factor 
%       estimates from the cyclical component estimates.
%
%       hx_cyc_comp_'year'.mat: MAT file with cyclical component 
%       estimates using Hamilton & Xi's OLS detrending method.
% -------------------------------------------------------------------------
% Breakdown of the script
% 
% Part 1: Load and label data
%
% Part 2: Process data -- transform each series to be stationary with 
% OLS regression 
%
% Part 3: Estimate factors of the OLS-detrended data.
%
% Part 4: Plot the factors. 
%
% -------------------------------------------------------------------------
% Auxiliary functions
% List of auxiliary functions to be saved in same folder as this script.
%
%   create_yX() - create dependent variable y and regressor matrix X given
%   prediction horizon and number of lags
%
%   fhat() - estimates factors of stationary data 
%
%   lsq() - performs least square regression on non-missing values
% -------------------------------------------------------------------------
% Control parameters

p = 12; % number of lags included for ols detrending procedure
h = 24; % forecasting time horizon (in months)
k = 8; % Number of factors to be estimate

% Type of transformation performed on each series before factors are
% estimated
%   0 --> no transformation
%   1 --> demean only
%   2 --> demean and standardize
type = 2;

real_time = 0;
    % set real_time = 0 for baseline analysis
    % set real_time = 1 to estimate parameters only up to date end_real
end_real = '01-Dec-2014';  
     % end date for OLS estimation to be used when real_time = 1


% Specify starting date of sample period for factor estimation
startdate = '01-Jan-1960';

version = 2024;
% use version = 2015 for estimating factors with 2015-04.csv vintage
%      version = 1997 for 2015-04.csv vintage but use data only through 1996:12
%      version = 2023 for 2023-11.csv
%      version = 2024 for 2024-12.csv


% Specify ending date of sample period, name of data file, name of FRED-MD
% appendix file 
if version == 2015
   enddate = '01-Dec-2014';
   csv_in = '2015-04.csv';
   app_in = 'FRED-MD_historic_appendix.csv';
elseif version == 2023
   enddate = '01-Jun-2023';
   csv_in = '2023-11.csv';
   app_in = 'FRED-MD_updated_appendix.csv'; 
elseif version == 1997  
   enddate = '01-Dec-1996';
   csv_in = '2015-04.csv';
   app_in = 'FRED-MD_historic_appendix.csv';
elseif version == 2024
   enddate = '01-Sep-2024';
   csv_in = '2024-12.csv';
   app_in = 'FRED-MD_historic_appendix.csv';

else
   disp('This case not programmed')
   stop
end

% =========================================================================
% PART 1: LOAD AND LABEL DATA

% Load data from CSV file
dum = importdata(csv_in,','); 
app = importdata(app_in,',');

% Transformation code (provided by McCraken & Ng)
tcode=dum.data(1,:);

% Variable names
varnames = dum.textdata(1,2:end);

% Label variables with descriptions in appendix 
% Remark: later vintages add one more variable
if version <= 2015+3/12
    des = app.textdata(2:end-1,4); 
else
    des = app.textdata(2:end,4); 
end

% Numerical values of variables 
seriesdata = dum.data(2:end,:);

% Dates
dates = dum.textdata(3:end,1);

% Read and format date vector 
dates = datetime(dates,'InputFormat','MM/dd/uuuu');

% Find the location of the starting datetime
[~,Loc_start] = ismember(startdate,dates);
[~,Loc_end] = ismember(enddate,dates);
if real_time == 1
    [~,Loc_real] = ismember(end_real,dates);
end

% Specified sample period series data
seriesdata = seriesdata(Loc_start:Loc_end,:);
dates = dates(Loc_start:Loc_end,:);

% Month/year of final observation
final_datevec=datevec(dates(end,1));
final_month=final_datevec(2);
final_year=final_datevec(1);

% Month/year of first observation
first_datevec= datevec(dates(1,1));
first_month = first_datevec(2);
first_year = first_datevec(1);


% Numerical dates (monthly) are of the form YEAR+(MONTH-1)/12
% e.g. Jan 1963 is represented as 1963, March 1970 is represented as
% 1970+2/12
dates_num_hx = [first_year+first_month/12:1/12:final_year+final_month/12]';
% This follows McCracken-Ng convention that 1963 means Dec 1962
dates_num_hx = dates_num_hx - 1/12;
% This converts to 1963 means Jan 1963


% =========================================================================
% Part 2: Process data -- transform each series to be stationary with OLS regression
% 1) transformation code and 2) OLS regression respectively. 

vals = [];    % balanced detrended data set will be collected in vals
kept_names = []; % names of retained variables
N0 = size(seriesdata,2);
kept_series = zeros(1,N0);   % kept_series is indexes of series that are included in balanced panel
dropped_series = zeros(1,N0); % dropped_series is indexes of series that are dropped from balanced panel
for i = 1:N0

    % z = seriesdata(first_obs:last_obs,i);
    z = seriesdata(:,i);
    % ------------------------------------------
    % check if series contains missing values for proposed balanced sample
    if sum(isnan(z)) > 0
       dropped_series(i) = 1;
    else
       kept_series(i) = 1;
    end
 
    if kept_series(i) == 1
    % ------------------------------------------
        % identify whether log, level, or other transformation based on McCracken-Ng classifications
        if tcode(i) == 7
            z(2:end) = z(2:end)./z(1:end-1);
            z(1) = 1;
        elseif tcode(i) < 4
           % no transformation in this case  
        else
           z = 100*log(z);
        end

    % ------------------------------------------
    % detrend with OLS
    if real_time == 0
       [y,X] = create_yX(z,h,p); 
       [bhat,resids]=lsq(y,X); 
    elseif real_time == 1
       [y,X] = create_yX(z(1:Loc_real),h,p);
       [bhat,resids2]=lsq(y,X);
       [y,X] = create_yX(z,h,p);
       resids = y-X*bhat;
    end 
    vals= [vals resids];
    kept_names = [kept_names varnames(i)];
    end

end

dropped_summary = num2cell([dropped_series' kept_series']);
des = des(kept_series == 1);
% =========================================================================
% Part 3: Factor estimation

[lambda_hx,fhat_hx,eigenvalue_hx]=fhat(vals,k,type);

% normalize signs
lambda_hx(:,1) = -lambda_hx(:,1);
fhat_hx(:,1) = -fhat_hx(:,1);


% =========================================================================
% Part 4: Plot figures 

% remove the first h+p periods due to OLS detrending
dates_num_hx = dates_num_hx(h+p:end);
dates = dates(h+p:end);


% plot the first two factors based on OLS detrended series
subplot(2,1,1)
plot(dates,fhat_hx(:,1),'LineWidth',1,'Color','b')
hold on
xlim([min(dates) max(dates)])
title('PC1 of OLS-detrended series')
hold off
subplot(2,1,2)
plot(dates,fhat_hx(:,2),'LineWidth',1,'Color','b')
hold on
xlim([min(dates) max(dates)])
title('PC2 of OLS-detrended series')
hold off

% Part 5: Export results to Excel and MAT files

num_big = (final_year - first_year +1)*12;
datex = zeros(num_big,1);
   
 i = 1;
   for t = first_year:final_year
     for month = 1:12
        datex(i,1)=t+month/100;
        i=i+1;
     end
   end
datex = datex(h+p:end,1);
datex = datex(1:size(fhat_hx,1));

% -------------------------------------------------------------------------
% SAVE output 

if real_time == 0
   output_factor = strcat('hxfactor',string(version),'.csv');
   output_cyc = strcat('hx_cyc_comp_',string(version),'.xlsx');
   output_dropped = strcat('dropped_',string(version),'.xlsx');

   % save factor estimates to CSV file
   writematrix([datex fhat_hx],output_factor,'Delimiter',',')

   % save dropped variable summary to spreadsheet
   dropped_title = {'dropped', 'kept'};
   writecell(varnames',output_dropped,'Sheet',1,'Range','A2');
   writecell(dropped_summary,output_dropped,'Sheet',1,'Range','B2')
   writecell(dropped_title,output_dropped,'Sheet',1,'Range','B1')

   % save cyclical component estimates to spreadsheet
   cell_cyc = num2cell([datex vals]);
   writecell(kept_names,output_cyc,'Sheet',1,'Range','B1')
   writecell(cell_cyc,output_cyc,'Sheet',1,'Range','A2')

   % -------------------------------------------------------------------------
   % SAVE DATA TO .MAT FILE: 
   % convert to dating convention used in forecast program
   %dates_num_hx = dates_num_hx + 1/12;

   filename_f = strcat('hxfactor','_',string(version));
   save(filename_f,'fhat_hx','lambda_hx', 'dates_num_hx', 'dates') 
   filename_cyc = strcat('hx_cyc_comp_',string(version));
   save(filename_cyc, 'varnames','des','vals','dates_num_hx', 'dates')

elseif real_time == 1
   fhat_real_hx = fhat_hx;
   filename_f = strcat('hxfactor_real','_',string(version));
   save(filename_f,'fhat_real_hx','lambda_hx', 'dates_num_hx', 'dates') 
end
