clear
close all
clc

% =========================================================================
% DESCRIPTION 
% This script loads in a FRED-MD dataset, processes the dataset, and then
% estimates factors. This code replicates output of codes published on  
% https://research.stlouisfed.org/econ/mccracken/fred-databases. 

% INPUT FILES 
%       FRED-MD data in CSV format (published on 
%       https://research.stlouisfed.org/econ/mccracken/fred-databases).

% OUTPUT FILES
%       outilers'year'.csv: CSV file with data after transformation and
%       outlier removal. 
%       
%       n_'year'.csv: a CSV file recording the number of outliers removed 
%       for each series 
%
%       mnfactor'year'.csv: CSV file containing the first factor estimates 
%       and dates.
%
%       mnfactor_'year'.mat: MAT file that contains the first factor 
%       estimates, dates, and numerical dates. 
%
%       outdates'year'.csv: CSV file contains indicator of location of 
%       outliers.
% -------------------------------------------------------------------------
% BREAKDOWN OF THE SCRIPT
% 
% Part 1: Load and label FRED-MD data.
%
% Part 2: Process data -- transform each series to be stationary and remove
%         outliers.
%
% Part 3: Estimate factors and compute R-squared and marginal R-squared. 
%
% Part 4: Save outputs
% -------------------------------------------------------------------------
% AUXILIARY FUNCTIONS
% List of auxiliary functions to be saved in same folder as this script.
%
%   prepare_missing() - transforms series based on given transformation
%       numbers
%
%   remove_outliers() - removes outliers
%
%   factors_em() - estimates factors
%
%   mrsq() - computes R-squared and marginal R-squared from factor 
%       estimates and factor loadings
%
% -------------------------------------------------------------------------
% NOTES
% Authors: James Hamilton and Jin Xi
% Date: 12/27/2023
% Version: MATLAB 2022a
% Required Toolboxes: None
%
% -------------------------------------------------------------------------
% PARAMETERS TO BE CHANGED

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

% Information criterion used to select number of factors; for more details,
% see auxiliary function factors_em()
%   1 --> information criterion PC_p1
%   2 --> information criterion PC_p2
%   3 --> information criterion PC_p3
jj=2;

% Maximum number of factors to be estimated; if set to 99, the number of
% factors selected is forced to equal 8
kmax=8;

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

version = 2024;
% use version = 2015 for estimating factors with 2015-04.csv vintage
% use version = 2024 for 2024-12.csv
%      version = 2023 for 2023-11.csv
%      version = 1997 to calculate factors only using data through end of 1996

if version == 2015
   enddate = '01-Dec-2014';
   csv_in = '2015-04.csv';
   app_in = 'FRED-MD_historic_appendix.csv';
elseif version == 2023
   enddate = '01-Jul-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,',');

% variable names
series=dum.textdata(1,2:end);

% transformation numbers
tcode=dum.data(1,:);

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

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

% extract 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);

% extract sample period
seriesdata = seriesdata(Loc_start:Loc_end,:);
dates = dates(Loc_start:Loc_end,:);

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

% 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_mn= [first_year+first_month/12:1/12:final_year+final_month/12]';
% This follows McCracken-Ng convention that 1963 means Dec 1962
dates_num_mn = dates_num_mn - 1/12; % This converts to 1963 means Jan 1963

% =========================================================================
% PART 2: PROCESS DATA

% transform raw data to be stationary using auxiliary function
% prepare_missing()
yt=prepare_missing(seriesdata,tcode);

% reduce sample to usable dates: remove first two months because some
% series have been first differenced
T=size(dates,1);
yt=yt(3:T,:);
dates=dates(3:T,:);
dates_num_mn=dates_num_mn(3:T,:);

% remove outliers using auxiliary function remove_outliers(); see function
% or readme.txt for definition of outliers
%   data = matrix of transformed series with outliers removed
%   n = number of outliers removed from each series
[data,n]=remove_outliers(yt);

output_cyc = strcat('mn_cyc_comp_',string(version),'.csv');
writematrix(data,output_cyc,'Delimiter',',');

% =========================================================================
% PART 3: ESTIMATE FACTORS AND COMPUTE R-SQUARED 

% estimate factors using function factors_em()
%   ehat    = difference between data and values of data predicted by the 
%             factors
%   Fhat    = set of factors
%   lamhat  = factor loadings
%   ve2     = eigenvalues of data'*data
%   x2      = data with missing values replaced from the EM algorithm
[ehat,fhat_mn,lamhat,ve2,x2] = factors_em(data,kmax,jj,DEMEAN);

% compute R-squared and marginal R-squared from estimated factors and
% factor loadings using function mrsq()
%   R2      = R-squared for each series for each factor
%   mR2     = marginal R-squared for each series for each factor
%   mR2_F   = marginal R-squared for each factor
%   R2_T    = total variation explained by all factors
%   t10_s   = top 10 series that load most heavily on each factor
%   t10_mR2 = marginal R-squared corresponding to top 10 series
%             that load most heavily on each factor 
[R2,mR2,mR2_F,R2_T,t10_s,t10_mR2] = mrsq(fhat_mn,lamhat,ve2,series);

f1_mn= fhat_mn(:,1); % first factor estimates

% dates in the format year.month
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

% =========================================================================
% PART 4: SAVE OUTPUTS
output_outliers = strcat('outliers',string(version),'.csv');
output_n = strcat('n_outliers',string(version),'.csv');
output_outdates = strcat('outdates',string(version),'.csv');
output_fhat = strcat('mnfactor',string(version),'.csv');

writematrix(data,output_outliers,'Delimiter',',')
writematrix(n,output_n,'Delimiter',',')
datex = datex(3:size(data,1)+2,1);
datez = (datex*ones(1,size(data,2))).*isnan(data);
writematrix(datez,output_outdates,'Delimiter',',')
writematrix([datex f1_mn],output_fhat,'Delimiter',',')

% =========================================================================
% SAVE DATA TO .MAT FILE: 

output_mat = strcat('mnfactor','_',string(version));
save(output_mat, 'f1_mn', 'dates_num_mn', 'dates') 

output_mat_all = strcat('mnfactor_all','_',string(version));
save(output_mat_all, 'fhat_mn', 'dates_num_mn', 'dates') 
