The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
110,304 characters
Principal Component Analysis for a Mix of Stationary and Nonstationary Variables
\title{Principal Component Analysis for a Mix of Stationary and
Nonstationary Variables\thanks{
An earlier version of this paper was circulated under the title, "Principal
Component Analysis for Nonstationary Time Series." We thank Serena Ng, Mark
Watson, Ken West, and anonymous referees and editors for helpful comments.
Xi gratefully acknowledges financial
support from the National Natural Science Foundation of China (No. 72503234).
Authors are listed in alphabetical order; all contributed equally to this
work.}}
\author[]{James D. Hamilton\thanks{
Corresponding author. Email: \texttt{[email removed].}} \and Xinwei Ma\thanks{
Email: \texttt{[email removed].}} \and Jin Xi\thanks{
Email: \texttt{[email removed].
}}}
\date{May 28, 2022 \\
Revised: August 17, 2026}
\maketitle
\begin{abstract}
\doublespacing This paper develops a procedure for uncovering the common
cyclical factors that drive a mix of stationary and nonstationary variables.
The method does not require knowing which variables are nonstationary or the
nature of the nonstationarity. An application to the FRED-MD macroeconomic
dataset demonstrates that the approach offers similar benefits to those of
traditional principal component analysis with some added advantages.
\end{abstract}
\affil[1]{University of California at San Diego}
\affil[2]{University of
California at San Diego} \affil[3]{Center for Forecasting Science, Academy
of Mathematics and Systems Science, Chinese Academy of Sciences}
\bigskip \thispagestyle{empty} \vspace{2em}
\textbf{Keywords}: principal components, nonstationary, economic activity
indexes \bigskip
\clearpage
\pagenumbering{arabic}
\section{Introduction}
Principal component analysis (PCA) has become a key tool for building
dynamic models of vector time series with a large cross-sectional dimension.
A typical application first subtracts the sample mean from each variable and
divides the demeaned variables by their sample standard deviations. PCA
finds linear combinations of the standardized variables that have maximum
sample variance. These principal components are then used to build dynamic
models for each individual series. For surveys of PCA and its usefulness in
economics see \cite{BaiNg2008} and \cite{StockWatson2016Handbook}.
One difficulty with PCA is that many of the time series encountered in
economics and finance are nonstationary. For a nonstationary variable, the
population mean is undefined and the sample standard deviation diverges to
infinity as the number of time-series observations gets large. \cite
{onatski2021spurious} detailed some of the problems that can arise from
trying to apply PCA to nonstationary data. The typical solution to this
problem is for researchers to examine each series individually by hand to
determine the transformation that needs to be applied to that series before
calculating principal components of the set of variables.
This approach has two shortcomings. First, while for some variables it may
be fairly clear what transformation is necessary to achieve stationarity,
for others it is far from obvious. For example, interest rates have a strong
downward trend since 1980. Should yields be treated as stationary? If
nonstationary, should we take their first differences or deviations from a
time trend before performing PCA? Some authors work with first differences
of inflation and the unemployment rate while others leave those variables as
is. Many decisions like these have to be made before applying PCA to large
data sets. In this paper, we propose an automatic procedure that allows the
researcher to avoid these judgment calls.
A second fundamental problem is the appropriateness of the methodology
itself. Suppose we somehow overcame the first problem and knew for certain
the true nature of the trend in each individual series. Suppose for
illustration we knew correctly that the first variable $y_{1t}$ is a
stationary $AR(1)$ process with autoregressive coefficient ${\Greekmath 011A} =0.99$
while the second variable $y_{2t}$ is a random walk. The currently
prescribed procedure would instruct the researcher to use the first variable
as is and the second variable in first differences. But while the levels of $
y_{1t}$ and $y_{2t}$ exhibit very similar properties, the level of $y_{1t}$
and first difference of $y_{2t}$ are radically different. Should we expect
that there is some linear combination of $y_{1t}$ and $\Delta y_{2t}$ that
can summarize the common economic drivers behind the two variables? If
differencing is the appropriate transformation for a random walk, it seems
we should be using some similar transformation for an $AR(1)$ process whose
root is close but not quite equal to unity.
This paper proposes an approach to PCA that solves these problems. We first
note as in \cite{hamilton2018you} that it is possible to use OLS to estimate
an $h$-period-ahead forecast of the level of any variable as a linear
function of its own lags without knowing the nature of the nonstationarity.
Moreover, the errors from these linear forecasts are stationary for a broad
class of underlying nonstationary processes. Our proposal is to use PCA to
identify common factors behind the forecast errors. Specifically, we
estimate an OLS regression of $y_{it}$ on $(1,y_{i,t-h},y_{i,t-h-1},\dots
,y_{i,t-h-p+1})^{\prime }$, where $i=1,\dots ,N$ represents the index of the
variable, $h$ is the forecasting horizon, and $p$ is the number of lags used
for the forecast. We then calculate principal components of the residuals.
This approach solves both of the problems identified above. The procedure is
fully automatic and treats every variable in the same way. The
transformation of $y_{it}$ is a continuous function of the estimated
autoregressive coefficients and implies no discontinuity in the way that
persistent stationary series are treated as the largest autoregressive root
tends toward unity.
\cite{hamilton2018you} suggested that a two-year-ahead forecast error could
be interpreted as the stationary cyclical component of the series. Our
empirical application follows this suggestion and uses\thinspace $h=24$ for
monthly data to uncover the common cyclical factors behind a range of
economic and financial indicators. We demonstrate that using $h=24$ also
offers a practical solution to the problem of outliers which has posed a
severe challenge to using PCA on macroeconomic observations in 2020.
\cite{BaiNg2004} and \cite{barigozzi2021large} developed methods that are
suitable when the observed data are an unknown mix of $I(0)$ and $I(1)$
variables. Their approaches start by applying PCA to first-differences of
all of the variables whether stationary or not. \cite{BaiNg2004} estimated
factors for the original data by accumulating factors estimated from the
differenced data. \cite{barigozzi2021large} estimated factors by applying
the loadings estimated from differenced data to the detrended levels of the
original data. \cite{barigozzi2021large} showed that their procedure can
estimate the set of factors that account for common stochastic trends in the
original data. In contrast to these approaches, the goal of our analysis is
to uncover the cyclical components that are driving both the stationary and
nonstationary variables.
Our results are also useful for estimation of dynamic factor models. One
common approach to estimating dynamic factor models is to apply a
transformation that is believed to make the individual series stationary,
calculate principal components of the transformed data, and then fit a
vector autoregression to the estimated factors (\cite
{stock2002forecasting,StockWatson2016Handbook}; \cite{BaiNg2002,
BaiNgEcta2006}). In contrast, our approach allows the researcher to
implement the initial transformations without the need for subjective
judgments about each individual series.
Section \ref{sec_detrend} discusses the use of forecasting regressions to
isolate a stationary component of a possibly nonstationary time series.
Section \ref{sec_c_obs} describes how PCA could be used to uncover the
common factors behind the forecast errors if we somehow could know the true
process followed by each variable. Section \ref{sec_c_not_observed} analyzes
the case when PCA is applied to the residuals from estimated OLS regressions
and establishes consistency of the method. Simulations in Section \ref
{sec_simulations} report results for a variety of possible data-generating
processes. We find that using two-year-ahead forecast errors works well when
we have 50 years of data and that $h=1$ performs quite well even in much
shorter samples. Section \ref{sec_fred_md} illustrates the promise of our
approach in an empirical analysis using the FRED-MD large macroeconomics
data set. Proofs of theorems and additional simulations and empirical
results are provided in an online supplement.
\section{Isolating a stationary component from a nonstationary series}
\label{sec_detrend} Let $y_{it}$ denote the level of a possibly
nonstationary variable at time $t$. Collect its $p$ most recent values as of
date $t-h$ along with a constant term in a vector {${z}
_{i,t-h}=(1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})^{\prime }.$} The
population linear projection of $y_{it}$ on {${z}_{i,t-h}$} is given by {$
\mathbb{P}(y_{it}|{z}_{i,t-h})={\Greekmath 010B} _{i0}^{\prime }{z}_{i,t-h}$} where the
coefficient ${\Greekmath 010B} _{i0}$ is defined as the vector that minimizes the
expected squared error of a forecast of $y_{it}$ based on a linear function
of {${z}_{i,t-h}$}:
\begin{equation}
{\Greekmath 010B} _{i0}=\arg \underset{{\Greekmath 010B} }{\min }E(y_{it}-{\Greekmath 010B} ^{\prime }{z}
_{i,t-h})^{2}. \label{eqn_pop_lin_proj_definition}
\end{equation}
For example, if $y_{it}$ is covariance stationary with {$E({z}_{i,t-h}{z}
_{i,t-h}^{\prime })$} nonsingular, the population linear projection
coefficient is given by
\begin{equation}
{\Greekmath 010B} _{i0}=\left[ E({z}_{i,t-h}{z}_{i,t-h}^{\prime })\right] ^{-1}E({z}
_{i,t-h}y_{it}). \label{eqn_pop_lin_proj_stationary}
\end{equation}
If $y_{it}$ is ergodic for second moments, the population parameter ${\Greekmath 010B}
_{i0}$ can be consistently estimated by an OLS regression of $y_{it}$ on {${z
}_{i,t-h}$}:
\begin{equation*}
\hat{{\Greekmath 010B}}_{i}=\left[ \mathop{\textstyle \sum }\nolimits_{t=1}^{T}{z}_{i,t-h}{z}
_{i,t-h}^{\prime }\right] ^{-1}\left[ \mathop{\textstyle \sum }\nolimits_{t=1}^{T}{z}
_{i,t-h}y_{it}\right] =\left[ T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}{z}_{i,t-h}{z}
_{i,t-h}^{\prime }\right] ^{-1}\left[ T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}{z}
_{i,t-h}y_{it}\right] \overset{p}{\rightarrow }{\Greekmath 010B} _{i0}.
\end{equation*}
Note that there is no assumption in the definition (\ref
{eqn_pop_lin_proj_definition}) nor in its stationary solution (\ref
{eqn_pop_lin_proj_stationary}) that the process is linear or that it can be
characterized by an ARMA representation of known order. There may be better
forecasts of $y_{it}$ that could be obtained using a nonlinear function or
more lags like $y_{i,t-h-p},y_{i,t-h-p-1},...$ But there is an optimal
forecast within the class of linear forecasts that use only $p$ lags. For a
covariance-stationary process, the optimal forecast within that class is
characterized by (\ref{eqn_pop_lin_proj_stationary}).
A solution to (\ref{eqn_pop_lin_proj_definition}) also exists if $y_{it}$ is
nonstationary but its first difference $\Delta y_{it}$ is covariance
stationary. This can be seen from the accounting identity that holds for any
time series:
\begin{equation}
y_{it}=y_{i,t-h}+\Delta y_{i,t-h+1}+\Delta y_{i,t-h+2}+\cdots +\Delta y_{it}.
\label{eqn_I_1_identity}
\end{equation}
For an\thinspace $I(1)$ process,\footnote{
Note that the space spanned by $z_{i,t-h}$ is the same as the space spanned
by $(z_{i,t-h}^{\dag },y_{i,t-h})$ which is bigger than the space spanned by
$z_{i,t-h}^{\dag }.$ However, the coefficient on $y_{i,t-h}$ in the
projection $\mathbb{P}(\Delta y_{it}|{z}_{i,t-h}^{\dag },y_{i,t-h})$ is
zero, because multiplying the $I(1)$ variable $y_{i,t-h}$ by any nonzero
value would produce a forecast of $y_{it}$ that has infinite MSE. For this
reason $\mathbb{P}(\Delta y_{it}|z_{i,t-h})=\mathbb{P}(\Delta y_{it}|{z}
_{i,t-h}^{\dag }).$}
\begin{equation}
\mathbb{P}(y_{it}|{z}_{i,t-h})=y_{i,t-h}+\mathbb{P}(\Delta y_{i,t-h+1}|{z}
_{i,t-h}^{\dag })+\mathbb{P}(\Delta y_{i,t-h+2}|{z}_{i,t-h}^{\dag })+\cdots +
\mathbb{P}(\Delta y_{it}|{z}_{i,t-h}^{\dag }) \label{eqn_lin_proj_I_1}
\end{equation}
for ${z}_{i,t-h}=(1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})^{\prime },$ ${z}
_{i,t-h}^{\dag }=(1,\Delta y_{i,t-h},\Delta y_{i,t-h-1},...,\Delta
y_{i,t-h-p+2})^{\prime }$ and \linebreak $\mathbb{P}(\Delta
y_{i,t-s}|z_{i,t-h}^{\dag })=\left[ E(\Delta y_{i,t-s}{z}_{i,t-h}^{\dag
\prime })\right] \left[ E({z}_{i,t-h}^{\dag }{z}_{i,t-h}^{\dag \prime })
\right] ^{-1}{z}_{i,t-h}^{\dag }.$ Moreover, the population linear
projection coefficient ${\Greekmath 010B} _{i0}$ can again be consistently estimated by
an OLS regression of the level of $y_{it}$ on $(1,y_{i,t-h},y_{i,t-h-1},
\allowbreak ...,y_{i,t-h-p+1})^{\prime }.$ The intuition for this is that an
OLS levels regression minimizes the sample analog to the population
minimization in (\ref{eqn_pop_lin_proj_definition}):
\begin{equation}
\hat{{\Greekmath 010B}}_{i}=\arg \underset{{\Greekmath 010B} }{\min }T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}(y_{it}-{\Greekmath 010B} ^{\prime }{z}_{i,t-h})^{2}.
\label{eqn_OLS_objective_function}
\end{equation}
If the OLS estimate $\hat{{\Greekmath 010B}}_{i}$ picks up the unit root, the value of (
\ref{eqn_OLS_objective_function}) converges to a finite number as $T$ gets
large. If it does not, the value of (\ref{eqn_OLS_objective_function}) would
diverge to infinity. Thus in large samples, OLS estimation is dominated by
the incentive inherent in the objective function of OLS to remove any
nonstationary features of the data. More formally, the fitted values from
the regression (\ref{eqn_OLS_objective_function}) are numerically identical
to the fitted values from an OLS regression of $y_{it}$ on {$(y_{i,t-h},{z}
_{i,t-h}^{\dag ^{\prime }}).$} In the latter regression, the OLS coefficient
on the lagged level $y_{i,t-h}$ will converge in probability to one whenever
the dependent variable is an $I(1)$ process. In fact, the OLS estimate of
this parameter is superconsistent, converging at rate $T$ rather than at
rate $\sqrt{T}$, and the asymptotic distribution of the other coefficients
is identical to what would result if we regressed $\Delta y_{it}$ on {${z}
_{i,t-h}^{\dag };$} see \cite{hamilton2018you}.
\cite{hamilton2018you} showed that these results generalize to any
nonstationary process that is covariance stationary around a deterministic
polynomial in time $t$ provided that the order of the polynomial is no
greater than the number of lags $p,$ or if the process is $I(d)$ with the $d$
th difference of $y_{it}$ covariance stationary for some $d$ no greater than
$p.$ One does not need to know the order of the polynomial or the value of $
d $ in order to consistently estimate the population linear projection
coefficient ${\Greekmath 010B} _{i0}$ using a regression of the level of $y_{it}$ on a
constant and its $p$ most recent levels as of date $t-h.$
\cite{hamilton2018you} also showed that the residual from the population
linear projection, {$c_{it}=y_{it}-{\Greekmath 010B} _{i0}^{\prime }{z}_{i,t-h},$} is
covariance stationary for any of the above nonstationary processes. For
example, for an $I(1)$ process, we see from (\ref{eqn_I_1_identity}) and (
\ref{eqn_lin_proj_I_1}) that
\begin{eqnarray*}
c_{it} &=&y_{it}-\mathbb{P}(y_{it}|{z}_{i,t-h}) \\
&=&\left[ \Delta y_{i,t-h+1}-\mathbb{P}(\Delta y_{i,t-h+1}|{z}_{i,t-h}^{\dag
})\right] +\left[ \Delta y_{i,t-h+2}-\mathbb{P}(\Delta y_{i,t-h+2}|{z}
_{i,t-h}^{\dag })\right] + \\
&&\cdots +\left[ \Delta y_{it}-\mathbb{P}(\Delta y_{it}|{z}_{i,t-h}^{\dag })
\right]
\end{eqnarray*}
is covariance stationary for any finite $h.$ It is thus possible to isolate
a stationary component of a large range of nonstationary processes from the
residuals of an OLS regression of $y_{it}$ on $
(1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})^{\prime }$ without knowing
whether the variable $y_{it}$ is stationary or nonstationary.\footnote{
Specifically, if either: (i) $y_{it}$ is stationary around a deterministic
polynomial function of time of order $d_{i}\leq p$ satisfying
\begin{equation*}
T^{-1/2}\sum_{s=1}^{[Tr]}(y_{it}-{\Greekmath 010E} _{i0}-{\Greekmath 010E} _{i1}t-{\Greekmath 010E}
_{i2}t^{2}-\cdots -{\Greekmath 010E} _{i}t^{d_{i}})\overset{d}{\rightarrow } {\Greekmath 0121}
_{i}W_{i}(r)
\end{equation*}
where $[Tr]$ denotes the largest integer no greater than $Tr$, $W_{i}(r)$
denotes standard Brownian motion, and $\overset{d}{\rightarrow } $ denotes
convergence in distribution; or alternatively if (ii) $d_{i}$ differences of
$y_{it}$ are stationary for some $d_{i}\leq p$ satisfying
\begin{equation*}
T^{-1/2}\sum_{s=1}^{[Tr]}(\Delta ^{d_{i}}y_{it}-{\Greekmath 0116} _{i})\overset{d}{
\rightarrow } {\Greekmath 0121} _{i}W_{i}(r);
\end{equation*}
then \cite{hamilton2018you} showed that $c_{it}$ is stationary and that the
estimated coefficient $\hat{{\Greekmath 010B}}_{i}$ from an OLS regression gives a
consistent estimate of the population parameter ${\Greekmath 010B} _{i0}.$}
It is instructive to compare this approach with other common ways of
thinking about the trend in a nonstationary series. If the $i$th observed
variable is a deterministic function of time plus a zero-mean stationary
ergodic process, $y_{it}={\Greekmath 010E} _{i}^{\det }(t)+c_{it}^{\det }$, we
typically could describe the deterministic time trend as the limit of a
forecast that would have been made in the arbitrarily distant past:
\begin{equation*}
{\Greekmath 010E} _{i}^{\det }(t)=\underset{h\rightarrow \infty }{\,\lim }\underset{
p\rightarrow \infty }{\,\lim }E(y_{it}\mid
y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1}).
\end{equation*}
By contrast, if the first difference of the $i$th variable is a zero-mean
stationary process, \cite{beveridge1981new} suggested that we think of the
trend as the forecast of the variable in the arbitrarily distant future.
They decomposed $y_{it}={\Greekmath 010E} _{it}^{BN}+c_{it}^{BN}$ where
\begin{equation*}
{\Greekmath 010E} _{it}^{BN}=\underset{h\rightarrow \infty }{\,\lim }\underset{
p\rightarrow \infty }{\,\lim }E(y_{i,t+h}\mid
y_{it},y_{i,t-1},...,y_{i,t-p+1}).
\end{equation*}
While these concepts of trend have some appeal, they have the significant
practical drawback that both are based on the properties of forecasts at
infinite horizons. They thus depend on conjectures of what happens at
infinity, conjectures that are impossible to verify on the basis of a finite
sample of observations.
By contrast, forecasts at a two-year horizon are something that can be
reasonably investigated without auxiliary assumptions in samples of typical
size. If we choose $h$ to correspond to a two-year horizon, \cite
{hamilton2018you} argued that the decomposition
\begin{equation}
y_{it}=\mathbb{P}(y_{it}|1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})+c_{it}
\label{eqn_hamilton_regression}
\end{equation}
can be viewed as a way to implement traditional approaches to trend-cycle
decomposition that is practical and robust. Moreover, the primary reason we
would go wrong in making a two-year-ahead forecast of most economic time
series is due to unforeseen cyclical changes. For example, the value of $
y_{it}$ will be significantly below {$\mathbb{P}(y_{it}|{z}_{i,t-h})$} if
the economy goes into a recession after $t-h$ and significantly above the
forecast if recovery from a downturn is more robust than expected. For this
reason, \cite{hamilton2018you} proposed to use $h$ corresponding to a
two-year horizon as a way to isolate the stationary cyclical component of
the $i$th variable, for example, $h=24$ for monthly data.
Another important benefit of setting $h=24$ is that it offers a solution to
the huge practical problem of dealing with outliers. To see why this is the
case, suppose for illustration that the $i$th variable is characterized by a
random walk: $y_{it}=y_{i,t-1}+{\Greekmath 0122} _{it}.$ For a monthly random
walk, the residual from the two-year population linear projection is given
by $c_{it}=\mathop{\textstyle \sum }\nolimits_{s=0}^{23}{\Greekmath 0122} _{i,t-s}.$ From the Central
Limit Theorem, $c_{it}$ typically exhibits much less kurtosis than the
one-period-ahead forecast error ${\Greekmath 0122} _{it}.$\footnote{\label{fn_CLT}
For example, if $c_{it}^{(h)}= {\Greekmath 0122} _{it}+\ldots + {\Greekmath 0122}
_{i,t-h+1}$ and ${\Greekmath 0122} _{it}$ is i.i.d. with $E({\Greekmath 0122} _{it})=0,$
$E({\Greekmath 0122} _{it}^{2})={\Greekmath 011B} ^{2},$ and $E({\Greekmath 0122} _{it}^{4})={\Greekmath 0114}
,$ then $E\left[ (c_{it}^{(h)})^{2}\right] =h{\Greekmath 011B} ^{2}$ and $E\left[
(c_{it}^{(h)})^{4}\right] =h{\Greekmath 0114} +3h(h-1){\Greekmath 011B} ^{4}.$ The excess kurtosis
of $c_{it}^{(h)}$ is then
\begin{equation*}
\frac{E\left[ (c_{it}^{(h)})^{4}\right] }{\left\{ E\left[ (c_{it}^{(h)})^{2}
\right] \right\} ^{2}}-3=\frac{h{\Greekmath 0114} +3h(h-1){\Greekmath 011B} ^{4}}{h^{2}{\Greekmath 011B} ^{4}}
-3=\frac{\left[ ({\Greekmath 0114} /{\Greekmath 011B} ^{4})-3\right] }{h}.
\end{equation*}
If ${\Greekmath 0122} _{it}$ has positive excess kurtosis (${\Greekmath 0114} /{\Greekmath 011B}
^{4}>3$), then the excess kurtosis of $c_{it}^{(h)}$ declines monotonically
to zero as $h$ increases. For a numerical illustration, suppose ${\Greekmath 0122}
_{it}$ comes from a mixture of two Normals, where ${\Greekmath 0122} _{it}\sim
N(0,{\Greekmath 011B} _{1}^{2})$ with probability $p$ and $N(0,{\Greekmath 011B} _{2}^{2})$ with
probability $1-p.$ Then $c_{it}^{(h)}=\mathop{\textstyle \sum }\nolimits_{s=0}^{h-1}{\Greekmath 0122}
_{i,t-s}$ is a mixture of $h+1$ Normals, distributed
\begin{equation*}
N(0,k{\Greekmath 011B} _{1}^{2}+(h-k){\Greekmath 011B} _{2}^{2})\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ with probability }\left(
\begin{array}{l}
h \\
k
\end{array}
\right) p^{k}(1-p)^{h-k}
\end{equation*}
and variance $E\left[ c_{it}^{(h)}\right] ^{2}=V^{(h)}=h[p{\Greekmath 011B}
_{1}^{2}+(1-p){\Greekmath 011B} _{2}^{2}].$ For ${\Greekmath 011B} _{1}=1,$ ${\Greekmath 011B} _{2}=5,$ and $
p=0.9,$ the probability that $c_{it}^{(h)}$ exceeds four of its standard
deviations, that is, Prob$\left( \left\vert c_{it}^{(h)}\right\vert >4\sqrt{
V^{(h)}}\right) $, is 1.4\% for $h=1,$ 0.1\% for $h=12,$ and 0.005\% for $
h=24.$} We find in our empirical application in Section \ref{sec_fred_md}
that this feature is extremely helpful in dealing with outliers such as
those seen during the COVID-19 pandemic in 2020.
Notwithstanding, these benefits do not come without a cost. The larger the
value of $h,$ the more serial correlation there will be in $c_{it},$ and the
more prone PCA can be in small samples to the spurious factor problem
identified by \cite{onatski2021spurious}. Our simulations in Section \ref
{sec_simulations} suggest that if we have 50 years of data, using
two-year-ahead forecast errors works pretty well. For shorter samples, the
spurious factor problem can be more serious, and our recommendation for
those applications is to look for common factors in the one-year-ahead ($
h=12)$ or one-month-ahead $(h=1)$ forecasts. In the latter case, additional
correction for outliers is likely needed.
To summarize, our procedure is to estimate the same regression for every
variable $y_{it},$ regardless of whether we think it is stationary and
without making any conjecture about the nature of any nonstationarity. To
allow for persistent seasonal components in $y_{it}$, we recommend choosing $
p$ to be the number of observations in a year. Our procedure estimates the
following regression by OLS for every variable,
\begin{equation}
y_{it}=k_{i}+{\Greekmath 010B} _{i1}y_{i,t-h}+{\Greekmath 010B} _{i2}y_{i,t-h-1}+\cdots +{\Greekmath 010B}
_{ip}y_{i,t-h-p+1}+c_{it}, \label{eqn_recommended_regression}
\end{equation}
with $h=8$ and $p=4$ for quarterly data and $h=24$ and $p=12$ for monthly
data. We will refer to the residual from the estimated regression $\hat{c}
_{it}$ as the OLS residual for variable $y_{it}$ and the residual from the
population linear projection $c_{it}$ as the true cyclical component. The
value $\hat{c}_{it}$ is a consistent estimate of $c_{it},$ and the true
value $c_{it}$ is stationary as long as any nonstationarity in $y_{it}$ is
characterized by either a polynomial time trend of order $d_{i}$ or an $
I(d_{i})$ process with $d_{i}\leq p.$ Our procedure is to perform PCA on the
regression residuals $\{\hat{c}_{1t},...,\hat{c}_{Nt}\}.$
One practical decision is whether a nonlinear transformation of the raw data
is necessary for $\Delta ^{d_{i}}y_{it}$ to be stationary for some $d_{i}$.
If taking first differences of the log is the correct way to produce a
stationary series, then taking the change in the level would not produce a
stationary series. We recommend using logs for variables like output or
prices which are usually described in terms of growth rates. \ For such
variables we use the log of the level of the variable, $y_{it}=\log Y_{it}$
as the variable in the regression (\ref{eqn_recommended_regression}). \ For
variables like interest rates or the unemployment rate that are already
quoted in percentage terms, we use the raw data $y_{it}=Y_{it}$ in the
regression.
The true cyclical component $c_{it}$ has mean zero and is stationary for a
wide range of processes. However, the population value of ${\Greekmath 010B} _{i0}$ is
not known but must be estimated by regression. In Section \ref{sec_c_obs} we
characterize our assumptions about the factor structure that we hypothesize
describes the true values of $c_{it},$ and use standard results to establish
that these population factors could be consistently estimated if the true
values of $c_{it}$ were observed without error. Section \ref
{sec_c_not_observed} considers the case when we do not know the value of $
d_{i}$ for each series, do not know whether it is stationary or
characterized by a deterministic time trend or an $I(d_{i})$ process, and
the $c_{it}$ are not observed. In that section we analyze the consequences
of performing PCA on the estimated OLS residuals $\hat{c}_{it}.$
\section{Principal component analysis when the true cyclical component is
observed}
\label{sec_c_obs}In the previous section we defined the true cyclical
component $c_{it}$ to be the residual from a population linear projection of
$y_{it}$ on $(1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})^{\prime }$, and
noted that $c_{it}$ is stationary for a broad class of possible processes.
In this section we provide sufficient conditions under which the true
cyclical components for a collection of $N$ different variables would have a
factor structure that could be consistently estimated using PCA if we
observed the true value of $c_{it}$ for each variable. The set-up and
results in this section closely follow \cite{stock2002forecasting}.
\subsection{Assumed factor structure of the true cyclical components}
Collect the true cyclical components for the $N$ different series at time $t$
in an $(N\times 1)$ vector $C_{t}=(c_{1t},...,c_{Nt})^{\prime }.$ We
postulate that these are characterized by a factor structure of the form
\begin{equation}
\underset{(N\times 1)}{C_{t}}=\underset{(N\times r)}{\Lambda }\underset{
(r\times 1)}{F_{t}}+\underset{(N\times 1)}{e_{t}}.
\label{eqn_factor_equation}
\end{equation}
The number of latent factors $r$ is much less than the number of variables $
N,$ but the $r$ factors are assumed to account for most of the variance of $
C_{t}$ in a sense made formal below. Since the factors are unobserved, $
C_{t}=\Lambda H^{-1}HF_{t}+e_{t}$ would imply the identical observable model
as (\ref{eqn_factor_equation}). Thus some normalizations are necessary in
order to talk about consistently estimating the $j$th factor $f_{jt}.$ In
empirical estimation, practitioners typically resolve this ambiguity by
estimating the $j$th column of $\Lambda $ by the eigenvector associated with
the $j$th largest eigenvalue of $T^{-1}\sum\nolimits_{t=1}^{T}C_{t}C_{t}^{
\prime }.$ Note that such a procedure implies a normalization in which the
columns of $\Lambda $ are orthogonal to each other and the elements of $
F_{t} $ are uncorrelated with each other and ordered by the size of their
variance. We follow \cite{stock2002forecasting} in how to characterize these
conventions as the cross-section dimension $N$ and time-series dimension $T$
get large.\footnote{
See \cite{BaiNg2013} and \cite{StockWatson2016Handbook} for discussion of
alternative normalizations.}
\begin{assumption}[factor structure]
\label{ass_factor} \mbox{}
\begin{enumerate}
\item \label{ass_factor_a} $(\Lambda^{\prime }\Lambda/N) \rightarrow I_r$.
\item \label{ass_factor_b} $E[F_tF_t^{\prime }] = \Omega_{FF}$, where $
\Omega_{FF}$ is a diagonal matrix with ${\Greekmath 0121}_{ii} > {\Greekmath 0121}_{jj} > 0 $ for $
i < j$.
\item \label{ass_factor_c}$|{\Greekmath 0115} _{ij}|\leq \bar{{\Greekmath 0115}}<\infty $.
\item \label{ass_factor_d}$T^{-1}\sum_{t}F_{t}F_{t}^{\prime }\overset{p}{
\rightarrow }\Omega _{FF}$.
\end{enumerate}
\end{assumption}
In addition to implementing the property that eigenvectors of a symmetric
matrix are orthogonal, Assumption \ref{ass_factor}(\textit{\ref{ass_factor_a}
}) requires that each factor makes a nonnegligible contribution to the
average variance of $c_{it}$ across $i$. That is, if we were to imagine
adding more variables (increasing $N)$ with ${\Greekmath 0115} _{ij}=0$ for all $i$
greater than some fixed $N_{0},$ then Assumption \ref{ass_factor}(\textit{
\ref{ass_factor_a}}) could not hold. Likewise \ref{ass_factor}(\textit{\ref
{ass_factor_b}}) and \ref{ass_factor}(\textit{\ref{ass_factor_d}}) require
that each factor continues to matter as the number of time-series
observations $T$ grows. These conditions are consistent with serial
dependence of the factors, but rely on the fact that $C_{t}$ is stationary.
Let ${\Greekmath 010D} $ denote an $(N\times 1)$ vector and $\Gamma =\{{\Greekmath 010D} :{\Greekmath 010D}
^{\prime }{\Greekmath 010D} /N=1\}.$ Note that if ${\Greekmath 010D} $ were the $j$th column of $
\Lambda ,$ the scalar ${\Greekmath 010D} ^{\prime }\Lambda F_{t}/N$ would converge to $
f_{jt}$ and $(N^{2}T)^{-1}\sum\nolimits_{t=1}^{T}{\Greekmath 010D} ^{\prime }\Lambda
F_{t}F_{t}^{\prime }\Lambda ^{\prime }{\Greekmath 010D} \overset{p}{\rightarrow }{\Greekmath 0121}
_{jj}.$ The assumption that the idiosyncratic elements $e_{t}$ do not have a
factor structure requires that there is no value of ${\Greekmath 010D} $ for which the
analogous operation applied to $e_{t}$ would lead to anything other than
zero: $\underset{{\Greekmath 010D} \in \Gamma }{\sup }(N^{2}T)^{-1}\sum
\nolimits_{t=1}^{T}{\Greekmath 010D} ^{\prime }e_{t}e_{t}^{\prime }{\Greekmath 010D} \overset{p}{
\rightarrow }0.$ \cite{stock2002forecasting} used the following assumptions
to guarantee the absence of a factor structure in $e_{t}.$
\begin{assumption}[moments of the errors]
\label{ass_error} \mbox{}
\begin{enumerate}
\item $\underset{N\rightarrow \infty}{\lim}\sup_t
\sum_{s=-\infty}^{\infty}|E[e_t^{\prime }e_{t+s}/N]|<\infty$.
\item $\underset{N\rightarrow \infty}{\lim}\sup_t
N^{-1}\sum_{i=1}^{N}\sum_{j=1}^{N}|E[e_{it}e_{jt}]|<\infty$, where $e_{it}$
denotes the $i$th element of $e_t$.
\item $\underset{N\rightarrow \infty }{\lim }\sup_{t,s}N^{-1}\sum_{i=1}^{N}
\sum_{j=1}^{N}|cov[e_{is}e_{it},e_{js}e_{jt}]|<\infty $.
\end{enumerate}
\end{assumption}
Some might be concerned that we have simply postulated that the true
cyclical components are characterized by Assumptions \ref{ass_factor} and
\ref{ass_error}. But something very similar is done in traditional
applications that assume conditions like these characterize specified
stationary transformations of the original data. Indeed, insofar as the
cyclical components have a common primitive definition in terms of $h$
-period-ahead forecast errors, we would argue that these assumptions are
easier to defend in our application than in many others.
\subsection{Consequences of applying PCA to the true cyclical components}
Recall that the $(N\times 1)$ vector of true cyclical components $C_{t}$ is
stationary and has population mean zero. If $C_{t}$ was observed directly,
its estimated sample variance matrix would be $S=T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}C_{t}C_{t}^{\prime }$ and a linear combination ${\Greekmath 010D}
^{\prime }C_{t}$ for any $(N\times 1)$ vector ${\Greekmath 010D} $ would have sample
variance ${\Greekmath 010D} ^{\prime }S{\Greekmath 010D} .$ If $C_{t}$ was observed, the first
estimated principal component (denoted $\tilde{f}_{1t}=N^{-1}\tilde{{\Greekmath 0115}}
_{1}^{\prime }C_{t})$ would be defined as the linear combination that has
maximum sample variance subject to a normalization condition such as ${\Greekmath 010D}
\in \Gamma =\{{\Greekmath 010D} :{\Greekmath 010D} ^{\prime }{\Greekmath 010D} /N=1\}$:
\begin{equation}
\tilde{{\Greekmath 0115}}_{1}=\arg \underset{{\Greekmath 010D} \in \Gamma }{\sup }\tilde{R}
({\Greekmath 010D} ) \label{eqn def lambda1 tilde}
\end{equation}
\begin{equation}
\tilde{R}({\Greekmath 010D} )=(N^{2}T)^{-1}{\Greekmath 010D} ^{\prime
}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}C_{t}C_{t}^{\prime }{\Greekmath 010D} . \label{eqn_Rtilde}
\end{equation}
Note we are normalizing $\tilde{{\Greekmath 0115}}_{1}^{\prime }\tilde{{\Greekmath 0115}}
_{1}/N=1 $ as we did asymptotically for the columns of $\Lambda $ in
Assumption \ref{ass_factor}(\ref{ass_factor_a})$.$ We also divide the sample
variance of ${\Greekmath 010D} ^{\prime }C_{t}$ by $N^{2}$ in anticipation of the
result that $\tilde{R}(\tilde{{\Greekmath 0115}}_{1})$ will converge to a fixed
constant as $N$ and $T$ grow. The solution to (\ref{eqn def lambda1 tilde})
is obtained by setting $\tilde{{\Greekmath 0115}}_{1}$ proportional to the eigenvector
of $S=T^{-1}\sum\nolimits_{t=1}^{T}C_{t}C_{t}^{\prime }$ associated with the
largest eigenvalue. For example, if we calculated eigenvectors of this
matrix using code that normalizes eigenvectors to have unit length and
orders eigenvalues by decreasing size, $\tilde{{\Greekmath 0115}}_{1}$ would be $\sqrt{
N}$ times the first eigenvector. The largest eigenvalue of $S$ is equal to $
T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\tilde{f}_{1t}^{2},$ the sample variance of
the first principal component$.$ The $j$th principal component $N^{-1}\tilde{
{\Greekmath 0115}}_{j}^{\prime }C_{t}$ is found by maximizing $\tilde{R}({\Greekmath 010D} )$
subject to the constraint that ${\Greekmath 010D} $ is orthogonal to $\tilde{{\Greekmath 0115}}
_{1},...,\tilde{{\Greekmath 0115}}_{j-1}.$ The solution for $\tilde{{\Greekmath 0115}}_{j}$ is
proportional to the eigenvector of $S$ associated with the $j$th largest
eigenvalue.
Alternatively, if we observed the true factors $F_{t}$ and loadings $\Lambda
,$ we could calculate the component of the variance of ${\Greekmath 010D} ^{\prime
}C_{t}$ that is attributable to the $r$ factors alone:
\begin{equation}
R^{\ast }({\Greekmath 010D} )=(N^{2}T)^{-1}{\Greekmath 010D} ^{\prime
}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\Lambda F_{t}F_{t}^{\prime }\Lambda ^{\prime
}{\Greekmath 010D} . \label{eqn_R*}
\end{equation}
\cite{stock2002forecasting} showed that under Assumptions \ref{ass_factor}
and \ref{ass_error}, the maximum value for (\ref{eqn_Rtilde}) (which is
given by the largest eigenvalue of $S)$ and the supremum of (\ref{eqn_R*})
over all ${\Greekmath 010D} \in \Gamma $ converge in probability to the same number $
{\Greekmath 0121} _{11},$ which is the population variance of the first factor, and
that $\tilde{{\Greekmath 0115}}_{1}^{\prime }C_{t}/N$ gives a consistent estimate of $
f_{1t}$ up to a sign. If we were to estimate $k>r$ principal components, the
first $r$ would consistently estimate $f_{jt}$ up to a sign normalization
and the last $k-r$ would asymptotically have zero variance. Their results
are a special case of Theorem \ref{thm_convergence_general_M} developed in
the following section.
\section{Principal component analysis when the cyclical component must be
estimated}
\label{sec_c_not_observed}
In this section we assume that we do not observe the true cyclical component
$c_{it}$ of series $i$ but have an estimate $\hat{c}_{it}=c_{it}+\hat{v}
_{it}.$ Let $\hat{C}_{t}=(\hat{c}_{1t},...,\hat{c}_{Nt})^{\prime }$ and $
\hat{V}_{t}=(\hat{v}_{1t},...,\hat{v}_{Nt})^{\prime }.$ We investigate the
properties of principal components calculated from the estimated cyclical
components:
\begin{equation}
\hat{f}_{jt}=N^{-1}\hat{{\Greekmath 0115}}_{j}^{\prime }\hat{C}_{t}
\label{eqn_HX_f1_definition}
\end{equation}
\begin{equation*}
\hat{{\Greekmath 0115}}_{j}=\arg \underset{\{{\Greekmath 010D} \in \Gamma ,{\Greekmath 010D} ^{\prime }\hat{
{\Greekmath 0115}}_{1}=\cdots ={\Greekmath 010D} ^{\prime }\hat{{\Greekmath 0115}}_{j-1}=0\}}{\sup }\hat{R}
({\Greekmath 010D} )
\end{equation*}
\begin{equation*}
\hat{R}({\Greekmath 010D} )=(N^{2}T)^{-1}{\Greekmath 010D} ^{\prime }\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{C
}_{t}\hat{C}_{t}^{\prime }{\Greekmath 010D} .
\end{equation*}
We first state high-level sufficient conditions under which PCA applied to
the estimated cyclical components $\hat{C}_{t}$ gives consistent estimates
of the true factors $F_{t}.$ Let $\hat{v}_{it}=\hat{c}_{it}-c_{it}$ denote
the difference between the estimated and true cyclical component of series $
i $ at date $t$. {\ The conditions require }$N$ to grow with $T$ so as to
ensure that $N^{-1}\mathop{\textstyle \sum }\nolimits_{i=1}^{N}\hat{v}_{it}^{2}$ and $
T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}$ are $o_{p}(1).$
\begin{assumption}[high-level conditions on $\hat{v}_{it}$]
\label{ass_convergence_of_sum_vit2} The cross-section dimension $N$ grows
with the number of time series observations $T$ according to a function $
N(T) $ such that for all ${\Greekmath 010E} ,{\Greekmath 0122} >0$ there exists a $
T_{0}({\Greekmath 010E} ,{\Greekmath 0122} )$ such that for all $T>T_{0}({\Greekmath 010E} ,{\Greekmath 0122}
) $:
\begin{equation}
\mathnormal{Prob}\left\{ \frac{1}{N(T)}\sum_{i=1}^{N(T)}\hat{v}
_{it}^{2}>{\Greekmath 010E} \right\} <{\Greekmath 0122} \@ifnextchar*{\@tagstar}{\@tag}{i}
\label{eqn_high_assumption_i}
\end{equation}
\begin{equation}
\mathnormal{Prob}\left\{ \max_{1\leq i\leq N(T)}\frac{1}{T}\sum_{t=1}^{T}
\hat{v}_{it}^{2}>{\Greekmath 010E} \right\} <{\Greekmath 0122} . \@ifnextchar*{\@tagstar}{\@tag}{ii}
\label{eqn_high_assumption_ii}
\end{equation}
\end{assumption}
The following result establishes that if the error in estimating the
cyclical component satisfies Assumption \ref{ass_convergence_of_sum_vit2},
then the results $(\hat{f}_{jt},\hat{{\Greekmath 0115}}_{j},\hat{R}({\Greekmath 010D} ))$ of
applying PCA to the estimated cyclical components $\hat{C}_{t}$ give
consistent estimates of the magnitudes that characterize the true cyclical
components $C_{t}.$ For a proof see online Appendix A.
\begin{theorem}
\label{thm_convergence_general_M} Suppose that $C_{t}$ and $e_{t}$ in
equation (\ref{eqn_factor_equation}) satisfy Assumptions \ref{ass_factor}
and \ref{ass_error}. Let $\hat{C}_{t}=C_{t}+\hat{V}_{t}$ for $\hat{V}_{t}=(
\hat{v}_{1t},...,\hat{v}_{Nt})^{\prime }$ where $\hat{v}_{it}$ satisfy
Assumption \ref{ass_convergence_of_sum_vit2}. Let $\hat{f}_{1t},...,\hat{f}
_{kt}$ denote the first $k$ estimated principal components of $\hat{C}_{t}$ (
$\hat{f}_{jt}=\hat{{\Greekmath 0115}}_{j}^{\prime }\hat{C}_{t}/N)$ with $k\geq r$ and
let $\underset{(r\times 1)}{\hat{F}_{t}}=(\hat{f}_{1t},...,\hat{f}
_{rt})^{\prime }$ and\ $\underset{(N\times r)}{\hat{\Lambda}}=\left[
\begin{array}{ccc}
\hat{{\Greekmath 0115}}_{1} & \cdots & \hat{{\Greekmath 0115}}_{r}
\end{array}
\right] .$ Then
\begin{enumerate}
\item $T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{f}_{jt}^{2}=\hat{R}(\hat{{\Greekmath 0115}}
_{j})\overset{p}{\rightarrow }{\Greekmath 0121} _{jj}$ for $j=1,...,r;$
\item $T^{-1}\sum\nolimits_{t=1}^{T}\hat{f}_{jt}^{2}\overset{p}{\rightarrow }
0$ for $j=r+1,...,k;$
\item $\hat{\Xi}\hat{\Lambda}^{\prime }\Lambda /N\overset{p}{\rightarrow }
I_{r}$ where $\hat{\Xi}$ is a diagonal matrix whose row $j$ column $j$
element is $+1$ if $\hat{{\Greekmath 0115}}_{j}^{\prime }{\Greekmath 0115} _{j}>0$ and $-1$ if $
\hat{{\Greekmath 0115}}_{j}^{\prime }{\Greekmath 0115} _{j}<0;$
\item $\hat{\Xi}\hat{F}_{t}-F_{t}\overset{p}{\rightarrow }0$.
\end{enumerate}
\end{theorem}
Before presenting formal sufficient conditions for verifying Assumption \ref
{ass_convergence_of_sum_vit2}, we first discuss the intuition for why we
might expect it to hold. For $z_{i,t-h}=(1,y_{i,t-h},y_{i,t-h-1},...,
\allowbreak y_{i,t-h-p+1})^{\prime },$ the true cyclical component $c_{it}$
is the residual from a population linear projection of $y_{it}$ on $
z_{i,t-h} $ and $\hat{c}_{it}$ is the residual from the corresponding
estimated regression:
\begin{equation*}
c_{it}=y_{it}-{\Greekmath 010B} _{i0}^{\prime }z_{i,t-h}
\end{equation*}
\begin{equation*}
\hat{c}_{it}=y_{it}-\hat{{\Greekmath 010B}}_{i}^{\prime }z_{i,t-h}
\end{equation*}
\begin{equation*}
\hat{v}_{it}=\hat{c}_{it}-c_{it}=({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}_{i})^{\prime
}z_{i,t-h}
\end{equation*}
\begin{equation}
({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}_{i})=-\left(
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime }\right) ^{-1}\left(
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}c_{it}\right)
\label{eqn_alpha_minus_alphahat}
\end{equation}
\begin{equation}
\hat{v}_{it}^{2}=({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}_{i})^{\prime
}z_{i,t-h}z_{i,t-h}^{\prime }({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}_{i})
\label{eqn_vit2}
\end{equation}
\begin{equation}
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}=({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}
_{i})^{\prime }\left( \mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime
}\right) ({\Greekmath 010B} _{i0}-\hat{{\Greekmath 010B}}_{i}). \label{eqn_sum_vit2}
\end{equation}
Expression (\ref{eqn_sum_vit2}) will be recognized as the OLS Wald statistic
for testing the true null hypothesis $H_{0}:{\Greekmath 010B} _{i}={\Greekmath 010B} _{i0}$
multiplied by $\hat{{\Greekmath 011B}}_{i}^{2}=(T-k)^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{c}
_{it}^{2},$ the average squared regression residual. As we now sketch, this
statistic would be expected to be $O_{p}(1)$ for a wide class of stationary
and nonstationary processes. This would lead us to expect that the
individual terms $\hat{v}_{it}^{2}$ should be $o_{p}(1),$ with some
uniformity conditions across $i$ then guaranteeing Assumption \ref
{ass_convergence_of_sum_vit2}(\ref{eqn_high_assumption_i}). Likewise if we
divide (\ref{eqn_sum_vit2}) by $T$ we should again get an $o_{p}(1)$ random
variable, as required by Assumption \ref{ass_convergence_of_sum_vit2}(\ref
{eqn_high_assumption_ii}).
We now explore the intuition for why we would typically expect (\ref
{eqn_sum_vit2}) to be $O_{p}(1).$ Note that the expression can be written
\begin{eqnarray}
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2} &=&\left(
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}c_{it}z_{i,t-h}^{\prime }\right) \left(
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime }\right) ^{-1}\left(
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}c_{it}\right) \notag \\
&=&\left( \mathop{\textstyle \sum }\nolimits_{t=1}^{T}c_{it}z_{i,t-h}^{\prime }\Upsilon
_{iT}^{-1}\right) \left( \Upsilon
_{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime }\Upsilon
_{iT}^{-1}\right) ^{-1}\left( \Upsilon
_{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}c_{it}\right)
\label{eqn_sum_vit2_scaled}
\end{eqnarray}
where $\Upsilon _{iT}^{-1}$ could be any sequence of nonsingular matrices.
Note that the researcher does not need to know the value of $\Upsilon _{iT}.$
This matrix is just a device to calculate the asymptotic properties of the
left side of (\ref{eqn_sum_vit2_scaled}) (which does not depend on $\Upsilon
_{iT})$ under various possible forms of nonstationarity for $z_{it}.$
If $y_{it}$ were stationary, we would use $\Upsilon
_{iT}^{-1}=T^{-1/2}I_{p+1}$ to find the asymptotic distribution of $
\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}.$ In this case, $z_{i,t-h}c_{it}$
is a stationary random variable which by the definition of the true cyclical
component $c_{it}$ has population mean zero. In this case, a central limit
theorem could be used to establish that $\Upsilon
_{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}c_{it}=T^{-1/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}z_{i,t-h}c_{it}$ converges to a Normal distribution.
Likewise in the stationary case,
\begin{equation*}
\Upsilon _{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime
}\Upsilon
_{iT}^{-1}=T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime }\overset
{p}{\rightarrow }E(z_{i,t-h}z_{i,t-h}^{\prime }).
\end{equation*}
For this reason, $\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}$ would typically
be $O_{p}(1)$ when $y_{it}$ is stationary.
Consider next the case when $d_{i}=1$ and $\Delta y_{it}$ (denoted below by $
u_{it})$ is a zero-mean stationary process. Taking a regression with $p=3$
lags for illustration, an OLS regression of $y_{it}$ on $
z_{i,t-h}=(1,y_{i,t-h},y_{i,t-h-1},y_{i,t-h-2})^{\prime }$ has the identical
fitted values as an OLS regression of $y_{it}$ on $
R_{i}^{-1}z_{i,t-h}=(u_{i,t-h},u_{i,t-h-1},1,y_{i,t-h})^{\prime }.$ In this
case to calculate the asymptotic distribution of (\ref{eqn_sum_vit2_scaled})
we would take
\begin{equation*}
\Upsilon _{iT}^{-1}=\left[
\begin{array}{ll}
T^{-1/2}I_{3} & 0 \\
0 & T^{-1}
\end{array}
\right] R_{i}^{-1}
\end{equation*}
\begin{eqnarray*}
&&\Upsilon _{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime
}\Upsilon _{iT}^{-1} \\
&=&\left[
\begin{array}{llll}
T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}u_{i,t-h}^{2} & T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h}u_{i,t-h-1} & T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h} & T^{-3/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h}y_{i,t-h} \\
T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}u_{i,t-h-1}u_{i,t-h} & T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h-1}^{2} & T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h-1} & T^{-3/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h-1}y_{i,t-h} \\
T^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}u_{i,t-h} & T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}u_{i,t-h-1} & 1 & T^{-3/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}y_{i,t-h} \\
T^{-3/2}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}y_{i,t-h}u_{i,t-h} & T^{-3/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}y_{i,t-h}u_{i,t-h-1} & T^{-3/2}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}y_{i,t-h} & T^{-2}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}y_{i,t-h}^{2}
\end{array}
\right] .
\end{eqnarray*}
Under standard unit-root asymptotics (e.g., \citet[p. 506]{hamilton1994}),
we expect $\Upsilon
_{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}z_{i,t-h}^{\prime }\Upsilon
_{iT}^{-1}\overset{d}{\rightarrow }Q_{i}$ with
\begin{equation}
Q_{i}=\left[
\begin{array}{llll}
{\Greekmath 010D} _{i0} & {\Greekmath 010D} _{i1} & 0 & 0 \\
{\Greekmath 010D} _{i1} & {\Greekmath 010D} _{i0} & 0 & 0 \\
0 & 0 & 1 & {\Greekmath 0121} _{i}\int_{0}^{1}W_{i}(r)dr \\
0 & 0 & {\Greekmath 0121} _{i}\int_{0}^{1}W_{i}(r)dr & {\Greekmath 0121} _{i}^{2}\int_{0}^{1}\left[
W_{i}(r)\right] ^{2}dr
\end{array}
\right] \label{eqn_Q_example}
\end{equation}
for ${\Greekmath 010D} _{ij}=E(u_{it}u_{i,t-j}),$ ${\Greekmath 0121}
_{i}^{2}=\mathop{\textstyle \sum }\nolimits_{j=-\infty }^{\infty }{\Greekmath 010D} _{ij},$ and $W_{i}(r)$
standard Brownian motion. Similar tools can be used to establish that $
\Upsilon _{iT}^{-1}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}z_{i,t-h}c_{it}\overset{d}{
\rightarrow }q_{i}$ and $\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}\overset{d}{
\rightarrow }q_{i}^{\prime }Q_{i}^{-1}q_{i}\sim O_{p}(1).$
If instead $E(\Delta y_{it})={\Greekmath 0116} _{i}\neq 0$, the level $y_{it}$ would be
dominated asymptotically by a deterministic time trend ${\Greekmath 0116} _{i}t.$\footnote{
\label{fn_trend} That is, for any nonnegative integer ${\Greekmath 0117} $,
\begin{equation*}
T^{-({\Greekmath 0117} +1)}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}y_{it}^{{\Greekmath 0117} }-T^{-({\Greekmath 0117}
+1)}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}({\Greekmath 0116} _{i}t)^{{\Greekmath 0117} }\overset{p}{\rightarrow }0
\end{equation*}
as $T\rightarrow \infty $ and $T^{-({\Greekmath 0117} +1)}\mathop{\textstyle \sum }\nolimits_{t=1}^{T}({\Greekmath 0116}
_{i}t)^{{\Greekmath 0117} }\rightarrow {\Greekmath 0116} _{i}^{{\Greekmath 0117} }/({\Greekmath 0117} +1).$} In this case (and for
general $p\geq 2)$ we would use\ $R_{i}^{-1}z_{i,t-h}=(\Delta y_{i,t-h}-{\Greekmath 0116}
_{i},\Delta y_{i,t-h-1}-{\Greekmath 0116} _{i},...,\Delta y_{i,t-h+2}-{\Greekmath 0116}
_{i},1,y_{i,t-h})^{\prime }$ and
\begin{equation*}
\Upsilon _{iT}^{-1}=\left[
\begin{array}{ll}
T^{-1/2}I_{p} & 0 \\
0 & T^{-3/2}
\end{array}
\right] R_{i}^{-1}
\end{equation*}
to establish that $\mathop{\textstyle \sum }\nolimits_{t=1}^{T}\hat{v}_{it}^{2}\sim O_{p}(1).$
General cases for $d_{i}\leq p$ are examined in \cite{hamilton2018you}. Note
that we do not need to know the value of $d_{i}$ or $\Upsilon _{iT}$ to
estimate any magnitudes -- in every case we are talking about a regression
of $y_{it}$ on the lagged levels $
z_{i,t-h}=(1,y_{i,t-h},y_{i,t-h-1},...,y_{i,t-h-p+1})^{\prime }.$ Instead, $
\Upsilon _{iT}$ is just a device to analyze the properties if the levels
regression is applied to a variety of different stationary and nonstationary
processes.
We are now in a position to state formally sufficient conditions under which
Assumption \ref{ass_convergence_of_sum_vit2} would be satisfied, for which
we let $\left\Vert X\right\Vert $ denote the Frobenius/Euclidean norm (the
square root of the sums of squares of the elements of a vector or matrix $X)$
.\footnote{
An earlier version of this paper avoided the requirement that $T$ grows
faster than $N$ by assuming that the true cyclical component $c_{it}$ is
uniformly bounded.}
\begin{assumption}[sufficient conditions for Assumption \protect\ref
{ass_convergence_of_sum_vit2}]
\label{ass_sufficient} For each $i$, there exist some nonsingular matrix $
\Upsilon_{iT}$ and a possibly random matrix $Q_i$, such that \mbox{}
\begin{enumerate}
\item \label{ass_sufficient_a} For some positive constant $\mathfrak{c}_{1}$
,
\begin{align*}
\max_{1\leq i\leq N}\max_{1\leq t\leq T}E( \Vert \sqrt{T}
\Upsilon_{iT}^{-1}z_{i,t-h}\Vert^2 ) \leq \mathfrak{c}_1\quad \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ and}
\quad \max_{1\leq i\leq N}E( \Vert\sum_{t=1}^T
\Upsilon_{iT}^{-1}z_{i,t-h}c_{it}\Vert^2 ) \leq \mathfrak{c}_1.
\end{align*}
\item \label{ass_sufficient_b} For some positive constants $\mathfrak{c}_2$,
$\mathfrak{c}_3$, $\mathfrak{c}_4$, and $\mathfrak{c}_5$,
\begin{align*}
\max_{1\leq i\leq N}\mathnormal{Prob}\big({\Greekmath 0125}_{\min}(Q_i) < {\Greekmath 010E}\big)
\leq \mathfrak{c}_2 \exp\big\{-\mathfrak{c}_3\big(\frac{1}{{\Greekmath 010E}}\big)^{1/
\mathfrak{c}_4}\big\}\quad \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ for all ${\Greekmath 010E}\leq \mathfrak{c}_5$},
\end{align*}
where ${\Greekmath 0125}_{\min}$ denotes the minimum eigenvalue. In addition, $N\log^{2
\mathfrak{c}_4}(N)/T \rightarrow 0$.
\item \label{ass_sufficient_c} The cross-section dimension $N$ grows with
the number of time series observations $T$ according to a function $N(T) $,
such that
\begin{align*}
\mathnormal{Prob}\big(\max_{1\leq i\leq N} \big\Vert \sum_{t=1}^T
\Upsilon_{iT}^{-1}z_{i,t-h}z_{i,t-h}^{\prime }\Upsilon_{iT}^{-1} - Q_i
\big\Vert > \mathfrak{c}_{6}(N)\big) \rightarrow 0
\end{align*}
for some positive sequence $\mathfrak{c}_{6}(N)= o(\log^{-\mathfrak{c}
_4}(N)) $.
\end{enumerate}
\end{assumption}
{Condition \ref{ass_sufficient}(\ref{ass_sufficient_a}) is quite mild, as it
only requires finite moments on the transformed regressors and the cyclical
components after suitable normalization. For example, when }$\Delta y_{it}$
is a zero-mean stationary process, the first inequality requires $\Delta
y_{it}$ and $y_{it}/\sqrt{T}$ to have finite variance. {Conditions \ref
{ass_sufficient}(\ref{ass_sufficient_b}--\ref{ass_sufficient_c}) are key to
our theoretical developments. Specifically, they posit that the matrix that
gets inverted in (\ref{eqn_sum_vit2_scaled}) can be uniformly approximated
by a possibly random limit $Q_{i}$, and that the minimum eigenvalue of the
limiting matrix exhibits an exponential-type tail at zero. }
\begin{theorem}
\label{thm_main_result} Assumption \ref{ass_sufficient} implies Assumption
\ref{ass_convergence_of_sum_vit2}.
\end{theorem}
{Next we relate Assumption \ref{ass_sufficient} to conditions for general
linear processes.}
\begin{assumption}[linear process]
\label{ass_linear_process} There exist positive constants $\mathfrak{c}_{7}$
--$\mathfrak{c}_{11}$, such that for each $i\in \{1,2,\dots ,N\}$ there
exist a ${\Greekmath 0116} _{i}$ and an integer $d_{i}$ for which $
u_{it}=(1-L)^{d_{i}}y_{it}-{\Greekmath 0116} _{i}=\sum_{l=0}^{\infty }{\Greekmath 0120} _{il}{\Greekmath 0111}
_{i,t-l}$, and
\begin{enumerate}
\item \label{ass_linear_process_a} $\max_{1\leq i\leq N}|{\Greekmath 0120} _{il}|\leq
\mathfrak{c}_{7}l^{-\mathfrak{c}_{8}}$ with $\mathfrak{c}_{8}>2$.
\item \label{ass_linear_process_a1} $\min_{1\leq i\leq N}|\sum_{l=0}^{\infty
}{\Greekmath 0120}_{il}|\geq \mathfrak{c}_{9}>0$.
\item \label{ass_linear_process_b} ${\Greekmath 0111} _{it}$ is an i.i.d. process with
zero mean, unit variance, and $\max_{1\leq i\leq N}E({\Greekmath 0111} _{it}^{4})\leq
\mathfrak{c}_{10}$.
\item \label{ass_linear_process_d} Let $\tilde{Q}_{i}$ be the covariance
matrix of $(u_{it},\ u_{i,t-1},\ \dots ,\ u_{i,t-p+d_{i}+1})$. Then ${\Greekmath 0125}
_{\min }(\tilde{Q}_{i})\geq \mathfrak{c}_{11}>0$.
\end{enumerate}
\end{assumption}
Conditions \ref{ass_linear_process}(\ref{ass_linear_process_a}) and \ref
{ass_linear_process}(\ref{ass_linear_process_a1}) would be satisfied if $
u_{it}$ is any zero-mean stationary $AR(g)$ process for some $g$,
\begin{equation*}
u_{it}={\Greekmath 011E} _{i1}u_{i,t-1}+\cdots +{\Greekmath 011E} _{ig}u_{i,t-g}+{\Greekmath 0111} _{it}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ \
with }\left\Vert z\right\Vert >1\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ for any }z\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ satisfying\ }
\left\vert 1-{\Greekmath 011E} _{i1}z-\cdots -{\Greekmath 011E} _{ig}z^{g}\right\vert =0.
\end{equation*}
Alternatively, \ref{ass_linear_process}(\ref{ass_linear_process_a}) and \ref
{ass_linear_process}(\ref{ass_linear_process_a1}) would also be satisfied
for any $MA(g)$ process that has a representation that is not arbitrarily
close to being noninvertible,
\begin{equation*}
u_{it}={\Greekmath 0111} _{it}+{\Greekmath 0112} _{i1}{\Greekmath 0111} _{i,t-1}+\cdots +{\Greekmath 0112} _{ig}{\Greekmath 0111} _{i,t-g}
\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ \ with }\left\vert 1+{\Greekmath 0112} _{i1}+\cdots +{\Greekmath 0112} _{ig}\right\vert >
\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{$\mathfrak{c}$}_{9}.
\end{equation*}
Condition \ref{ass_linear_process} (\ref{ass_linear_process_a1}) defines the
necessary degree of differencing $d_{i}.$ For example, if $y_{it}={\Greekmath 0120}
_{i}(L){\Greekmath 0111} _{it}$ with ${\Greekmath 0120} _{i\ell }$ satisfying ${\Greekmath 0120} _{i}(1)>\mathfrak{c
}_{9},$ then the $MA(\infty )$ coefficients $\tilde{{\Greekmath 0120}}_{i\ell }$ of the
first difference $(1-L)y_{it}=(1-L){\Greekmath 0120} _{i}(L){\Greekmath 0111} _{it}=\tilde{{\Greekmath 0120}}
_{i}(L){\Greekmath 0111} _{it}$ would not satisfy (\ref{ass_linear_process_a1}) because $
\tilde{{\Greekmath 0120}}_{i}(1)=(1-1){\Greekmath 0120} _{i}(1)=0.$
\begin{theorem}
\label{thm_main_result_linear_process} Suppose that $p\geq 1$ and Assumption
\ref{ass_linear_process} is satisfied for every $i\in\{1,2,\dots,N\}$ for
some $d_i\in\{0,1\}$, and that the true cyclical components $c_{it}$ satisfy
Assumptions \ref{ass_factor} and \ref{ass_error}. If $N \log^{12}(N) / T
\rightarrow 0$, then results (i)--(iv) of Theorem \ref
{thm_convergence_general_M} hold.
\end{theorem}
We now briefly outline our proof technique and main technical contribution.
It is helpful to compare the current setting with the case where $N$ is
fixed and $T$ grows. The fixed-$N$ case is straightforward, as one can
establish either convergence in probability to a nonrandom $Q_i$ (under weak
stationarity) or convergence in distribution to a random $Q_i$ (as
functionals of Brownian motion). This aligns with the reasoning provided
earlier (see \eqref{eqn_Q_example} and the related discussion). However,
when $N$ is allowed to grow, the theoretical analysis requires uniform
control of the estimation error across $i=1, \dots, N$, as implied by
Assumption \ref{ass_convergence_of_sum_vit2}. The core of our proof is to
decompose the quantity of interest into a sum of i.i.d. components plus
negligible remainder terms. We then apply the coupling method of \cite
{komlos1976approximation} to establish a distributional approximation that
is uniformly valid for $i=1,\dots,N$.
While our proof of Theorem \ref{thm_main_result_linear_process} only covers
the case when $d_{i}\leq 1,$ we conjecture that a similar result could be
obtained for any $d_{i}\leq p.$ Extending the proof to this case requires
verifying the three conditions in Assumption \ref{ass_sufficient}.
Assumption \ref{ass_sufficient}(\ref{ass_sufficient_a}) requires only finite
moments of the regressors and the cyclical components after suitable
transformation and rescaling by the appropriate $\Upsilon _{iT}$. This is
true for processes satisfying Assumption \ref{ass_linear_process} with $
d_{i}=2,$ although the calculations demonstrating this would be more tedious
than the $d_{i}=1$ case we review. It is also possible to establish \ref
{ass_sufficient}(\ref{ass_sufficient_c}) for $d_{i}\geq 2$. The key step in
our proof builds on the coupling result of \cite{komlos1976approximation},
which provides a uniform (in both $1\leq i\leq N$ and $1\leq t\leq T$)
strong approximation of a random walk by sums of independent Gaussian
variables. This result can also be applied to provide a uniform (in $1\leq
i\leq N$) approximation of the $Q_{iT}$ matrix using $d_{i}$-fold integrals
of Brownian motion. As another extension, we show in Appendix A.4 that the
strong approximation result also holds for series that are local to unity.
The main technical challenge in extending Theorem \ref
{thm_main_result_linear_process} to the case when $d_{i}\geq 2$ and the
local-to-unity scenario is verifying Assumption \ref{ass_sufficient}(\ref
{ass_sufficient_b}). When $d_{i}\geq 2$ and ${\Greekmath 0116} _{i}=0,$ the limiting random
matrix $Q_{i}$ involves $n$-fold integrals of Brownian motion. Based on the
literature of small ball probabilities for $n$-fold integrals of Brownian
motion \citep{chen2003quadratic,gao2003integrated}, we conjecture that \ref
{ass_sufficient}(\ref{ass_sufficient_b}) holds with $\mathfrak{c}
_{4}=2\max_{1\leq i\leq N}d_{i}$, though we were unable to locate a concrete
demonstration of this in the literature. Regarding the local-to-unity case,
existing results on small ball probabilities for Ornstein-Uhlenbeck
processes \citep{ai2016note} suggest that \ref{ass_sufficient}(\ref
{ass_sufficient_b}) holds with $\mathfrak{c}_{4}=2$. Thus, although our
current Theorem \ref{thm_main_result_linear_process} applies only to
stationary and $I(1)$ processes, we expect that it extends to series with
any order of integration and those that are local to unity, with the only
change that the exponent on the log of $N$ may differ.
Finally, we note that no additional proof is necessary to extend Theorem \ref
{thm_main_result_linear_process} to the case when $2\leq d_{i}\leq p$ if the
number of series for which $d_{i}\geq 2$ does not increase with $N$ or $T$.
\section{Results from simulations}
\label{sec_simulations}
In this section we report results from applying our method in a variety of
different settings. In these simulations we take the number of cross-section
variables to be $N=100$ and vary the number of time-series observations $T$
from 100 to 1000. For each generated sample, we calculate the $(N\times N)$
correlation matrices of: (1) the raw data; (2) the OLS regression residuals
for that data set; and (3) the true cyclical components for that data set
implied by the particular process that was used to generate that sample. For
each of these correlation matrices, let $\hat{{\Greekmath 0118}}_{j}$ denote the $j$th
largest eigenvalue of the correlation matrix. The fraction of the variance
of the sample explained by the $j$th principal component is
\begin{equation}
R_{j}^{2}=\hat{{\Greekmath 0118}}_{j}/N. \label{eqn_R2}
\end{equation}
We also calculated the number of factors that would be selected for that
sample based on the $IC_{p2}$ criterion of \cite{BaiNg2002} recommended by
\citet[p. 436]{StockWatson2016Handbook},
\begin{equation}
r^{\ast }=\arg \left\{ \underset{r\in \{0,1,...,r_{0}\}}{\min }\log \left( 1-
\frac{\mathop{\textstyle \sum }\nolimits_{j=0}^{r}\hat{{\Greekmath 0118}}_{j}}{N}\right) +r\left( \frac{N+T}{NT
}\right) \log(\min \{N,T\})\right\} , \label{eqn_bic}
\end{equation}
with $\hat{{\Greekmath 0118}}_{0}$ defined to be 0. We took $r_{0}=10$ and for each of the
different cases generated 1,000 different samples.
\subsection{Mix of unrelated stationary and nonstationary variables}
For our first example, half the variables are random walks and the other
half are white noise,
\begin{equation*}
y_{it}=\left\{
\begin{array}{cc}
y_{it-1}+{\Greekmath 0122} _{it}={\Greekmath 0122} _{it}+{\Greekmath 0122} _{i,t-1}+\cdots
+{\Greekmath 0122} _{i1} & \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{for }i=1,2,...,N/2 \\
{\Greekmath 0122} _{it} & \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{for }i=(N/2)+1,...,N
\end{array}
\right.
\end{equation*}
for $t=1,...,T.$ The innovations ${\Greekmath 0122} _{it}\sim N(0,1)$ are
independent across all $i$ and $t.$ Thus each of the $N$ variables is
completely independent of the others and there is no factor structure in the
true data-generating process.
The first two columns of Table \ref{tab_random_walk} report the results from
calculating principal components of the raw data. In a sample of $T=100$
observations, the first three principal components seem to account for 38\%
of the variance of the full set of $N=100$ variables. This result is
entirely spurious, and illustrates the cautions raised by \cite
{onatski2021spurious} about using PCA when some of the variables are
nonstationary. The criterion (\ref{eqn_bic}) would always lead us
incorrectly to conclude that there is more than one factor in a sample of
size $T=100.$ This problem in the apparent number of factors gets even worse
when the sample size increases. The latter is the expected result, since (
\ref{eqn_bic}) is in the class of criteria for which
\citet[p.
602]{onatski2021spurious} demonstrated that the number of factors selected
diverges as $N$ and $T$ go to infinity.
\begin{table}[tbph]
\caption{Mixture of independent random walks and white noise}
\label{tab_random_walk}
\begin{tabular}{c|cc|cc|cc|cc|cc|cc}
\hline
& \multicolumn{2}{c|}{Raw data} & \multicolumn{2}{c|}{$c_t$ ($h$=24)} &
\multicolumn{2}{c|}{$\hat{c}_t$ ($h$=24)} & \multicolumn{2}{c|}{$\hat{c}_t$ (
$h$=12)} & \multicolumn{2}{c|}{$\hat{c}_t$ ($h$=8)} & \multicolumn{2}{c|}{$
\hat{c}_t$ ($h$=1)} \\
\cmidrule{2-3} \cmidrule{4-5} \cmidrule{6-7} \cmidrule{8-9} \cmidrule{10-11}
\cmidrule{12-13} $j$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ & $R^2$
& $r^*$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ \\
& (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12)
\\ \hline
$T=100$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 0 & --- & 1 & --- & 13 & --- & 96 & --- & 100 \\
1 & 22.6 & 0 & 16.0 & 0 & 11.9 & 26 & 10.3 & 48 & 8.4 & 4 & 4.1 & 0 \\
2 & 9.9 & 67 & 12.0 & 6 & 9.4 & 54 & 8.5 & 35 & 7.1 & 0 & 3.8 & 0 \\
3 & 5.6 & 33 & 8.3 & 94 & 7.3 & 19 & 6.9 & 4 & 6.1 & 0 & 3.6 & 0 \\ \hline
$T=200$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 0 & --- & 0 & --- & 78 & --- & 100 & --- & 100 \\
1 & 22.7 & 0 & 9.4 & 0 & 10.2 & 0 & 6.6 & 22 & 5.2 & 0 & 2.9 & 0 \\
2 & 9.4 & 13 & 7.8 & 3 & 8.2 & 8 & 5.7 & 0 & 4.6 & 0 & 2.7 & 0 \\
3 & 5.3 & 87 & 6.6 & 97 & 6.5 & 92 & 4.9 & 0 & 4.2 & 0 & 2.6 & 0 \\ \hline
$T=400$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 23 & --- & 3 & --- & 100 & --- & 100 & --- & 100 \\
1 & 22.5 & 0 & 5.9 & 51 & 6.5 & 34 & 4.3 & 0 & 3.5 & 0 & 2.2 & 0 \\
2 & 9.4 & 0 & 5.2 & 22 & 5.5 & 47 & 3.8 & 0 & 3.2 & 0 & 2.1 & 0 \\
3 & 5.1 & 100 & 4.5 & 4 & 4.7 & 16 & 3.5 & 0 & 3.0 & 0 & 2.1 & 0 \\ \hline
$T=600$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 99 & --- & 79 & --- & 100 & --- & 100 & --- & 100 \\
1 & 22.4 & 0 & 4.6 & 1 & 5.0 & 21 & 3.4 & 0 & 2.9 & 0 & 1.9 & 0 \\
2 & 9.3 & 0 & 4.1 & 0 & 4.4 & 0 & 3.1 & 0 & 2.7 & 0 & 1.9 & 0 \\
3 & 5.2 & 100 & 3.7 & 0 & 3.9 & 0 & 2.9 & 0 & 2.5 & 0 & 1.8 & 0 \\ \hline
$T=800$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 100 & --- & 100 & --- & 100 & --- & 100 & --- & 100 \\
1 & 22.7 & 0 & 4.0 & 0 & 4.2 & 0 & 3.0 & 0 & 2.6 & 0 & 1.8 & 0 \\
2 & 9.3 & 0 & 3.6 & 0 & 3.7 & 0 & 2.7 & 0 & 2.4 & 0 & 1.7 & 0 \\
3 & 5.0 & 100 & 3.3 & 0 & 3.4 & 0 & 2.6 & 0 & 2.3 & 0 & 1.7 & 0 \\ \hline
$T=1000$ & & & & & & & & & & & & \\
0 & --- & 0 & --- & 100 & --- & 100 & --- & 100 & --- & 100 & --- & 100 \\
1 & 22.6 & 0 & 3.5 & 0 & 3.7 & 0 & 2.7 & 0 & 2.4 & 0 & 1.7 & 0 \\
2 & 9.3 & 1 & 3.2 & 0 & 3.3 & 0 & 2.5 & 0 & 2.2 & 0 & 1.7 & 0 \\
3 & 5.1 & 99 & 3.0 & 0 & 3.0 & 0 & 2.4 & 0 & 2.1 & 0 & 1.6 & 0 \\
\hline\hline
\end{tabular}
\vspace{5pt} \newline
Notes to Table \ref{tab_random_walk}. $R^2$ indicates the percentage of
total variance accounted for by the $j$th principal component for $j=1,2$ or
3. $r^*$ indicates the percentage of samples for which the criterion (\ref
{eqn_bic}) selects the number of factors to be $j=0,1,2,$ or $\ge 3$. In
every case, the true number of factors is $r=0$ and the cross-section
dimension is $N=100$.
\end{table}
The next two columns of Table \ref{tab_random_walk} report what the results
would be if we somehow knew the true cyclical component of each variable.
For this case, we applied PCA to a sample of $T-h$ observations for which
the $i$th observed variable for $t=h+1,h+2,...,T$ is given by
\begin{equation}
c_{it}=\left\{
\begin{array}{cc}
y_{it}-y_{i,t-h}={\Greekmath 0122} _{it}+{\Greekmath 0122} _{i,t-1}+\cdots +{\Greekmath 0122}
_{i,t-h+1} & \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{for }i=1,2,...,N/2 \\
{\Greekmath 0122} _{it} & \relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{for }i=(N/2)+1,...,N
\end{array}
\right. . \label{eqn_true_cycle_random_walk}
\end{equation}
Note that the cyclical components in (\ref{eqn_true_cycle_random_walk}) can
be serially correlated, but this autocorrelation vanishes for observations
separated by more than $h$ periods. The cyclical components $c_{it}$ thus
satisfy by construction the conditions under which \cite{BaiNg2002}
demonstrated that (\ref{eqn_bic}) would asymptotically select the correct
number of factors. We find in our simulations that (\ref{eqn_bic}) does
indeed correctly conclude there is no factor structure for these data sets
provided the number of time-series observations is 600 or larger. For
smaller $T$ it is less reliable. The reason is that there is a small-sample
version of the \cite{onatski2021spurious} spurious factor problem that
arises from the serial correlation of some of the variables that is induced
by the definition of the true cyclical component in equation (\ref
{eqn_true_cycle_random_walk}). If $T$ is large enough, this problem goes
away, but for smaller $T$ it can make a difference.
Columns (5) and (6) examine the case where the analysis is based on the
residuals from running an OLS regression on the raw data $Y_{it}$ for all
variables $i=1,...,N$ without making any judgments about which variables are
stationary and which are not. For large samples, the results are similar to
those that we would obtain if we somehow knew the exact correct
transformation to use for every variable.
The small-sample problem in columns (3)-(6) results from the serial
correlation that is a consequence defining the cyclical component to be the
error from a 24-period-ahead forecast. Columns (7) and (8) report results if
we instead were to look for common factors in the 12-period-ahead forecast
errors. This typically would reach the correct conclusion even in a sample
of only $T=200$ observations. Columns (9) and (10) consider 8-period-ahead
forecast errors, such as our suggested cyclical calculation would use for
quarterly data. The results indicate that if we have more than 50 years of
data ($T=600$ for monthly data or $T=200$ for quarterly data), conducting
PCA on the two-year-ahead OLS forecast residuals ($h=24$ for monthly data or
$h=8$ for quarterly data) is reasonably reliable. With less than 50 years of
data, some researchers might want to use a smaller value for $h$ or place
less reliance on (\ref{eqn_bic}) as a criterion for selecting the number of
factors.
The last two columns of Table \ref{tab_random_walk} examine looking for
common factors in the one-period-ahead forecast errors. In these
simulations, this reaches the correct conclusion 100\% of the time that
there are zero factors in these data sets even for a sample of $T=100$
observations. Thus our proposed method appears to be quite reliable if the
interest is in identifying common factors behind one-period-ahead forecast
errors. However, one-period-ahead forecast errors are more sensitive to
outliers. This is an important consideration, as will be demonstrated in our
analysis of actual data in Section \ref{sec_fred_md}.
\subsection{Mix of unrelated stationary variables with differing persistence}
In our second example, for $i=1,2,...,N/2$ the variables are generated by a
stationary but persistent AR(1) process:
\begin{equation*}
y_{i1}\sim N(0,1/(1-{\Greekmath 011A} ^{2}))
\end{equation*}
\begin{equation}
y_{it}={\Greekmath 011A} y_{i,t-1}+{\Greekmath 0122} _{it}\relax\protect\ifmmode\expandafter\text@\else\expandafter\mbox\fi{ \ \ for }t=2,3,...,T.
\label{eq_simulation_AR1}
\end{equation}
The remaining $N/2$ variables are white noise $(y_{it}={\Greekmath 0122} _{it}$
for $i=(N/2)+1,...,N).$ Our example uses ${\Greekmath 011A} =0.99,$ so all the variables
are stationary but half of them are highly persistent. The innovations $
{\Greekmath 0122} _{it}\sim N(0,1)$ are independent across all $i$ and $t,$ so
there is no factor structure in the true data-generating processes.
Columns (1) and (2) of Table \ref{tab_stationary_coint}\ report the results
from applying PCA to the raw data. Note that for this example, the raw data
themselves satisfy the \cite{BaiNg2002} conditions for asymptotic validity
of PCA. Nevertheless, even in a sample of size $T=1000,$ the first principal
component alone appears to explain a third of the data, and the criterion in
(\ref{eqn_bic}) would always conclude incorrectly that there is at least one
factor. This is a small-sample manifestation of the \cite
{onatski2021spurious} spurious factor phenomenon. In a sufficiently large
sample, this problem would go away. But $T=1000$ is not large enough for
persistence characterized by ${\Greekmath 011A} =0.99.$
\begin{table}[tbph]
\caption{Monte Carlo results in other settings}
\label{tab_stationary_coint}
\begin{tabular}{c|cc|cc|cc|cc|cc|cc|cc|}
\hline
& \multicolumn{6}{c|}{${\Greekmath 011A}=0.99$} & \multicolumn{4}{c|}{Cointegrated} &
\multicolumn{4}{c|}{Mixed} \\
\cmidrule{2-7} \cmidrule{8-11} \cmidrule{12-15} & \multicolumn{2}{c|}{Raw
data} & \multicolumn{2}{c|}{$\hat{c}_t$ ($h$=24)} & \multicolumn{2}{c|}{$
\hat{c}_t$ ($h$=1)} & \multicolumn{2}{c|}{Raw data} & \multicolumn{2}{c|}{$
\hat{c}_t$ ($h $=24)} & \multicolumn{2}{c|}{Raw data} & \multicolumn{2}{c|}{$
\hat{c}_t$ ($h $=24)} \\
\cmidrule{2-3} \cmidrule{4-5} \cmidrule{6-7} \cmidrule{8-9} \cmidrule{10-11}
\cmidrule{12-13} \cmidrule{14-15} $j$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ & $R^2$
& $r^*$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ & $R^2$ & $r^*$ \\
& (1) & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) & (10) & (11) & (12) &
(13) & (14) \\ \hline
$T=100$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 0 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0 \\
1 & 43.0 & 98 & 15.8 & 3 & 4.1 & 0 & 46.1 & 100 & 40.3 & 100 & 27.5 & 0 &
23.8 & 0 \\
2 & 3.7 & 2 & 10.6 & 97 & 3.8 & 0 & 2.8 & 0 & 3.4 & 0 & 18.1 & 0 & 11.2 & 1
\\
3 & 2.8 & 0 & 5.5 & 0 & 3.6 & 0 & 2.5 & 0 & 3.1 & 0 & 8.5 & 100 & 8.7 & 99
\\ \hline
$T=200$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 0 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0 \\
1 & 42.5 & 83 & 10.8 & 0 & 2.9 & 0 & 48.1 & 100 & 45.7 & 100 & 25.1 & 0 &
22.6 & 0 \\
2 & 3.3 & 17 & 8.7 & 1 & 2.7 & 0 & 2.1 & 0 & 2.3 & 0 & 18.0 & 0 & 9.8 & 0 \\
3 & 2.3 & 0 & 7.0 & 99 & 2.6 & 0 & 2.0 & 0 & 2.2 & 0 & 8.7 & 100 & 7.9 & 100
\\ \hline
$T=400$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 5 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0 \\
1 & 40.5 & 52 & 6.3 & 33 & 2.2 & 0 & 48.8 & 100 & 46.5 & 100 & 23.8 & 0 &
21.1 & 0 \\
2 & 3.4 & 47 & 5.4 & 45 & 2.1 & 0 & 1.8 & 0 & 1.8 & 0 & 17.9 & 0 & 6.4 & 0
\\
3 & 2.2 & 1 & 4.8 & 18 & 2.1 & 0 & 1.7 & 0 & 1.7 & 0 & 8.8 & 100 & 5.4 & 100
\\ \hline
$T=600$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 91 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0
\\
1 & 38.6 & 38 & 4.9 & 9 & 1.9 & 0 & 49.3 & 100 & 47.0 & 100 & 23.3 & 0 & 20.8
& 0 \\
2 & 3.5 & 58 & 4.3 & 0 & 1.9 & 0 & 1.6 & 0 & 1.6 & 0 & 17.9 & 0 & 4.9 & 4 \\
3 & 2.3 & 5 & 3.9 & 0 & 1.8 & 0 & 1.5 & 0 & 1.6 & 0 & 8.9 & 100 & 4.3 & 96
\\ \hline
$T=800$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 100 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0
\\
1 & 37.0 & 33 & 4.1 & 0 & 1.8 & 0 & 49.4 & 100 & 47.0 & 100 & 23.2 & 0 & 20.5
& 33 \\
2 & 3.5 & 60 & 3.7 & 0 & 1.7 & 0 & 1.5 & 0 & 1.5 & 0 & 17.7 & 0 & 4.1 & 43
\\
3 & 2.4 & 8 & 3.3 & 0 & 1.7 & 0 & 1.5 & 0 & 1.5 & 0 & 9.2 & 100 & 3.7 & 24
\\ \hline
$T=1000$ & & & & & & & & & & & & & & \\
0 & --- & 0 & --- & 100 & --- & 100 & --- & 0 & --- & 0 & --- & 0 & --- & 0
\\
1 & 35.5 & 33 & 3.6 & 0 & 1.7 & 0 & 49.5 & 100 & 47.1 & 100 & 22.8 & 0 & 20.4
& 86 \\
2 & 3.5 & 59 & 3.3 & 0 & 1.7 & 0 & 1.5 & 0 & 1.5 & 0 & 17.6 & 0 & 3.7 & 14
\\
3 & 2.5 & 8 & 3.0 & 0 & 1.6 & 0 & 1.4 & 0 & 1.4 & 0 & 9.2 & 100 & 3.3 & 0 \\
\hline\hline
\end{tabular}
\vspace{5pt} \newline
Notes to Table \ref{tab_stationary_coint}. $R^2$ indicates the percentage of
total variance accounted for by the $j$th principal component for $j=1,2$ or
3. $r^*$ indicates the percentage of samples for which the criterion (\ref
{eqn_bic}) selects the number of factors to be $j=0,1,2,$ or $\ge 3$. In
columns (1)-(6), the true number of factors is $r=0$. In columns (7)-(14),
the true number of factors is $r=1$. In every case, the cross-section
dimension is $N=100$.
\end{table}
In columns (3)-(4) we apply PCA to the residuals from a 24-period-ahead
forecasting regression. Again we estimated the same regression for all
variables, whether persistent or not. And again using regression residuals
solves the problem pretty reliably in samples larger than $T=600.$ If we
look for a factor structure in the one-period-ahead regression residuals as
in columns (5)-(6), the problem is solved 100\% of the time even in a sample
of $T=100.$
\subsection{Results for other data-generating processes}
We also report simulations for a number of other data-generating processes.
For details of these and other simulations see online Appendix C. We first
considered a dataset in which half the variables are cointegrated with a
single common factor and the other are independent white noise. Columns (7)
and (8) of Table \ref{tab_stationary_coint}\ reproduce the well-known result
that PCA on raw nonstationary cointegrated data can correctly pick out the
single common factor; for more discussion see \cite{harris1997} and \cite
{Onatski_cointegration_2018}. We find that this result is preserved when all
$N$ series are prefiltered as proposed here, as seen in columns (9)-(10) of
Table \ref{tab_stationary_coint}.
We also considered an example where the common factor is stationary while
the idiosyncratic components are a mix of stationary and nonstationary
processes. When PCA is applied directly to the observed data $y_{it},$ the
familiar \cite{onatski2021spurious} result is observed in column (11) of
Table \ref{tab_stationary_coint}: higher-order factors spuriously appear to
explain a large amount of the variance of the data. When PCA is applied
instead to the $h=24$-period-ahead regression forecast residuals $\hat{c}
_{it},$ the contribution of higher-order factors $\hat{F}_{jt}$ for $j\geq 1$
is substantially lower, though the selection criterion (\ref{eqn_bic}) would
still typically incorrectly conclude that $r>2$ for $T\leq 600.$ The correct
conclusion $(r=1)$ would be reached 86\% of the time when $T=1000.$
Moreover, we found that the average correlation between the true realization
of $F_{t}$ for a particular simulation and the first estimated factor $\hat{F
}_{1t}$ of the OLS regression residuals from that simulation is $0.98$ for $
T\geq 600$. In other words, the first principal component of the OLS
regression residuals accurately uncovers the true single common feature of
these data.
Table C1 in online Appendix C examines local-to-unity and fractionally
integrated processes. The results in these simulations turn out to be in
between the unit-root example in Table \ref{tab_random_walk} and the
stationary-but-persistent example in columns (1)-(6) of Table \ref
{tab_stationary_coint}. Table C2 compares our approach with the PANIC method
of \cite{BaiNg2004}. They proposed to take first-differences of all the
original data, apply PCA to the changes, and then accumulate the resulting
principal components. We find PANIC performs very similarly to our approach
for these examples when we set $h=1.$
Overall, these results confirm the asymptotic theory that applying PCA to
OLS regression residuals is a promising approach to handling both
nonstationarity and stationary persistence of unknown form provided that the
time-series dimension $T$ is reasonably large.
\section{Characterizing a large macroeconomic data set}
\label{sec_fred_md}
The use of large macroeconomic data sets was pioneered by \cite
{stock1999forecasting}. Their goal was to use the information of 168
different macroeconomic variables to produce better forecasts of inflation.
They found that the first principal component of indicators of real economic
activity produced the best inflation forecasts over the period 1959:1 to
1997:9. Their findings led to the development of the Chicago Fed National
Activity Index (CFNAI).
\cite{mccracken2016fred} developed the FRED-MD database whose 2015:4 vintage
covered 134 macroeconomic variables. \ These include monthly measures in
eight broad categories: (1) output and income; (2) labor market; (3)
housing; (4) consumption, orders, and inventories; (5) money and credit; (6)
interest and exchange rates; (7) prices; and (8) stock market. This data set
offers benefits of continuity and continuous updating and is the basis for
the analysis in this paper.
In previous applications of PCA to large macroeconomic data sets, each of
the variables needed to be transformed using a detrending method that was
selected individually for each series. For details of how this has been done
for the CFNAI see \cite{cfnai_background} and for FRED-MD see the data
appendix to \cite{mccracken2016fred}. Figure \ref{fig_fred_selected_series}
illustrates these transformations for three important macroeconomic
indicators. The first column plots the raw data, while the second column
plots the data as transformed by \cite{mccracken2016fred} using the same
data set as in their original paper. Everyone agrees that industrial
production (row 1) is nonstationary, and all previous researchers have used
first differences of the log of industrial production shown in panel (1,2).
While there is little doubt that this is a valid way to generate a
stationary series for this variable, monthly growth rates of industrial
production exhibit a lot of high-frequency fluctuations around the dominant
cyclical patterns. For the unemployment rate (row 2), it is less clear
whether the series should be regarded as stationary. \cite{mccracken2016fred}
used first differences of unemployment, which behave quite differently from
the level. The purchasing managers composite index from the Institute of
Supply Management (row 3) appears to be stationary, and \cite
{mccracken2016fred} entered this series directly into PCA without any
transformation.
\begin{figure}[tbp]
\caption{Level, transformed value, and cyclical component of industrial
production, unemployment, and PMI Composite, 1962:3 to 2014:12}
\label{fig_fred_selected_series}\centering
\includegraphics[width=\textwidth, height = 9.5cm]{selected_graphs.png}
\end{figure}
The third column of Figure \ref{fig_fred_selected_series} plots the cyclical
components of industrial production, unemployment, and PMI as estimated by
the residuals of the OLS regression (\ref{eqn_recommended_regression}) with $
h=24$ and $p=12$.\footnote{
For those series that \cite{mccracken2016fred} transformed using logs, first
differences of logs, or second differences of logs (their transformations
4-6), we simply took the log of the variable before performing the
regression. Thus for example the series plotted in the upper left panel of
Figure \ref{fig_fred_selected_series} is 100 times the natural logarithm of
the industrial production index. For those series that they used as is, as
first differences, or second differences (their transformations 1-3), we
simply used the variable as is. They employed a special transformation (7)
for nonborrowed reserves. One would have expected to take logs of a variable
like this, but the variable took on negative values in 2008. For this series
their transformation was $y_{it}=\Delta (x_{t}/x_{t-1}-1.0)$ and we used $
y_{it}=$ $x_{t}/x_{t-1}.$} PMI is almost impossible to predict two years in
advance, and our cyclical component is almost identical to the original
series. Thus both our method and the traditional approach use this variable
essentially as is. There is some but not much predictability of the
unemployment rate at the two-year horizon, so for this variable our
transformation much more closely resembles the original series than it does
the first-difference transformation. For industrial production, our approach
takes out the broad trend while retaining the essential cyclical behavior
observed in the raw data. The three variables in the third column, unlike
those in the second column, all share a common characterization of what is
happening over the business cycle. Consistent with a long tradition in
business cycle research, when plotted this way PMI appears as a leading
indicator, industrial production as a coincident indicator, and unemployment
as a coincident or lagging indicator, with all three clearly following the
same cycle.
The top panel of Figure \ref{fig_macro_fred_2015} plots the first principal
component of the transformed series arrived at by \cite{mccracken2016fred}.
\footnote{
We generated this figure using the exact data and code posted at
https://research.stlouisfed.org/econ/mccracken/fred-databases/. Note that we
have multiplied the series by $-1$ in order to give it the property that the
factor declines in recessions, and that their Figure 3 plots the
accumulations $(s_{1t}=\mathop{\textstyle \sum }\nolimits_{j=1}^{t}\tilde{f}_{1t})$ whereas our
graph shows $\tilde{f}_{1t}$ itself.} This inherits some of the
high-frequency fluctuations seen in the (1,2) and (2,2) panels of Figure \ref
{fig_fred_selected_series}. Indeed, \cite{mccracken2016fred} regarded this
series as too volatile to reliably identify business cycles and turning
points, and instead plotted in their Figure 3 the accumulation of this
series. The CFNAI (shown in panel 2 of Figure \ref{fig_macro_fred_2015}) is
very similar to the first principal component of the FRED-MD macro data set.
\begin{figure}[tbp]
\caption{First PC of FRED-MD variables as transformed by \protect\cite
{mccracken2016fred}, the Chicago Fed National Activity Index, and first PC
of cyclical components of FRED-MD variables, 1962:3 to 2014:12}
\label{fig_macro_fred_2015}\centering
\includegraphics[width=\textwidth, height = 9.5cm]{factor_cfnai_2015.png}
\end{figure}
The first principal component of the estimated cyclical components of the
variables in the data set is plotted in the bottom panel of Figure \ref
{fig_macro_fred_2015}. Unlike the CFNAI, this provides a very clean summary
of historical business cycles. There is another interesting difference
between the third panel and the first two. The NBER defines the
business-cycle trough (the end point of the shaded regions) as the low point
in the \textit{level} of overall economic activity. For example, in the
first month of a new expansion, the unemployment rate is still very high,
but it has started to come down. Our series in the bottom panel of Figure
\ref{fig_macro_fred_2015} captures this feature very well, reaching a trough
at exactly the point identified by the NBER Business Cycle Dating Committee.
By contrast, the low point in the series plotted in the first two panels
typically comes more towards the middle of the recession. This is because
the \textit{rate of decline }of real output (a common raw input in the
variables as usually transformed) starts to ease well before the recession
has ended. Both in terms of the cleanness of the series and its timing, we
would suggest that our approach offers a better characterization of the
state of the U.S. business cycle over this sample period.
\subsection{Choice of $h$}
Our method identifies common components in $h$-period-ahead forecast errors.
Any value of $h\geq 1$ is sufficient to produce a stationary series. Which
value of $h$ to use depends on which common components the researcher is
interested in.
McCracken and Ng applied the \cite{BryBoschan} algorithm for dating
business-cycle turning points to the first principal component that resulted
from their analysis. They found that the resulting series was in agreement
with NBER recession dates only 65\% of the time, meaning that 65\% of NBER
recession periods fell within the peak-to-trough phases of the series. We
were able to reproduce this result by adapting the code in \cite
{StockWatsonTurningPoints}.\footnote{
As discussed in \citet[f.n. 2]{StockWatsonTurningPoints}, the first step in
applying the Bry and Boschan algorithm is to take a centered moving average
of the series for which turning points are to be assigned, which in this
case is a moving average of the first principal component. Specifying nma =
12 seems to come closest to reproducing the numbers reported in
\citet[p.
583]{mccracken2016fred}. Our numbers, however, were not in all instances
identical. This appears to be in part because McCracken and Ng were
reporting results for a data set that has some slightly different series
from those in 2015-04.csv. In the case of the percent agreement with NBER
recession dates, our calculation (65\%) exactly reproduces theirs. Other
numbers reported here in the text are based on our reproduction using the
2015-04.csv data, which sometimes differ from those reported in
\citet[p.
583]{mccracken2016fred}.} We found a similar result (67\% agreement) when
Bry and Boschan is applied to the CFNAI.
The top panel of Figure \ref{fig_Bry_Boschan} plots the percent agreement
with NBER recessions that results from our method for every value of $h$
between 1 and 36 months. The common component of $h=1$- or 2-month-ahead
forecast errors has less correspondence with NBER recession dates than
either the McCracken and Ng series or the CFNAI. We attribute this to the
fact that replacing stationary series like the PMI in the bottom row of
Figure \ref{fig_fred_selected_series} with a 1-month-ahead forecast error
removes what we would normally think of as the business-cycle indicator
provided by the series. By contrast, for every $h\geq 3,$ our approach has a
closer correspondence to NBER recessions than either McCracken-Ng or the
CFNAI. The maximal agreement (99\%) is obtained using $h=25$. Choosing $h=24$
is very similar (98\%).
\begin{figure}[tb]
\caption{Agreement with NBER-dated recessions and number of series
exhibiting outliers as a function of forecast horizon $h$, 1967:2-2014:12}
\label{fig_Bry_Boschan}
\begin{center}
\includegraphics[width=\textwidth, height = 9.5cm]{Fig_Bry_Boschan-eps-converted-to.pdf}
\end{center}
\par
{\small {\ \emph{Notes to Figure \ref{fig_Bry_Boschan}.} Horizontal axis:
forecast horizon $h$. Top panel: percent of NBER-dated recession periods
that fall within Bry-Boschan-designated downturns. Bottom panel: number of
series exhibiting outliers. Horizontal lines denote results for the first
principal component found by \cite{mccracken2016fred}.} }
\end{figure}
McCracken and Ng noted that agreement of the original $\hat{f}_{t}^{MN}$
with NBER expansion dates was even weaker. We found that only 51\% of the
expansion dates based on $\hat{f}_{t}^{MN}$ were also designated as
expansion by NBER and still only 61\% agreement using $\hat{F}_{t}^{MN}$.
Our series with $h=24$ does a little better than either of these with 65\%
agreement.
\subsection{Outliers}
Previous users of large macro data sets devoted a lot of attention to
outliers and implemented procedures to mitigate their influence. Prior to
the COVID-19 recession of 2020, the CFNAI discarded observations that were
more than six times the interquartile range, as did \cite
{stock1999forecasting} in some of their analysis. \cite{mccracken2016fred}
discarded observations that were more than ten times the interquartile
range. This criterion identifies 79 different observations on 22 different
variables as outliers in the 1960:3 to 2014:12 data set.
We noted in footnote \ref{fn_CLT} that as a result of the Central Limit
Theorem, we should expect fewer outliers when we calculate principal
components of $h$-period-ahead forecasts for larger values of $h$. To
calculate outliers for our method, we constructed forecast errors from
leave-one-out regressions and designated values exceeding ten times the
interquartile range as outliers.\footnote{
That is, we calculated $\tilde{c}_{it}=y_{it}-\tilde{{\Greekmath 010B}}_{i,t}\tilde{z}
_{i,t-h}$ with $\tilde{{\Greekmath 010B}}_{i,t}=\left( \mathop{\textstyle \sum }\nolimits_{s=1,s\neq t}^{T}
\tilde{z}_{i,s-h}\tilde{z}_{i,s-h}^{\prime }\right) ^{-1}\left(
\mathop{\textstyle \sum }\nolimits_{s=1,s\neq t}^{T}\tilde{z}_{i,s-h}y_{is}\right) $ estimated
separately for each $i$ and $t$ and then divided $\tilde{c}_{i.t}$ by its
observed interquartile range.} The number of series exhibiting outliers is
plotted as a function of $h$ in the bottom panel of Figure \ref
{fig_Bry_Boschan}. Using $h=1$ results in 31 series exhibiting outliers,
substantially more than in calculating the \cite{mccracken2016fred}
principal components. Many of the additional outliers come from interest
rate spreads. \cite{mccracken2016fred} used these as is without filtering,
which turns out to be close to our definition of the cyclical component of
these variables. By contrast, using $h=1$-month-ahead forecast errors
highlights some of the unusual behavior of interest rates during the Volcker
monetary contraction in 1980-81. The number of outliers steadily decreases
as $h$ is increased. Only two series exhibit outliers for $h\geq 21.$
\footnote{
These are total bank reserves and nonborrowed reserves in the Fed's response
to the financial crisis.}
Our recommended procedure is to use two-year-ahead regression residuals and
make no corrections for outliers. The series that we have plotted in the
bottom panel of Figure \ref{fig_macro_fred_2015} is the unadjusted first
principal component of the full set of OLS residuals $\hat{c}_{it}.$
Outliers are an even bigger issue when data for 2020 are included. For the
2024:12 vintage of FRED-MD, the McCracken-Ng procedure would identify 40 of
the 126 variables as all being outliers in the single month of 2020:4.
Despite dropping all of these 40 observations, the first principal component
calculated using their algorithm shows an enormous decline in this month.
Indeed, in order to include the 2020 observations in the top panel of Figure
\ref{fig_macro_fred_2022}, the scale must be so large that it makes all the
previous cyclical fluctuations barely noticeable. The CFNAI modified its
procedure for dealing with anomalous observations to handle these
observations. Even so the CFNAI still displays an unprecedented drop in
2020, as seen in the second panel.
\begin{figure}[tbh]
\caption{First PC of FRED-MD variables as transformed by \protect\cite
{mccracken2016fred}, the Chicago Fed National Activity Index, and first and
second PC of cyclical components of FRED-MD variables, 1962:3 to 2024:9}
\label{fig_macro_fred_2022}\centering
\includegraphics[width=\textwidth, height = 11cm]{PC1_and_PC2-eps-converted-to.pdf}
\end{figure}
By contrast, only two variables are identified as outliers for 2020:4 for
purposes of our approach, these being new claims for unemployment insurance
and the number unemployed for less than 5 weeks. The result of applying our
procedure to the FRED-MD database available as of December of 2024 with no
corrections for outliers is displayed in the third panel of Figure \ref
{fig_macro_fred_2022}. Note that, unlike the top two panels, our series
describes the downturn in 2020 on a comparable scale as earlier recessions,
although our series indicates that the speed of the downturn was
unprecedented, as was the growth in the first two months of the recovery.
Also in contrast to the first two panels, our series indicates (correctly,
in our view) that the economy did not fully recover from the COVID-19 shock
until September of 2021. The sharp spike up in our series in April of 2022
reflects the fact that most macro variables were substantially higher in
April 2022 than one would have predicted based on observations in April
2020. This conclusion is also consistent with the aggressive actions of
policy makers in the spring of 2020. Our series further indicates that
economic activity remained unusually strong through the fall of 2023.
Another issue comes from discontinued, newly added, or missing variables.
\cite{mccracken2016fred} adapted the \cite{stock2002forecasting} algorithm
for unbalanced panels, though they found in their original data set that the
results are essentially identical if one simply drops variables as needed to
create a balanced panel. For our application, we have simply calculated
principal components of $\hat{c}_{it}$ on a balanced panel, though there is
no obstacle to applying the \cite{stock2002forecasting} algorithm to an
unbalanced panel of $\hat{c}_{it}.$
\subsection{Uses of macroeconomic cyclical factors}
A key use of PCA is to summarize the statistical information in a large
cross section of indicators. The movement in variable $i$ that is captured
by the $j$th factor alone is given by $\hat{{\Greekmath 0115}}_{ij}\hat{f}_{jt}.$
Since $\hat{c}_{it}$ is normalized to have unit variance, the fraction of
the variance of the stationary component of variable $i$ that is explained
by the $j$th macro factor is given by $\hat{{\Greekmath 0115}}_{ij}^{2}T^{-1}\mathop{\textstyle \sum }
\nolimits_{t=1}^{T}\hat{f}_{jt}^{2}.$
Table B1 in the online appendix reports the $R^{2}$ explained by the first
macroeconomic cyclical factor for each of the 119 variables used in our
analysis of the 2023 vintage data set. The first factor alone accounts for
almost 2/3 of the variance of typical indicators of real output or income
and more than half of the variance of typical indicators of labor-market
conditions. The first cyclical factor is far less successful at describing
financial indicators and nominal prices. It is interesting that when we add
the second cyclical factor, the $R^{2}$ for the median price indicator rises
to 66\%. The first factor thus seems mainly to capture real economic
conditions and the second characterizes nominal prices and interest rates.
The fourth panel in Figure \ref{fig_macro_fred_2022} plots the second
cyclical factor. This by construction is orthogonal to the first, and often
continues to fall even as the recovery in real economic activity is
beginning. This is consistent with the view that nominal variables may
respond sluggishly to business-cycle developments. It describes events
beginning in 2022 as a third big U.S. inflation wave, though less dramatic
than the big inflations of 1973-74 and 1979-81.
We next explore the use of the cyclical factors in forecasting. \cite
{stock1999forecasting} demonstrated that the first principal component of a
large data set of real macroeconomic variables could be very helpful for
forecasting inflation. Their finding gave rise to the Chicago Fed National
Activity Index (CFNAI), a PCA-based indicator that is still widely used
today. We compare the usefulness for forecasting of the CFNAI (denoted $\hat{
f}_{t}^{CF}),$ the first or second principal component calculated using the
algorithm and data set of \cite{mccracken2016fred} (denoted $\hat{f}
_{t}^{MN1}$ and $\hat{f}_{t}^{MN2},$ respectively$),$ or the first or second
principal component of the forecasting residuals (denoted $\hat{f}
_{t}^{HMX1} $ or $\hat{f}_{t}^{HMX2})$ calculated from the FRED-MD data set
using equation (\ref{eqn_HX_f1_definition}). Our approach to comparing
different forecasts is similar to that used by \cite{stock1999forecasting}
and \cite{mccracken2016fred}.
A particular model $m$ uses a set of variables $\mathbf{x}_{t}^{m}$ that are
observed at date $t$ to try to forecast the value of a variable of interest $
y_{t+h}^{h}$ that will not be observed until $t+h$:
\begin{equation}
y_{t+h}^{h}=\boldsymbol{{\Greekmath 0119} }^{m\prime }\mathbf{x}_{t}^{m}+u_{t+h}^{m,h}.
\label{eqn_forecast_regression}
\end{equation}
We estimated the value of $\boldsymbol{{\Greekmath 0119} }^{m}$ by OLS regression on the
subsample $t=T_{0},T_{0}+1,...,T_{1}-h-1$ and used these coefficients to
forecast $y_{T_{1}}^{h}.$ We then augmented the sample by one observation,
estimating the regression over $t=T_{0},T_{0}+1,...,T_{1}-h$ and using those
coefficients to forecast $y_{T_{1}+1}^{h}.$ We repeated this for an
evaluation period $T_{1}$ to $T_{2}$ and calculated the average squared
forecast error over this evaluation period. We performed this analysis using
three different evaluation periods. The first evaluation period is specified
by $T_{1}=1970$:1 to $T_{2}=$ 1996:12, which was the evaluation period in
the original study by \cite{stock1999forecasting}. The second evaluation
period is $T_{1}=$ 1997:1 to $T_{2}=$ 2014:12, which corresponds to the new
data used by \cite{mccracken2016fred} that were not available to \cite
{stock1999forecasting}. The third evaluation period is $T_{1}=$ 2015:1 to $
T_{2}=$ 2024:9, which is the new data available since publication of \cite
{mccracken2016fred}. The models we considered were a pure autoregressive
model, $\mathbf{x}_{t}^{AR}=(1,y_{t}^{1},y_{t-1}^{1},...,y_{t-5}^{1})^{
\prime },$ and models that add to the autoregressive model six lags of one
of the principal components estimates. For example, $\mathbf{x}_{t}^{CF}=(
\mathbf{x}_{t}^{AR\prime },\hat{f}_{t}^{CF},\hat{f}_{t-1}^{CF},...,\hat{f}
_{t-5}^{CF})^{\prime }.$ This differs a little from the forecast evaluations
performed by \cite{stock1999forecasting} and \cite{mccracken2016fred} in
that these authors used $BIC$ to select different lag lengths for the
autoregressive and principal components and for each subsample, whereas we
set the lag length to six for every evaluation. Also, since $\hat{f}
_{t}^{CF} $ is only available beginning in 1967:3, we used $T_{0}=$ 1967:9
as the first date for estimation of all models. In every case, for each $T$
we re-estimated the coefficients $\boldsymbol{\hat{{\Greekmath 0119}}}^{m}$ for the
forecasting regression (\ref{eqn_forecast_regression}) using an expanding
data set ending $h$ periods before the variable being forecast.
A separate question is the data set used to estimate the factors $\hat{f}
_{t} $ themselves. Insofar as the factors are only identified up to sign,
the meaning of coefficients multiplying $\hat{f}_{t}$ could change across
expanding samples. We formed data sets to calculate the factors $\hat{f}_{t}$
for each of the three methods as follows. (i) For the explanatory variable $
\hat{f}_{t}^{MN},$ for the first two evaluation samples we followed \cite
{mccracken2016fred} in calculating the factor $\hat{f}_{t}^{MN}$ using the
full historical vintage of the FRED-MD database available as of 2015:4.
\footnote{
Our series for $\hat{f}_{t}^{MN1}$ for this subsample is almost (but not
quite) identical to the series analyzed by \cite{mccracken2016fred}. We have
been unable to identify the source of the small discrepancies.} For the
third evaluation sample, we estimated $\hat{f}_{t}^{MN}$ using the full
database available as of 2024:12. (ii) For the explanatory variable $\hat{f}
_{t}^{CF}$ we used the value of the CFNAI as it was reported in 2024 for all
historical dates.\footnote{
We downloaded $\hat{f}_{t}^{CF}$ on January 13, 2025 from the FRED database
at https://\allowbreak fred.stlouisfed.org/series/CFNAI.} (iii) For the
explanatory variable $\hat{f}_{t}^{HMX},$ we need to specify both (a) the
data set used to estimate the coefficients $\hat{{\Greekmath 010B}}_{i}$ for the
detrending regression (\ref{eqn_recommended_regression}) and (b) the data
set used to calculate principal components of the forecast errors $\hat{c}
_{it}=y_{it}-\hat{{\Greekmath 010B}}_{i}^{\prime }z_{i,t-h}.$ For (a), for all three
evaluation samples, we estimated $\hat{{\Greekmath 010B}}_{i}$ from the 2015:4 vintage
data set. For (b), for the first two evaluation samples, we calculated
principal components of $\hat{c}_{it}$ using the 2015:4 vintage data set.
For the last evaluation sample, we calculated principal components of $\hat{c
}_{it}$ using the 2024:12 data set. Thus our forecasting exercise for the
third subsample is not affected by the potential concern that the detrending
regression used information that was not available at the time the forecast
was made. This is because the coefficients for the detrending regression
were all known as of 2015. We also repeated the exercise using instead for
(b) trend coefficients estimated using all the data through 2024. This
alternative specification produced similar results (not reported here).
\begin{table}[tbph]
\caption{Mean squared forecast errors for different models}
\label{tab_mse}
\begin{tabular}{cccccccc}
\hline
& & \multicolumn{6}{c}{Consumer Price Index} \\ \cline{3-8}
sample & horizon & AR & CF & MN1 & HMX1 & MN2 & HMX2 \\ \hline
& & & & & & & \\
1970-1996 & h=1 & 7.91 & 1.00 & 0.99 & 1.03 & 0.98 & 0.90 \\
& h=6 & 4.26 & 0.81 & 0.77 & 0.80 & 0.95 & 0.88 \\
& h=12 & 5.32 & 0.70 & 0.62 & 0.74 & 1.02 & 1.33 \\
& & & & & & & \\
1997-2014 & h=1 & 12.26 & 1.03 & 1.04 & 1.02 & 0.97 & 1.09 \\
& h=6 & 6.08 & 1.23 & 1.23 & 1.23 & 0.96 & 1.11 \\
& h=12 & 4.21 & 1.22 & 1.22 & 1.28 & 0.95 & 1.17 \\
& & & & & & & \\
2015-2024 & h=1 & 8.07 & 1.51 & 1.27 & 1.14 & 0.94 & 1.04 \\
& h=6 & 3.70 & 1.95 & 1.48 & 1.76 & 1.03 & 1.01 \\
& h=12 & 3.64 & 1.71 & 1.30 & 1.49 & 1.09 & 0.96 \\
& & & & & & & \\ \hline
& & \multicolumn{6}{c}{Industrial Production} \\ \cline{3-8}
sample & horizon & AR & CF & MN1 & HMX1 & MN2 & HMX2 \\ \hline
& & & & & & & \\
1970-1996 & h=1 & 76.73 & 0.97 & 0.94 & 0.96 & 1.01 & 1.02 \\
& h=6 & 38.66 & 0.92 & 0.93 & 0.83 & 0.71 & 0.79 \\
& h=12 & 27.19 & 1.01 & 1.06 & 1.21 & 0.49 & 0.87 \\
& & & & & & & \\
1997-2014 & h=1 & 58.90 & 0.85 & 0.83 & 0.98 & 1.04 & 1.00 \\
& h=6 & 22.61 & 0.93 & 0.94 & 1.05 & 1.24 & 1.12 \\
& h=12 & 20.11 & 0.96 & 1.01 & 1.06 & 1.26 & 1.11 \\
& & & & & & & \\
2015-2024 & h=1 & 506.67 & 1.73 & 0.93 & 1.10 & 1.02 & 1.00 \\
& h=6 & 87.37 & 2.29 & 1.16 & 1.03 & 1.05 & 0.94 \\
& h=12 & 38.10 & 2.31 & 1.14 & 0.85 & 1.03 & 0.93 \\
& & & & & & & \\ \hline\hline
\end{tabular}
\vspace{5pt} \newline
Notes to Table \ref{tab_mse}. AR column reports simulated out-of-sample mean
squared forecast error for purely autoregressive model evaluated over three
different out-of-sample periods. CF column reports the relative MSE when
lags of the Chicago Fed National Activity Index are added to the
autoregression, with a value less than one indicating the variable is useful
for forecasting. MN1 column reports the relative MSE when lags of the first
principal component calculated using the procedures in \cite
{mccracken2016fred} are used in place of the CFNAI. HMX1 reports relative
MSE when lags of the first principal component of the estimated cyclical
components are used in place of CFNAI. MN2 and HMX2 report results when the
second principal component is used instead of the first.
\end{table}
For our first set of evaluations we set $y_{t+h}^{h}$ to be the average
inflation rate between month $t$ and $t+h,$ quoted at an annual rate,
\begin{equation*}
y_{t+h}^{h}=(1200/h)\log (CPI_{t+h}/CPI_{t}),
\end{equation*}
where $CPI_{t}$ denotes the level of the consumer price index in month $t.$
\footnote{
Here again our evaluation design differs slightly from that in \cite
{stock1999forecasting} and \cite{mccracken2016fred} in that those authors
took the object of interest to be to forecast the \textit{change }in the
inflation rate as a function of lagged changes:
\begin{equation*}
y_{t+h}^{h}=(1200/h)\log (CPI_{t+h}/CPI_{t})-1200\log (CPI_{t}/CPI_{t-1}).
\end{equation*}
} The column labeled $AR$ in Table \ref{tab_mse} reports the simulated
out-of-sample mean squared error of a purely autoregressive model for each
of the three evaluation samples and for forecast horizons of $h=1,$ 6, or 12
months.\footnote{\cite{stock1999forecasting} reported results for $h=12$
months whereas \cite{mccracken2016fred} reported results for $h=1,$ 6, and
12.} The first panel reproduces the finding of \cite{stock1999forecasting}
that an index like CFNAI significantly improves forecasts for longer
horizons over the 1970-1996 period. The alternative measures $\hat{f}
_{t}^{MN1}$ or $\hat{f}_{t}^{HMX1}$ offer similar improvements. The indexes
offer little or no improvement for one-month-ahead forecasts over this
period, but again are similar to each other. All three indexes are
outperformed by simple autoregressive forecasts over either of the later two
evaluation periods. The observation that inflation has become much harder to
forecast in data since 1996 has been reported by a number of other
researchers, including \cite{Atkeson}, \cite{Fisher}, \cite{StockWatson2007}
, and \cite{StockWatson2008}. Interestingly, the second cyclical factor $
\hat{f}_{t}^{HMX2}$ does better than any of the other four indexes at
forecasting inflation at the one-month horizon for the 1970-1996 sample and
at the 6-12 month horizons for the 2015-2024 sample.
Table \ref{tab_mse} also reports forecasts of industrial production, setting
\begin{equation*}
y_{t+h}^{h}=(1200/h)\log (IP_{t+h}/IP_{t}),
\end{equation*}
for $IP_{t}$ the level of the industrial production index in month $t.$ The
indexes $\hat{f}_{t}^{CF},$ $\hat{f}_{t}^{MN1},$ and $\hat{f}_{t}^{HMX1}$
all help forecast industrial production over near horizons in the first two
evaluation periods. The CFNAI does particularly poorly at forecasting either
inflation or industrial production at any horizon for 2015-2024. Both $\hat{f
}_{t}^{HMX1}$ and $\hat{f}_{t}^{HMX2}$ do significantly better than CFNAI in
every case over this period.
We conclude that our approach offers similar benefits to conventional PCA
when evaluated in terms of simulated out-of-sample forecasts, and does much
better than measures like the CFNAI for recent data. We share the conclusion
of the earlier literature that the usefulness of any
principal-component-based measure for purposes of forecasting depends on the
variable, evaluation period, and horizon of the forecast.{}
\section{Conclusion}
Calculating principal components of medium-horizon forecast errors is a
viable approach to identifying the common cyclical factors that drive a
large collection of potentially nonstationary economic indicators. This
avoids the need to decide how to detrend each individual series and is much
more promising than approaches such as the Chicago Fed National Activity
Index for handling data that include the large outliers of 2020.
\newpage
\bibliographystyle{apalike}
\bibliography{HXrefs.bib}