Plain and Target-PCA factors for mixed-frequency Data
using the MATLAB Nowcasting Suite
This example uses various pre-processing fucntions, available in the Nowcasting Suite, with the goal of constructing a set of monthly target-relevant PCA factors, setting the quarterly US GDP growth rate as the target variable. The high-frequency factors can then be used in a MIDAS-type predictive regression to backcast, nowcast and/or forecast the low-frequency target, in real-time.
The Nowcasting Suite has been tailored for macroeconomic forecasting handling real-time mixed-frequency datasets, using MATLAB's native time-stamped data format, called timetable. The pcaF() function can implement target-PCA with either hard- or soft-thresholding ala Bai-Ng (2008). The example presented below uses hard-thresholding. Furthermore, pcaF() can create target-factors for both single and mixed-frequency data. That is, it can construct monthly target-relevant factors for a quarterly target by down-sampling the transformed high-frequency regressors into the target's lower frequency, conduct the pre-selection step at the target's sampling frequency, and then, collect the selected predictors at their highest available frequency.
Contents
1. Downloading Data in Real-time
2. Transforming the target variable
3. Time-Dummies for Target Variable (to control for irregular periods)
4. Pre-processing the Monthly Dataset
4.1.a. Log Transformations
4.1.b. Stationarity Transformations
4.2. Outlier Detection and Adjustment
4.3. Fill Missing Values using EM ala SW(2002)
5. Estimate Plain- and Target-PCA factors
5.1. Plain-PCA monthly factors
5.2. Target-PCA monthly factors with Hard-thresholding
5.3. Target-PCA quarterly factors with Hard-thresholding
6. Extracting Leads & Lags (MIDAS & single-frequency)
6.1 Leads for Quarterly Predictors (X's)
6.2 MIDAS Leads & Lags for Monthly Factors (F's)
6.3 Legendre MIDAS for Monthly Factors (F's)
List of Functions
1. fred_api()
2. transLHS(), StepsAhead4Nowcast(), trans()
3. LHStimeDummies()
4.
4.1.b. trans2I0(), trans2I0_ADF(), DFGLS()
4.2. outl_adj()
4.3. MissFill()
5. pcaF()
6.
6.1 DownSample(), trans2MIDAS()
6.2 AFTERtrans2MIDAS(), flatify()
6.3 trans2MIDASLegendre()
A short description of selected functions included in the Nowcasting Suite, together with a brief overview of their features, can be found here.
%Add the Nowcasting Suite folder to MATLAB's search-path
addpath(genpath(fullfile(pwd, 'MaFunctions')));
warning('off','all') %Suppress all MATLAB warnings
1. Downloading Data in Real-time
Let's start by downloading a set of mixed-frequency (30 monthly & the quarterly target) series in real-time, using the FRED API, accessed via the fred_api() function. The series mnemonics to be downloaded are stored in codes.M and codes.Q.
%Dataset: 20 Macro and Financial Series and 1 Quarterly
%To change the length of the series to train/estimate the models,
strtDate = datetime(2010,1,1);
%%% 1) MONTHLY DATA Codes %%%
codes.M = ["AAAFFM"; "AWHNONAG"; "BAAFFM"; "BOGMBASE"; "IPFINAL"; "CES9093000001"; ...
"CPIAUCSL"; "CPIULFSL"; "CUSR0000SAD"; "DMANEMP"; "DTCOLNVHFNM"; "IPDCONGD"; ...
"GS5"; "INDPRO"; "IPDMAT"; "LNS12032194"; "LNS13023557"; "LNS13023569"; ...
"LNS14000026"; "NDMANEMP"; "IPFPNSS"; "PPIIDC"; "T1YFFM"; "T5YFFM"; ...
"TB3MS"; "TB6SMFFM"; "TOTRESNS"; "UEMP5TO14"; "USPBS"; "USTPU"];
%%% 2) QUARTERLY DATA Codes %%%
codes.Q = ["GDPC1"; "A014RE1Q156NBEA"; "B021RE1Q156NBEA"; "CNCF"; "COMPRMS"; ...
"DHCERG3Q086SBEA"; "DHUTRG3Q086SBEA"; "DONGRG3Q086SBEA"; "COMPRNFB"; ...
"SLCE"; "OPHMFG"; "OPHNFB"; "OPHPBS"; "OUTBS"; "PCESV"; "RCPHBS"; "TLBSNNCB"; ...
"TNWBSNNB"; "ULCNFB"; "Y033RC1Q027SBEA"];
%%% Categories of Downloaded Variables %%%
%- Monthly (codes.M)
cat.M = ["Interest Rates"; "Labor"; "Interest Rates"; "Money and Credit"; "Industrial Production"; ...
"Labor"; "Prices"; "Prices"; "Prices"; "Labor"; "Money and Credit"; "Industrial Production"; ...
"Interest Rates"; "Industrial Production"; "Industrial Production"; "Labor"; "Labor"; ...
"Labor"; "Labor"; "Labor"; "Industrial Production"; "Prices"; "Interest Rates"; "Interest Rates"; ...
"Interest Rates"; "Interest Rates"; "Money and Credit"; "Labor"; "Labor"; "Labor"];
%- Quarterly (codes.Q)
cat.Q = ["Output, income, consumption"; "Output, income, consumption"; "Output, income, consumption"; ...
"Non-Household Balance Sheets"; "Earnings & Productivity"; "Prices"; "Prices"; "Prices"; ...
"Earnings & Productivity"; "Output, income, consumption"; "Earnings & Productivity"; ...
"Earnings & Productivity"; "Earnings & Productivity"; "Output, income, consumption"; ...
"Output, income, consumption"; "Earnings & Productivity"; "Non-Household Balance Sheets"; ...
"Non-Household Balance Sheets"; "Earnings & Productivity"; "Output, income, consumption"];
The table below shows the distribution of the downloaded variables into the various economic categories, which closely follow those of McCracken and Ng in the FRED-MD dataset.
%ALL Categories (i.e. the calssification system used consists of 11 categories)
tabulate(categorical([cat.M; cat.Q]))
%-------------------%
%%% Download Data %%%
%-------------------%
%%% 1) MONTHLY DATA %%%
[X.M, Meta, release] = fred_api(codes.M, 'StartDate',strtDate);
%The reference date in FRED API is the 1st of the Month. Change it to LAST day of month.
X.M.Time = dateshift(X.M.Time,'end','month',0);
RelDates = table(Meta.SeriesID, release.lastupdate', release.lastdate', Meta.Frequency, 'VariableNames',["SeriesID","LastUpdate","LastOb","Freq"]);
The following two tables contain selected metadata for the monthly panel. Specifically, the first reports the FRED mnemonic and the description of each series, the assigned economic category, and its measurement units. The second table relates to the publication lags and the resulted ragged-edge structure of the monthly dataset. It reports the last available observation for each series, as well as the release date of that last observation.
[~,loc] = ismember(codes.M, Meta.SeriesID); %loc: index in RIGHT for each value in LEFT
[Meta(loc,"SeriesID") table(cat.M,'VariableNames',"Category") Meta(loc,["Title"])]
[Meta(:,["SeriesID","Units"]) RelDates(:,["LastUpdate","LastOb"])]
%%% 2) QUARTERLY DATA %%%
[X.Q, Meta_q, release_q] = fred_api(codes.Q, 'StartDate',strtDate);
%The reference date in FRED API is the 1st day of the Quarter. Change it to the LAST day.
X.Q.Time = dateshift(X.Q.Time,'end','quarter',0);
Meta = [Meta; Meta_q];
RelDates_q = table(Meta_q.SeriesID, release_q.lastupdate',release_q.lastdate', Meta_q.Frequency, 'VariableNames',["SeriesID","LastUpdate","LastOb","Freq"]); %Keep to see which series have been updated since last time
RelDates = [RelDates; RelDates_q];
clear Meta_q release_q RelDates_q release loc
Xraw = X; %Rename the dataset
2. Transforming the target variable
Function transLHS() applies the desired transformation on the target variable. The specs passed (in struct 'LHS') below will transform the GDP from levels to annualized quarter-on-quarter (QoQ) discretely compounded growth rates.
Specifically, setting tcode=11 calculates the percentage changes, while tcode=1 would take the 1st differences, tcode=2 the 2nd differences, and tcode=21 would transform by taking the first differences of the percentage changes.
Setting XoX=LP determines the base period, with LP corresponding to 'last period' (QoQ, MoM), as opposed to SPLY or 'same period last year' (YoY). Apart from the non h-specific transformations (QoQ/MoM/YoY), transLHS() can also accommodate h-specific transformations, such as the h-quarters/months changes or percentage changes.
LHS.name = 'GDPC1'; %GDP growth (QoQ% Annualized; simple-rates)
LHS.f = 'Q';
LHS.tcode = 11; LHS.hSpecificY = false; LHS.Percent = true;
LHS.Annualize = true; LHS.log = false; LHS.XoX = 'LP';
LHS.log4lags = true; %Transformation for the AR-lags (output: LHS.y4lags)
LHS.DaysDelay = 65;
LHS.flow = true; %Stock/flow type referes to the variable in LEVELS (not the transformed version)
%Description (for plotting):
LHS.P.Title = 'US Real GDP growth';
LHS.P.Des = 'Annualized QoQ%';
The GDP series will normally contain provisional 1st and 2nd estimates on the first and second months of the following quarter, respectively, before the final figure is released, on the last month of the quarter following the reference quarter. There are 2 modelling choices for those preliminary figures:
1) Keep them among the realized values, and simply forecast the next period.
2) Delete them from the sample, and forecast them.
transLHS() allows to automate the detection and removal of provisional values with the optional input 'RmvProvisional'. Below, we are taking the 2nd route by making use of 'RmvProvisional'. By setting DaysDelay=65, observations released within the first 65 days after the closing of the reference quarter are labeled as provisional and deleted from the resulted transformed target (LHS.y). This ensures that only GDP values released in the 3rd month of the quarter remain in the sample. To follow the 1st route, we can simply set RmvProvisional=[] or DaysDelay=[].
% Treating Provisional/Preliminary/Advance figures.
RelDates(RelDates.SeriesID==LHS.name,:) %The release dates of the LAST of ob of the series
lastObRelDate = RelDates{RelDates.SeriesID==LHS.name,"LastUpdate"};
rmvpro.DateReleased = lastObRelDate; %We have set LHS.DaysDelay=65
%NB: If I use input 'RmvProvisional' to remove Provisional obs from the LHS,
%we also need to remove these obs from the inputted matrix X for consistency
%(- that's the purpose of the 2nd output of transLHS()).
hmax = 3;
[LHS, Xraw.(LHS.f)] = transLHS(LHS,Xraw.(LHS.f),'hall',1:hmax,'RmvProvisional',rmvpro);
clear lastObRelDate rmvpro hmax
%Display the last 5 obs
LHS.pr.y(end-4:end,1);
LHS.y(end-4:end,1); %<- The last (provisional) observation has been removed
LHS.y4lags(end-4:end,1); %MoM change for AR-lags (-don't multiply by 100)
%Plot the transformed LHS series
p = stackedplot(LHS.y,{LHS.y.Properties.VariableNames(1)});
p.Title = [LHS.P.Title ' (' LHS.P.Des ')'];
p.DisplayLabels = {''};
clear p
Apart form the transformed LHS, function transLHS() also returns a separate transformed version of the target (in LHS.y4lags), to be used as autoregressive (AR) terms. For target variables that are transformed in 1st or 2nd differences (not levels), regardless of the selected LHS transformation, the transformation for the AR lags is always:
(1) with respect to the previous period (i.e. MoM/QoQ, even if the user selected YoY for the LHS),
(2) non-annualized, and
(3) expressed as percentages in decimal form (not in '%' form, i.e. without dividing by 100).
tail(LHS.y4lags,5) %Transformation for AR terms
To handle the mismatch of frequencies, conveniently transLHS() also stores the date corresponding to the last available observation of the LHS, 
in LHS.T.
in LHS.T.LHS.T %Last LHS observation
Next, we use function StepsAhead4Nowcast() to define the number of h-steps ahead required to predict the current period (-nowcasting). The resulted 'h' measures in terms of the target-variable's sampling frequency. Therefore, here 'h' corresponds to quarters-ahead.
%today's date
asofDate = datetime('today', 'Format','dd/MM/yyyy');
%Find the 'h' needed in order to calculate the NOWCAST.
%NB: If instead of the NOWCAST, you want the 2-quarters ahead FORECAST,
% then simply do: h_now=StepsAhead4Nowcast(), and then set h=h_now+2;
h = StepsAhead4Nowcast(LHS.T,'Now',asofDate,'Frequency',LHS.f);
%fcastDate: The DATE to be forecasted (i.e. the DATE corresponding to h-steps ahead)
fcstDate = dateshift(LHS.T,'end','quarter', h);
msg = ['As of ' char(asofDate) ', the last available observation for the ' LHS.P.Title ' is for ' char(LHS.T) '.' newline 'The nowcast (' char(fcstDate) ') is obtained by setting h=' num2str(h) '.'];
disp(msg)
Function transLHS() makes an internal call to another Nowcasting Suite function, trans() which applies the desired transformation. Here is an example of using trans(), to convert headline and core CPI indicies to MoM% annualized inflation in discretely compounding rates. Similarily to transLHS(), setting XoX='SPLY' would transform CPI to YoY% inflation, while setting Log=true, would calculate MoM% annualized inflation in continuously compounding rates.
x0 = fred_api(["CPIAUCSL", "CPILFESL"], 'StartDate',datetime(2020,1,1));
x1 = trans(x0,11,'M','h',[],'hSpecificY',false,'Log',false,'Percent',true,'Annualize',true,'XoX','LP');
p = stackedplot(x1(:,:),{["CPIAUCSL", "CPILFESL"]});
p.DisplayLabels = {''};
p.Title = 'US CPI Inflation';
p.AxesProperties(1).LegendLabels = ["Headline","Core"];
clear p x0 x1 msg
3. Time-Dummies for Target Variable (to control for irregular periods)
In order to control for irregular periods (structural breaks, outliers) in the target variable, we can create a set of time-dummies capturing such periods, and subsequently, control for these by adding the dummies as additional control variables in the correpsonding predictive regressions. The same control dummies can accompany the thresholding regressions of the targeted-factors in the predictor pre-selection step. Function LHStimeDummies() uses outlier detection based on IQR to create a different set of time-dummies for each h-step ahead version of the LHS. Automating the time-dummy creation process, allows to easily employ different target-variables.
The table that follows shows that setting the threshold for the detection of irregular periods to 1.8, the function labels a total of 7 observations as 'irregular'. Those 7 irregular observations are clustered by LHStimeDummies() into 4 periods of 'irregular' activity, each made of consecutive 'irregular' quarters, 2 of which span multiple quarters. Specifically, over the period 2010-2023 covered in this sample, the 4 periods where GDP growth took extreme values are: 2014Q1, 2020Q1-2020Q3, 2021Q2 and 2021Q4-2022Q1. The set of dummies created by function LHStimeDummies() will consist of 4 controlling dummies (- some covering consecutive quarters), one for each of these periods.
%Detect irregular periods for Y
[Dh, PeriodsDetected, outliersY] = LHStimeDummies(LHS,'DetectOptions',struct('Threshold',1.8));
[outliersY(:,h) LHS.y(:,h)] %See the outliers (for the nowcast period)
Dh{h}(15:50,:)
clear PeriodsDetected outliersY
4. Pre-processing the Monthly Dataset
4.1. Log and Stationarity [I(0)] Transformations
4.2. Outlier Detection and Adjustment (i.e. set to NaN)
4.3. Fill Missing Values using EM ala SW(2002)
%%%%%%%%%%%%%%%%%%
% Pre-processing %
%%%%%%%%%%%%%%%%%%
f='M'; %Changing this to 'Q', would run the same pre-preprocessing for panel Q.
Xt = Xraw.(f);
[~,locmeta] = ismember(Xt.Properties.VariableNames, Meta.SeriesID); %index in RIGHT for each value in LEFT
MetaF = Meta(locmeta,:); %Frequency-specific 'Meta'
clear locmeta Xraw
4.1. Log and Stationarity Transformations
4.1.a. Log Transformations
The first pre-processing step is to choose in which of the predictors to apply a log-transformation before differencing, to ensure stationarity. Following Stock and Watson (2002), we take the natural logarithm in all non-negative series (-rule 1), that are not already in rates or percentage units (-rule 2).
To see which variables are already in rates or percentage units, we will need to read through the table containing the variable descriptions. Alternatively, to automate the 2nd rule's selection, we can search and exclude all the series that contain the word 'percent' in the Units variable of the metadata table of the monthly dataset.
%1a) Log Transformations
%Rule 1: Series that are non-negative (+ive or 0)
small = 1e-6; % Value close to zero
canlog = (min(Xt{:,:}) > small); %series NOT close to 0
%Rule 2: Series that are not already in rates or percentage units
MetaF(:, ["SeriesID","Units"])
notpct = ~contains(MetaF.Units, "percent",'IgnoreCase',true)'; %tilda = negation operation
%Combine Rules 1 & 2
what2log = notpct & canlog;
%Check the decisions and which X's will be logged.
array2table([string(Xt.Properties.VariableNames); notpct; canlog; what2log]','VariableNames',["Series","Non-Rate","Non-negative","What2Log"])
tcode.(f).log = what2log;
clear small canlog pos idx notpct MetaF
4.1.b. Stationarity Transformations
The vector containing the decisions for the series over which to take the natural logarithm, is passed as input in function trans2I0() which implements an iterative unit root testing procedure, in order to determine the apporpriate number of differences required to achieve stationarity, for each of the series in the monthly panel. Prior to transforming the series, trans2I0() applies the log-transformations.
The unit root testing procedure of trans2I0() is based on the DF-GLS test of Elliott, Rothenberg & Stock (1996). The function allows to choose the optimal lag-length k* in the underlying ADF regressions, with 3 different methods:
1. SIC,
2. Modified Akaike information criterion (MAIC), ala Ng & Perron (2001; Econometrica),
3. Sequential-t algorithm ala Ng & Perron (1995; JASA),
The stationarity transformations for the monthly panel below are based on the SIC, while the rejection of the null is set to the 5% significance level.
%1b) Stationarity [I(0)] Transformations
%Transform to Stationary based on iterative Unit Root (UR) tests
%(DFGLS with optimal truncation lag order selected via SIC, considering 0-Kmax lags)
opDFGLS = struct('Trend',1,'Kmin',0,'Kmax',[],'MethADF',1,'MethCVal',{'CL'},'MethLagOrder',{'SIC'},'SignifLevel',5);
[Xst, objDFGLS] = trans2I0(Xt,'How','URTest','What2Log',what2log,'TestOptions',opDFGLS,'MaxDiffs',3,'Diagnostics',true);
tcode.(f).D = objDFGLS.TransCodes;
clear what2log Xt
Function trans2I0() contains a ceiling on the number of times a series is allowed to be differenced until it achieves stationarity, with the maximum being set at 3. The choice of the maximum differences is in line with the stationarity transformation codes of McCracken and Ng (2016). The first column of the following table shows how many times each of the series was differenced (0=levels; 1=1st diffs; 2=2nd diffs etc.), and the second column shows whether the series managed to reject the null hypothesis of a unit root.
dfgls = array2table([objDFGLS.TransCodes' objDFGLS.URrejected'],'VariableNames',["TransCodes","URrejected"],'RowNames',Xst.Properties.VariableNames)
Moreover, trans2I0() saves the output of the DFGLS test for each series in levels, as well as, after each consecutive differencing, allowing an in-depth look of what is happening in the series, with each subsequent application of the test. The test's output closely follows the output of the 'dfgls' command in Stata. The next 4 tables report the DFGLS test results, one after each differencing applied on the 2nd series (AWHNONAG), until it achieved stationarity. The corresponding DFGLS test statistic is reported on the 3rd column of the test output. Recall that the desicion rule for the rejection of the null for a unit root, was set at the 5% level of significance with the optimal truncation lag order selected via the SIC (-1st row of the test output table). We can see that in the 4th table, corresponding to the series differenced 3 times, the unit root null is rejected even at the 1%.
series = 2;
objDFGLS.TestResults{series}{:}
clear series
%CHECK: Stationarity using SIMPLE ADF test (NOT DFGLS)
[Xtest, objADF] = trans2I0_ADF(Xst,'What2Log',[],'Diagnostics',false); %NB: Do NOT log here
hasUR = (objADF.TransCodes ~= 0);
find(hasUR);
dfgls(hasUR,:);
%objADF.TransCodes>0 => series did NOT have its UR rejected/removed by the DFGLS
chck = array2table([objADF.TransCodes', objDFGLS.URrejected'],'VariableNames',["ADF_TransCodes","DFGLS_URrejected"],'RowNames',Xst.Properties.VariableNames);
chck(hasUR,:);
cases2check = ~(objADF.TransCodes==0 & objDFGLS.URrejected==1);
chck(cases2check,:);
dfgls(cases2check,:);
Xst(:,hasUR) = Xtest(:,hasUR); %Replace the UR-series with the ADF-stationary series
tcode.(f).D(hasUR) = objDFGLS.TransCodes(hasUR) + objADF.TransCodes(hasUR); %Final tcodes of series in Xst
clear Xtest hasUR opDFGLS cases2check chck objADF objDFGLS dfgls
transLHS() makes an internal call to another function of the Nowcasting Suite, the DFGLS() function. DFGLS() conducts the unit root test in individual series and returns the test output. Let's run an example using the GDP series in continuously compounding (that is, log-) growth rates. The settings selected below replicate Stata's default behaviour for its 'dfgls' command. Specifically:
- Setting Trend=1, the series is GLS-transformed (de-trended) using a linear trend prior running the ADF regressions. This sets the alternative hypothesis (Ha) to the series being stationary around a linear trend. Alternatively, setting Trend=0, the GLS-tranformation only de-means the variable and the Ha is set to Yt being stationary with no linear time trend.
- The minimum number of lags on the ADF regression is set to 1, while the maximum (Kmax=[]) is set to the value proposed by Schwert (1989): 
- The critical values for the DF-GLS (aka ERS) test statistic are based on Cheung & Lai (1995). The test output also reports the Elliott, Rothenberg & Stock (1996) critical values.
%(A)DF OLS regressions: dYt = rho*Yt-1 + c1*dYt-1 + c2*dYt-2 +...+ ckmax*dYt-kmax + et
%When k=0 the regression equation becomes the DF regression (i.e. no dYt lags).
%where k is the # of dYt-j terms in the (A)DF regression.
x0 = X.Q(:,"GDPC1"); %GDP (in Levels)
x1 = trans(x0,11,'Q','h',[],'hSpecificY',false,'Log',true,'Percent',false,'Annualize',false,'XoX','LP');
dfglsout = DFGLS(x1,'Trend',1,'Kmin',1,'Kmax',[],'MethADF',0,'MethCVal','CL')
clear dfglsout x0 x1
The empty cells at the bottom of the output table that refer to the Ng and Perron (1995) sequential-t algorithm, imply that the 'last' dYt-k term at the different iterations for the lag order, was never statstically significant at the pre-selected 10% threshold.
4.2. Outlier Detection and Adjustment
Function outl_adj() detects outliers in the individual series based on each series' Interquartile Range (IQR), and adjusts them by either replacing with the median, or setting them to missing (NaN) to be filled at later pre-processing steps. A data point (x) is labeled as an outlier if its absolute deviation from the median is larger than a user-selected multiple (-threshold) of the series' IQR:
The threshold is normally set to 10.
Optional input 'IncludeUp2' allows the user to select the desired date up to which to perform outlier detection and correction. This feature of outl_adj() can be useful if one wants to leave intact the observations at the bottom end of the sample, that will be used to calculate the forecast. Such practise could lead to superior forecasting accuracy, when extreme economic conditions are unravelling, and extreme values appear at the bottom end of the sample.
For example, setting IncludeUp2=dateshift(LHS.T,'end','month',-12) will exclude from the outlier adjusting process any 'leading' observations, as well as any lagging observations pretaining to the 12 months preceding the last available observation of the target. Assuming that the forecast is estimated using V=12 lags of the monthly predictors (one year's worth of past information, plus any available leading observation), this will exclude any observations that will be used to calculate the forecast.
Alternatively, setting IncludeUp2=0 searches and fixes outliers to the entire series, which is the standard practise. That's what is done below, in this example.
%------------------------------%
% Outlier Detection/Adjustment %
%------------------------------%
%Step 1: Detect outliers & replace with missing
inclup2 = 0;
thr = 10; %IQR multiple
adjm = 0; %replace detected outliers with NaNs
[Xst, outlier] = outl_adj(Xst, thr, adjm,'IncludeUp2',inclup2);
%row-wise sum to check # of problematic obs per period
array2timetable(sum(outlier{:,:},2), 'RowTimes',outlier.Time,'VariableNames', "Outliers")
%column-wise sum to check # of problematic obs per indicators
array2table(sum(outlier{:,:},1)','VariableNames',"Outliers",'RowNames',Xst.Properties.VariableNames)
clear outlier adjm thr inclup2
%Check ROWS with too many NaNs
array2timetable(sum(isnan(Xst{:,:}),2),'RowTimes',Xst.Time,'VariableNames', "NaNs");
%NB: We will have a spike in the number of outliers/NaNs around the COVID-lockdowns period & in the BOTTOM because of the ragged edge
X0.(f) = Xst; %X0.M = Xst;
clear Xst X
4.3. Fill Missing Values using EM ala SW(2002)
In this section, we use function MissFill() to implement the EM algorithm of SW(2002) in order to fill any missing values that are due to:
1) the restricted history of the series,
2) asynchronous publication releases and the implied ragged-edges, and
3) any outliers detected and set to missing by outl_adj(), in the steps above.
Before filling the missing values of the transfromed monthly panel (X0.M), lets take a look on its structure.
R = 25;
tail(X0.M(:, 1:R),10)
%%% Fill the Missing Values (outliers & top NaNs) using EM ala SW(2002) %%%
%Step 2: Correct Outliers & the NaNs
%I am filling up to the LHS i.e. excluding the ragged-edge (fillup2=LHS.T)
opK = struct('Method','BN2002','Kmax',8,'Criterion',2);
X.M = MissFill(X0.M,'Method','EM','Options',opK,'FillUp2',LHS.T);
%Fill the monthly panel that will be used to construct the monthly PCA factors.
%Setting FillUp2=0 will fill the entire matrix. This will include any
%available obs ahead of the forecast period (LHS.T+h).
X4pca.M = MissFill(X0.M,'Method','EM','Options',opK,'FillUp2',0);
clear opK
%NB: The 'X4pca' dataset will be 'balanced' (meaning with no missing values), but X.M will have ragged-edges (unbalanced).
Both of the tables below show the bottom end of the monthly predictors' matrix. The first table shows a version of the predictors matrix where we have filled the missing values (ragged-edges and any detected outliers) all the way up to the last available observation of the LHS, partly leaving the ragged-edge intact.
The first dataset (X.M) can be used to assess the (leading) informational content of the series that become timely available, while the second dataset (X4pca.M) is used below to apply targeted-PCA and extract a set of monthly factors. The first table will therefore show the ragged-edge of the selected predictors, while the second has its ragged-edge filled.
tail(X.M(:, 1:R),10)
tail(X4pca.M(:, 1:R),10)
It might be sensible for one to need to exclude from the sample any timely observations that are available for the period ahead of the period to be forecasted (> LHS.T + h). In other words, if the task is to backcast 2024Q1, you might want to exclude monthly observations for 2024M4, 2024M5 etc., from the sample. That would make the predictors' matrix horizon-specific. For the predictors' matrix (X.M) we can achieve that by simply deleting the observations ahead of the forecast period (LHS.T+h), i.e. by executing X.M(X.M.Time > fcstDate,:)=[]. Similarily, if we don't want the estimated PCA factors to include observations ahead of the forecast period, we can do so when filling the matrix over which we will apply PCA (X4pca.M) by setting MissFill(.,FillUp2=fcstDate).
MissFill() can also fill missing values using iterated (as opposed to direct) multistep univariate AR(P) models for each of the series, where the optimal lag order P* can be selected with a number of methods, including BIC, AIC, as well as K-fold CV ala Bergmeir et al (2018). Because of the time-series regression nature of the underlying model, this method can only fill missing values in the middle (due to outliers) and at the bottom (due to publication lags), but not at the top of the series.
5. Estimate Plain- and Target-PCA factors
In this section we use the filled monthly panel (X4pca.M) to create two sets of monthly factors. First, we estimate a standard set of factors, and then, we use hard-thresholding, ala Bai-Ng (2008), to select target-relevant predictors over which PCA is subsequently applied. The mismatch in the frequencies between the target and the predictors is handled by temporally-aggregating the stationary monthly panel to the sampling frequency of the LHS. The target-relevant predictors selected during the pre-selection step, are kept at their highest available frequency, over which PCA is applied resulting into the 2nd set of monthly target-relevant factors.
In both PCA applications we use the 
information critirion of Bai-Ng (2002) to determine the optimal number of factors, setting the maximum allowed number of factors to 
.
information critirion of Bai-Ng (2002) to determine the optimal number of factors, setting the maximum allowed number of factors to .5.1. Plain-PCA monthly factors
clear X0
%Set the options for selecting the optimal number of factors (K*)
opK = struct('Method','BN2002','Kmax',8,'Criterion',2);
%- MONTHLY plain-Factors
opTrgt = []; %<- Do NOT implement Target-factors
[F.p, objF.p] = pcaF(X4pca.M,'treatX','M','KstarOptions',opK,'KstarAllMethods',false,'Groups',cat.M,'Diagnostics',true,'addGraph',true,'ExcludefromX',[],'Target',opTrgt);
Function pcaF() creates a number of post-estimation graphs and tables, that can help better understand the underlying structure of the estimated factors, both in terms of their explanatory power and in terms of factor interpretability.
The first graph is a Pareto chart showing the incremental (and cummulative) variance explained by each (additional) factor. Since factors are organized in decreasing order of their respective eigenvalues, the underlying ordering of the factors reflects their importance in explaining the set of predictors X. As such, the 1st factor explains more of the variability of X compared to the 2nd factor, and the 2nd more than the 3rd, and so on.
The second graph provides interpretability for each factor based on the economic categories that each series has been assigned, and the explanatory power of that particular factor on the series, as measured by its mR2. Specifically, mR2(i,k) measures the incremental explanatory power of the k-th factor for the i-th series.
The two tables that follow (PostPCA.T2 & PostPCA.T3) are the table versions of graphs 1 & 2, respectively. The first shows the variance explained by each factor, as well as the cummulative explained variance. The second table reflects each factor's interpretation. Specifically, for each factor Fk, it calculates the contributions of the various categories by summing the mR2's of the regressors assigned to the same category, and then dividing by the sum of all mR2(i,k)'s for that k.
objF.p.PostPCA.T2
objF.p.PostPCA.T3
The next post-estimation table (PostPCA.T1) sorts the series wrt their mR2, for each of the factors. The table below shows the 10 most important series for each of the factors.
head(objF.p.PostPCA.T1,10)
Finally, pcaF() stores the NxK* table with the factor loadings (Lambda), and two TxN timetables containing the standardized (Uz) and un-standardized (Ux) idiosyncratic components of the Factor model, 
.
.objF.p.Lambda %Table with the estimated loadings
tail(objF.p.Ux,10) %Timetable of the idiosyncratic components
5.2. Target-PCA monthly factors with Hard-thresholding
Target-relevant factors are obtained by using a pre-selection step which determines the predictors that are most relevant to the variable that we want to predict. Since standard PCA is an unsupervised learning technique, this additional pre-selection step is done in an effort to imporve the forecasting accuracy of models based on standard PCA factors, such as the factor-augmented autoregressive model (F-AR), which is also commonly referred to as the Diffusion Index (DI) model, popularized by Stock-Watson (2002).
Hard-thresholding, was introduced by Bai-Ng (2008). The pre-selection step involves estimating an OLS regression for each predictor, and picking the predictors that are above a pre-defined statistical significance threshold. pcaF() allows the user to:
1) introduce optional control variables, such as time-dummies and AR-lags in the pre-selection regressions, and to
2) estimate horizon-specific pre-selection regressions, so that the selected predictors are better suited for the period one wants to forecast.
The OLS regression of the hard-thresholding pre-selection step takes the form:
on AR-lags + time-dummies + Contemporaneous 
The Nowcasting Suite works as a ecosystem of related functions, where the output of one function can directly become the input of another. Conveniently, outputs LHS.y4lags of transLHS() and Dh of LHStimeDummies(), containing the AR-lags and time-dummies corresponding to the LHS of interest, can directly be passed to pcaF(), to be added as control variables in the pre-selection, without further processing, or the need to re-align/time-sync the data matricies, manually.
In hard-thresholding, the threshold corresponds to the critical value that defines the rejection area for the hypothesis of the insignificance of the i-th predictor. For example, it is common to set the threshold to 1.96, which is the critical value of the Standard Normal (and Student-t with df→∞) distribution corresponding to upper-tail (1-side) probability of 2.5%. This in turn will pick the X's with a t-statistic higher that 1.96 or lower than -1.96, implying that the PCA will be applied to all those predictors that are statistically significant at the 5% level of significance, or in other words, those that reject the 
against a 2-sided 
, with 95% confidence. Finally, pcaF() will include the LHS by default, among the pre-selected predictors over which it applies PCA. To exclude the LHS, we can simply add it in 'ExcludefromX'.
against a 2-sided , with 95% confidence. Finally, pcaF() will include the LHS by default, among the pre-selected predictors over which it applies PCA. To exclude the LHS, we can simply add it in 'ExcludefromX'.%- MONTHLY target-Factors
%Thresholding settings for target-PCA
prescrn = struct('How','Individual','Statistic','t');
topx = struct('Method','Threshold','k',1.96);
P=2; %Add 2 AR-lags as controlling variables
rollwin = []; %Set estimation window for the pre-selection OLS regression
opTrgt = struct('LHS',LHS,'h',h,'YLags',P,'Dummies',{Dh},'RollingWindow',rollwin,'What2Screen',prescrn,'Screen',topx);
[F.t, objF.t] = pcaF(X4pca.M,'treatX','M','KstarOptions',opK,'KstarAllMethods',false,'Groups',cat.M,'Diagnostics',true,'addGraph',true,'ExcludefromX',[],'Target',opTrgt);
After implementing target-PCA with hard-thresholding, pcaF() stores a table containing the t-statistics and the corresponding p-values for all the predictors in the monthly panel.
objF.t.Target.stat
Another table stored by pcaF() is the table with only the pre-selected predictors, ranked descendingly based on their absolute t-statistics from the most to the least target-relevant.
objF.t.Target.seriesrnk
Before extracting the leads and lags for the selected factors, let us see in a plot how the two sets of factors compare. Below we plot the 1st principal component of the monthly plain- and target-PCA factors, together with the GDP growth. The figure also contains shaded areas representing the recession periods, as determined by the NBER Business Cycle Dating committee. To add in a signle plot the quarterly GDP growth together with the monthly factors, we fill the months within the quarter by repeating the corresponding value.
%Fill the quarterly series by repeating within-period values
y2m = retime(LHS.y(:,h), F.t.Time, 'next'); %Fill within-quarter months with next value
%Print the last 10 rows
tail(y2m(:, :),10)
%Standardize the GDP to have comparable scale with the F's
y2m{:,:} = zscore4nan(y2m{:,:});
F.t.Properties.VariableNames = F.t.Properties.VariableNames + "T";
F.p.Properties.VariableNames = F.p.Properties.VariableNames + "P";
TTm = [y2m F.t F.p]; %Combine GDP & Factors in a single timetable
%----------------------%
% NBER recession Dates %
%----------------------%
%To add recession shaded areas in the graph
%0) Download the (monthly) NBER recession dates from FRED (mnemonic: USREC)
%1) Create an event table
%2) Pass the event table as a property to the timetable (to be ploted).
% When you display the timetable, you can also see the event labels.
%3) Passing the timetable (with the attached event table) to stackedplot()
% will plot events as shaded regions. (since R2023b)
recessions = fred_api("USREC", 'StartDate',datetime(1980,1,1));
recessions{:,:} = dum_me2(recessions{:,:}); % dum_me2() numerates the recessions
rvec = recessions{:,:};
recessionStrt = recessions.Time(arrayfun(@(z) find(rvec==z,1,'first'),1:max(rvec)));
recessionStrt = dateshift(recessionStrt,'end','month',0); %Convert to END-of-month
recessionEnd = recessions.Time(arrayfun(@(z) find(rvec==z,1,'last'),1:max(rvec)));
recessionEnd = dateshift(recessionEnd,'end','month',0);
recDates = array2table([recessionStrt recessionEnd],'VariableNames',["Start","End"])
EventDates = eventtable(recDates{:,1},'EventEnds',recDates{:,2}, 'EventLabels',"Recession"); %event table
TTm.Properties.Events = EventDates; %Pass the events table as a property to the timetable
clear recessions rvec recessionStrt recessionEnd recDates EventDates
%Plot the 1st target- and plain-Factor, against the GDP growth series
p = stackedplot(TTm,{["F1P" "F1T" LHS.name+"_h"+h]});
p.AxesProperties(1).LegendLabels = ["Plain-F1","Target-F1","GDP growth"];
p.AxesProperties(1).LegendLocation = 'SouthWest';
p.DisplayLabels = {''};
p.Title = 'US GDP growth & Monthly PCA Factors';
p.LineProperties(1).LineStyle = {'-','-', '-.'};
p.LineProperties(1).Marker = {'none','none','none'};
p.LineProperties(1).LineWidth = [0.5 0.5 1.2];
p.LineProperties(1).Color = [0.96, 0.49, 0.04; 0.06 0.55 0.91; 0 0 0]; %1=orange; 2=blue; 3=black
The following two tables report the correlation coefficients for the 3 series that were ploted above. Since those are based on repeating GDP values, correlation coefficients with the GDP, are going to be low.
vnm = ["GDPC1_h2","F1P","F1T"];
array2table(corr(TTm{:,vnm},'rows','pairwise'), 'VariableNames',vnm)
exclcovid = (TTm.Time<=datetime(2020,1,1) | TTm.Time>=datetime(2022,7,1));
array2table(corr(TTm{exclcovid,vnm},'rows','pairwise'), 'VariableNames',vnm)
clear p TTm y2m exclcovid vnm prescrn topx
5.3. Target-PCA quarterly factors with Hard-thresholding
Finally, pcaF() stores a copy of the temporally-aggegated (that is, the monthly-averaged-to-quarterly) predictors matrix used in the pre-selection step. We can use these to create a set of quarterly target-factors. The table below shows the quarterly versions of the pre-selected predictors.
Xq_targt = objF.t.Target.X4ps(:,objF.t.Target.series); %The pre-selected predictors in quarterly frequency
tail(Xq_targt,10)
%Apply PCA to the quarterly version of the predictors selected by the previous pcaf() call
opTrgt = [];
[~,loc] = ismember(Xq_targt.Properties.VariableNames, codes.M); %loc: index in RIGHT for each value in LEFT
[Fq, objFq] = pcaF(Xq_targt,'treatX','Q','KstarOptions',opK,'KstarAllMethods',false,'Groups',cat.M(loc),'Diagnostics',true,'addGraph',true,'ExcludefromX',[],'Target',opTrgt);
%Combine GDP & Quarterly Factors in a single timetable
yz = zscore4nan(LHS.y(:,h)); %Standardize the GDP
TTq = synchronize(yz, Fq, 'union','fillwithmissing');
vnm = ["GDPC1_h2","F1"];
array2table(corr(TTq{:,vnm},'rows','pairwise'), 'VariableNames',vnm)
p = stackedplot(TTq,{["F1" LHS.name+"_h"+h]});
p.AxesProperties(1).LegendLabels = ["Target-F1","GDP growth"];
p.AxesProperties(1).LegendLocation = 'SouthWest';
p.DisplayLabels = {''};
p.Title = 'US GDP growth & Quarterly PCA Factor';
p.LineProperties(1).LineStyle = {'-', '-.'};
p.LineProperties(1).LineWidth = [0.5 1.2];
p.LineProperties(1).Color = [0.96, 0.49 0.32; 0 0 0];
clear loc opTrgt objFq TTq p yz vnm Xq_targt opK Fq
6. Extracting Leads & Lags (MIDAS & single-frequency)
6.1 Leads for Quarterly Predictors (X's)
6.2 MIDAS Leads & Lags for Monthly Factors (F's)
6.3 Legendre MIDAS for Monthly Factors (F's)
In the nowcasting literature, 'leads' refer to the extracted time-series lags pertaining to 'leading periods', after applying the time-series lag operator to series that become more timely available compared to the target-variable's publication release schedule. These lags comprised of leading periods, contain ahead-in-time information that can potentially contribute to better predicting the relatively back-in-time target variable.
Using the trans2MIDAS() funtion of the Nowcasting Suite, the forecaster can distinguish the leads and lags of each predictor, allowing to experiment with different structures of the MIDAS polynomials.
Optional inputs allow the researcher to choose between two methods for extracting the leads & lags for each variable: (1) The standard method of extracting in total V leads + lags (i.e. applying the MIDAS time-series operator to obtain the first V lagged series of the high-frequency covariate); (2) Considering ALL leads + V lags, ala Andreou, Ghysels and Kourtellos (2013, JBES). The latter (AGK2013) method could be argued to be preferable, compared to the standard method of creating the MIDAS conditioning set, for the following reasons:
- If for an Xi we have many leads, using the first method might result in including little or NO lagged information at all, for that Xi.
- The latter method ensures that the lagging information is the same for all the X's, meaning that it covers the exact same period, past the last available LHS observation. For instance, assuming a quarterly-monthly mix, 1 quarter of past monthly information (prior to the last observed LHS) will be included if we set V=3, an entire year if V=12, two years if V=24 etc. To illustrate this point, let’s assume that the target variable's last observed value corresponds to Q4 (i.e., 31st Dec 2022). Also, say that for X1 we have 0 leads (i.e., last observed value for X1 is in 31/12/2022), while for X2 we have 5 leading months, meaning that we have 5 monthly observations ahead of the last observed LHS value (i.e., last observed value for X2 is in 31/5/2023). Extracting the last V=12 leads/lags (using the standard method) will have a year's worth of information for both X1 and X2, but only for X1 that information will correspond to the entire last year of the LHS variable.
For example, assuming that we have predictors at monthly and daily frequency (and a quarterly target), one might consider adding to the same model a restricted lead and lag (e.g. Almon) polynomial for the daily indicator to ensure parsimony, and an unrestricted polynomial for the leads of the monthly indicator, while excluding any monthly lags altogether. Furthermore, being able to distinguish between what comprises leading and what lagging information allows us to see whether potential forecasting gains are due to timely information flowing into the model or if it's just the outcome of flexibility in the temporal-aggregation functions that summarize the lagged information.
trans2MIDASLegendre() works in the same manner as trans2MIDAS(), but employs the newly introduced MIDAS-Legendre polynomials further reducing the dimensionality of covariates, while also adding all the benefits of having orthogonal covariates.
Finally, the two functions also synchronize/align the different conditioning variables (X,D,LY) to each other, and to the target variable (y), so that they correspond to the appropriate periods for the estimation of the (direct-forecasting) model producing the desired h-steps ahead forecast.
%adding code soon
References
Andreou E., Ghysels E., and Kourtellos A., (2013) "Should Macroeconomic Forecasters Use Daily Financial Data and How?", Journal of Business & Economic Statistics, 31(2): 240-251.
Bai, J., and Ng, S. (2002). "Determining the number of factors in approximate factor models." Econometrica, 70(1), 191-221.
Bergmeir, C., Hyndman, R. J., and Koo, B. (2018). "A note on the validity of cross-validation for evaluating autoregressive time series prediction." Computational Statistics & Data Analysis, 120: pp. 70-83.
Cheung and Lai (1995) Lag order and critical values of a modified Dickey–Fuller test. Oxford Bulletin of Economics and Statistics 57: 411–419.
Elliott, Rothenberg and Stock (1996) Efficient Tests for an Autoregressive Unit Root, Econometrica, Vol. 64, No. 4, pp. 813-836.
McCracken, M. and Ng, S., (2016), "FRED-MD: A Monthly Database for Macroeconomic Research," Journal of Business & Economic Statistics, 34(4): pp. 574-589.
Ng and Perron (1995) Unit root tests in ARMA models with data-dependent methods for the selection of the truncation lag. Journal of the American Statistical Association 90: 268–281.
Ng and Perron (2001) Lag length selection and the construction of unit root tests with good size and power. Econometrica 69: 1519–1554.
Schwert (1989) Tests for unit roots: A Monte Carlo investigation. Journal of Business and Economic Statistics 2: 147–159.
Stock, J. H., & Watson, M. W. (2002). "Macroeconomic forecasting using diffusion indexes." Journal of Business & Economic Statistics, 20(2), 147-162.