EconBase
← Back to paper

Diffusion index forecasts under weaker loadings: PCA, ridge regression, and random projections

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

55,025 characters · 16 sections · 54 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Diffusion index forecasts under weaker loadings: PCA, ridge regression, and random projections

\spacingset{1.8}

abstractWe study the accuracy of forecasts in the diffusion index forecast model with possibly weak loadings. The default option to construct forecasts is to estimate the factors through principal component analysis (PCA) on the available predictor matrix, and use the estimated factors to forecast the outcome variable. Alternatively, we can directly relate the outcome variable to the predictors through either ridge regression or random projections. We establish that forecasts based on PCA, ridge regression and random projections are consistent for the conditional mean under the same assumptions on the strength of the loadings. However, under weaker loadings the convergence rate is lower for ridge and random projections if the time dimension is small relative to the cross-section dimension. We assess the relevance of these findings in an empirical setting by comparing relative forecast accuracy for monthly macroeconomic and financial variables using different window sizes. The findings support the theoretical results, and at the same time show that regularization-based procedures may be more robust in settings not covered by the developed theory.\\ \noindentKeywords: diffusion index forecasting, principal component analysis, ridge regression, random projections.\\ JEL codes: C22, C38, C53.

Introduction

The diffusion index forecasting model developed by stock1998diffusion,stock2002forecasting provides a parsimonious approach to use information from large macroeconomic data sets such as the monthly FRED-MD and quarterly FRED-QD mccracken2015fred,mccracken2020fred. To establish consistency of diffusion index forecasts for the conditional mean, the classical assumption is that the factors are `pervasive' in the sense that these factors drive most variation in the available predictors. When there are $r$ factors, this pervasiveness is reflected in a sharp drop in the eigenvalues of the predictors' covariance matrix after the $r$th largest eigenvalue. In the absence of such a drop, PCA fails to consistently estimate the factors onatski2012asymptotics, and hence, yields inconsistent forecasts.

Recently, bai2023approximate provide an analysis of PCA-based factor estimation in the intermediate regime where the drop in the eigenvalues is `small'. In this paper, we utilize their results to analyze the consistency of diffusion index forecasts across three different estimation methods. First, we examine forecasts based on factors estimated by PCA. When the factor loadings are strong, the results naturally coincide with those in stock2002forecasting. Under weaker loadings, we show that the convergence rate depends on the loadings' strength as well as on the relative growth of the cross-sectional dimension ($N$) and the time dimension ($T$). The convergence rate tends to drop when $N$ increases faster than $T$. We also show that a substantial improvement in the convergence rate can be obtained when the idiosyncratic errors in the factor model are temporally independent.

Next, we explore two alternative regularization-based forecast methods: ridge regression forecasts and random projection based forecasts. The former were analyzed under a factor model structure by de2008forecasting as a form of Bayesian shrinkage. The latter can also be expected to work well when the data exhibits a factor structure boot2019forecasting. We provide a unified analysis of both methods by deriving bounds on the implicit regularization of the forecast induced by the choice of the ridge penalty and the random projection dimension. Both methods effectively soft-threshold the eigenvalues of the predictor matrix, unlike PCA-based methods that apply hard-thresholding.

We find that appropriate selection of regularization parameters ensures the consistency of both regularization-based forecasts for the conditional mean under the same conditions as the PCA-based forecast. Under strong loadings, the results improve the convergence rates for ridge regression forecasts established by de2008forecasting. We also show that the convergence rate can be negatively impacted by regularization bias under weaker loadings, particularly when the number of predictors grows faster than the sample size. In this case, the convergence rate can be slower relative to PCA.

We provide two sets of simulations. The crucial difference between these simulations is the presence of serial correlation in the idiosyncratic errors. Without serial correlation, the simulations are in line with the asymptotic results even for small samples. If there is autocorrelation in the idiosyncratic errors, then for weak loadings and small sample sizes, the accuracy of PCA seems to falter, while that of the regularization methods stays remarkably stable. Part of the deterioration of PCA can be explained by the fact that in small samples the nominal number of factors may be lower than the optimal number of factors. When we increase the dimensions of the data, the accuracy of PCA again improves over ridge regression and random projections.

The simulations tie in with our empirical application on forecasting monthly and quarterly macroeconomic and financial series from the FRED-MD and FRED-QD mccracken2015fred,mccracken2020fred. We estimate a diffusion index forecast model and manipulate the ratio of the sample size of the number of predictors in two ways while keeping the number of variables fixed: by choosing different estimation windows and by comparing across frequencies. In line with the developed theory, we find that as we decrease the moving window size, PCA forecasts become relatively more accurate compared to the regularization methods. At the same time, as in boot2019forecasting, and our simulations results for small samples with serial correlation in the idiosyncratic errors, random projections and ridge frequently outperform the PCA-based forecast. PCA yields better forecasts in about 10 percent of the cases when $T>N$, while PCA-based forecasts are more accurate than ridge and random projections for close to 40 percent of the variables when $T=\tfrac{1}{3}N$. In addition, in the monthly dataset, PCA-based forecasts are more accurate in 7 to 13 percent of the cases, but this increases to 19 percent when using quarterly data.

\paragraph{Related literature} The possible presence of weak loadings, for example in macroeconomic and financial applications, has triggered an active literature. onatski2012asymptotics considers an asymptotic regime where the eigenvalues are not well-separated, as in the spiked covariance matrix model of johnstone2009consistency and paul2007asymptotics. There, the ratio of the leading eigenvalues of the covariance matrix of the predictors with the smaller eigenvalues is assumed to be bounded. This constitutes an even weaker setting relative to the one considered in bai2023approximate that we follow in this paper.

Weak loadings as defined in the current paper can arise through a sparse structure of the loading matrix. Estimation and inference in this setting is considered in uematsu2022estimation and uematsu2022inference. fan2020learning consider a setting with weak loadings that are not necessarily sparse. The approach is somewhat similar to the random projection forecast as they use a cross-sectional projection to estimate the factors. The results presented there require the weights of the projections to be bounded (Assumption 2.1.i in fan2020learning). Also, instead of using a single projection, we consider averaging over multiple random projections to eliminate additional variance from the fact that we project onto a random subspace. karabiyik2020forecasting consider simple averages of predictors within a pre-specified group to recover group-specific factors. In the context of the FRED-MD, the groups would be the predefined variable groups. The random projection forecast we consider, takes a randomly weighted average across the entire predictor set.

Initial convergence rates for ridge regression in factor models as in de2008forecasting and carrasco2016sample suggested that ridge forecasts may converge at a slower rate compared to PCA-based forecasts. de2024asymptotic show the asymptotic equivalence of PCA-based forecast and ridge regression in settings with possible weak factors. However, this requires additional assumptions that do not fit the factor model of bai2023approximate that we consider in the present paper. We find that the equivalence between ridge forecasts and PCA-based forecasts holds under serially uncorrelated errors and a cross-section dimension not much larger than the time dimension. However, in other cases the regularization-based methods attain a slower convergence rate.

In a variant of the diffusion index forecast model where the idiosyncratic noise in the factor equation enters the outcome equation through a random effects specification, he2024ridge shows the optimality of ridge regression. Importantly, the paper shows the validity of a bias-corrected cross-validation procedure to select the optimal penalty in practice. Our theoretical results determine the optimal penalty parameter up-to-scale, and in practice we set it based on historical forecast performance.

Whereas ridge regression has a long history in forecasting applications, the use of random projections is relatively recent. Random projections are generally motivated from the Johnson-Lindenstrauss lemma johnson1984extensions, which can be used to show that they maintain a substantial fraction of the information, see for instance achlioptas2003database. Random projections can be applied to discrete choice models chiong2016random, to forecast product sales schneider2016forecasting, to reduce the dimension of vector autoregressions koop2016bayesian, to macroeconomic forecasting boot2019forecasting, and in high-dimensional hypothesis testing liu2023random. The equivalence to ridge regression, in terms of the convergence rate, appears new.

\paragraph{Notation} Denote $\delta_{NT} = \min(\sqrt{N},\sqrt{T})$ and the $j$th eigenvalue of a matrix $\bs A'\bs A$ as $\mu_{j}(\bs A'\bs A)$. $\mathcal{O}(T)$ denotes the space of orthogonal $T\times T$ matrices. Let, $||\bs A||^2=\text{tr}(\bs A'\bs A)$ denote the Frobenius norm. The vector $\bs b_{t}$ is the standard basis vector with a 1 in position $t$. Throughout, $M$ denotes a generic positive constant independent of any of the involved dimensions. For random variables $X$ and $Y$, we write $X\overset{(d)}{=}Y$ if $X$ and $Y$ have the same distribution. We write $T\asymp N$ if both $T/N=O(1)$ and $N/T = O(1)$.

Model and forecasts

We consider the diffusion index forecast model by stock2002forecasting,stock2002macroeconomic where a small number of factors $\bs f_{t}\in \mathbb{R}^{r}$ drives the variable of interest, and the same factors underlie a high-dimensional set of predictors $\bs x_{t}\in \mathbb{R}^{N}$,

eqnarray[eqnarray omitted — 184 chars of source]

Let $t=0,\ldots,T$, $\bs X=(\bs x_{0},\ldots,\bs x_{T-h})'$, $\bs F=(\bs f_{0},\ldots,\bs f_{T-h})'$, $\bs \Lambda = (\bs \lambda_{1}, \ldots,\bs\lambda_{N})'$, and $\bs y = (y_{h}, \ldots, y_{T})'$. One can add lagged outcome variables to (ref), but we omit those from the analysis to keep the discussion focused on the effects of a weak factor structure. Throughout, we write $\bar{y}_{t+h}$ for the conditional mean of (ref), i.e. $\bar{y}_{t+h} = \bs f_{t}'\bs\gamma$.

We follow bai2023approximate in defining the rescaled matrix $\bs Z = \bs{X}/\sqrt{NT}$ that has the singular value decomposition

equation[equation omitted — 123 chars of source]

Here $\bs{U}\in \mathcal{O}(T)$, $\bs{D}$ is $T\times N$ with the $\delta_{NT}^2$ nonzero singular values of $\bs Z$ on the diagonal of its upper left $\delta_{NT}^2\times \delta_{NT}^2$ submatrix, and $\bs{V}\in\mathcal{O}(N)$. Moreover, $\bs D_{r}\in\mathbb{R}^{r\times r}$ is a diagonal matrix with the $r$ largest singular values on its diagonal, the columns of $\bs U_{r}\in\mathbb{R}^{T\times r}$ and $\bs V_{r}\in\mathbb{R}^{N\times r}$ contain the corresponding left- and right-singular vectors. The matrices $\bs U_{-r}$ and $\bs V_{-r}$ contain the remaining left- and right-singular vectors.

The standard assumption on the loadings accompanying (ref) is that ${\bs \Lambda'\bs\Lambda/N\rightarrow_{p}\bs\Sigma_{\Lambda}}$ where $\bs\Sigma_{\Lambda}$ is a positive definite matrix. As in bai2023approximate, this paper considers the case where $\bs \Lambda'\bs\Lambda/N^{\alpha}\rightarrow_{p}\bs\Sigma_{\Lambda}$ for some $\alpha\in (0,1]$. When $\alpha=1$, we say that the loadings are `strong', while when $\alpha<1$, we refer to the loadings as `weak'.

PCA-based forecasts

PCA-based forecasts rely on the following estimates of the factors and loadings,

equation[equation omitted — 100 chars of source]

Denote by $\hat{\bs\gamma}$ the least squares estimator of $\bs\gamma$ in (ref) with $\bs f_{t}$ replaced by $\tilde{\bs f}_{t}=\tilde{\bs F}'\bs b_{t}$. The resulting PCA-based forecast is given by,

equation[equation omitted — 96 chars of source]

where

equation[equation omitted — 239 chars of source]

stock2002forecasting Theorem 2 shows that $\tilde{\bs f}_{T}'\hat{\bs\gamma}-\bs f_{T}'\bs\gamma\rightarrow_{p}0$ and the limiting distribution was established in bai2006confidence. Under an approximate factor model with strong loadings, the convergence rate is $\delta_{NT}$.

\paragraph{Remark} There is a slight subtlety in the construction of the forecasts. Suppose that $\bs e_{t}$ is independent across $t$. In that case, in the final term of (ref), $\bs e_{T}$ is independent of the error in estimating the loadings. This independence can be used to derive an improved convergence rate. An alternative way to construct the PCA-based forecast is to perform an eigendecomposition on the matrix formed by $(\bs x_{0}',\ldots,\bs x_{T}')$. However, in (ref), this would lead to dependence between $\bs e_{T}$ and $\tilde{\bs\Lambda}$. With currently available results, this yields a slower convergence rate.

Regularization-based forecasts

When the number of predictors is large, we can regularize the (potentially ill-defined) inverse of $\bs X'\bs X$ by adding a ridge penalty hoerl1970ridge. This leads to the following forecast.

equation[equation omitted — 142 chars of source]

This forecast is identical to a Bayesian forecast with a $N(\bs 0,\frac{k}{NT} \bs I_{N})$ prior on the regression coefficients. de2008forecasting analyze the consistency of this forecast for the conditional mean $\bar{y}_{t+h}$ and find that it may converge at a slower rate compared to the PCA-based forecast (ref), even when loadings are strong. We will show in (ref) that the convergence rate of (ref) actually coincides with that of (ref) under strong loadings.

boot2019forecasting consider the following alternative forecasting model that is based on random projections of the available variables,

equation[equation omitted — 94 chars of source]

where $\bs R\in\mathbb{R}^{N\times k}$ and the element of $\bs R$ are independent standard normal random variables. We emphasize that (ref) is interpreted as a model that defines a forecast. The actual data generating process is always taken to be (ref) and (ref).

Denote by $\hat{\bs\theta}_{R}$ the least squares estimator of $\bs\theta_{R}$ in (ref) and denote the associated forecast as $\hat{y}_{T+h|T,R}^{rp}=\bs x_{T}'\bs R\hat{\bs\theta}_{R}$. The random projection introduces additional noise in the forecast, which can be reduced by averaging over multiple draws of the random matrix $\bs R$. The random projection forecast then approximates,

equation[equation omitted — 193 chars of source]

boot2019forecasting show that if we take $O(N\log N)$ draws of the projection matrix $\bs R$, the mean squared forecast error approaches that of $\hat{y}_{T+h|T}^{rp}$.

It is not a coincidence that the regularization in both (ref) and (ref) is determined by the parameter $k$, as it is scaled such that the same $k$ (up to a constant scaling factor) determines the convergence rate to the conditional mean. In the subsequent sections, we use the notation $\hat{y}^{reg}_{T+h|T}$ to refer to both the ridge regression forecast and the random projection forecast.

Regularizing the empirical eigenvalues

For all three forecasts, the leading term in the expansion of the forecasts that captures the conditional mean $\bs f_{T}'\bs\gamma$ can be written as $ \bs f_{T}'(\bs F'\bs F)^{-1}\bs F'\bs U\bs \Delta\bs U'\bs F\bs\gamma.$ Suppose that we know that the true number of factors is $r$, and the regularization parameter $k$ is selected as we derive in (ref). Then, $\bs \Delta$ is a diagonal matrix that satisfies the following.

center[center omitted — 238 chars of source]

This clarifies that PCA is hard-thresholding the eigenvalues, while ridge regression and random projections are soft-thresholding the eigenvalues. However, in the asymptotic limit where $(N,T)\rightarrow\infty$, ridge regression and random projection may converge to PCA under a suitably chosen of $k$. To establish this formally requires a more detailed analysis of other terms appearing in the expansion of the forecasts. We also see that for weaker loadings, the regularization parameter $k$ has to be set larger to capture the signal in the first $r$ factors. However, at the same time this leads to a slower decline to zero for the remaining eigenvalues as this part is independent of the loading strength $\alpha$. This will eventually lower the convergence rate under weak loadings.

Theory

We have the following assumptions.

assumption\begin{enumerate}[label=(\roman*)] • The errors in (ref) can be written as $\bs e = \bs R_{T}^{1/2}\bs a\bs G_{N}^{1/2},$ with $\bs R_{T}\in \mathbb{R}^{T\times T}$ and $\bs G_{N}\in \mathbb{R}^{N\times N}$ deterministic matrices, • Conditional on $\bs\Lambda$ and $\bs F$, $a_{it}$ is i.i.d.\ over $i$ and $t$, with $\mathbb{E}[a_{it}|\bs\Lambda,\bs F]=0$, $\mathbb{E}[a_{it}^4|\bs\Lambda,\bs F]\leq M$, • $\mu_{T}(\bs R_{T})\geq C$, $ \mu_{N}(\bs G_{N})\geq C$, • $\max_{t}\sum_{s=1}^{T}|\bs R_{T,ts}|\leq C$, $\max_{i}\sum_{j=1}^{N}|\bs G_{N,ij}|\leq C$. \end{enumerate}

This assumption is closely related to the assumptions used in onatski2010determining and ahn2013eigenvalue. Part (ref) slightly strengthens the usual assumption on the maximum eigenvalue of $\bs R_{T}$ and $\bs G_{N}$, so we can easily establish that (ref) implies Assumption A1 and A3 in bai2023approximate.

corollaryLet $\bs e_{t}$ be defined as before and denote by $\bs e_{(i)}$ the $i$th column of $\bs e$. Under (ref) we have the following. \begin{enumerate}[label=(\roman*)] • For all $t$, $\frac{1}{N\sqrt{T}}\|\bs e_{t}'\bs e'\| = O_{p}(\delta_{NT}^{-1})$ and for all $i$, $\frac{1}{T\sqrt{N}}\|\bs e_{(i)}'\bs e\| =O_{p}(\delta_{NT}^{-1})\|$, • For each $t$, $\mathbb{E}\|N^{-\alpha/2}\sum_{i=1}^{N}\bs\lambda_{i}e_{it}\|^2\leq C$ and $\frac{1}{NT}\bs e_{t}'\bs e'\bs F = O_{p}(\delta_{NT}^{-2})$, • For each $i$, $\mathbb{E}\|T^{-1/2}\sum_{t=1}^{T}\bs f_{t}e_{it}\|^2\leq M$ and $\frac{1}{N^{\alpha}T}\bs e_{(i)}'\bs e\bs\Lambda = O_{p}\left(\frac{1}{N^{\alpha}}\right)+O_{p}\left(\frac{1}{\sqrt{N^{\alpha}T}}\right)$, • $\bs\Lambda'\bs e'\bs F=\sum_{t=1}^{T}\sum_{i=1}^{N}\bs\lambda_{i}\bs{f}_{t}'e_{it}=O_{p}(\sqrt{N^{\alpha}T})$. \end{enumerate}

We continue with Assumption A2 from bai2023approximate.

assumption\begin{enumerate}[label=(\roman*)] • $\mathbb{E}[\|\bs f_{t}\|^4]\leq M$, $\plim_{T\to\infty} \frac{\bs F'\bs F}{T}=\bs\Sigma_{F}\succ 0$, • $\|\bs\lambda_{i}\|\leq M$, $\lim_{N\to\infty}\frac{\bs\Lambda'\bs\Lambda}{N^{\alpha}}=\bs\Sigma_{\Lambda}\succ 0$ for some $\alpha>0$ with $\alpha\in (0,1]$, • The eigenvalues of $\bs\Sigma_{\Lambda}\bs\Sigma_{F}$ are distinct. \end{enumerate}

We also repeat Assumption A4 from bai2023approximate that limits the expansion rate of the cross-section dimension relative to the time dimension.

assumptionAs $(N,T)\rightarrow\infty$, $\frac{N}{N^{\alpha}}\frac{1}{T}\rightarrow 0$, for $\alpha$ as before.

To derive results for the diffusion index model, we need to make assumptions on the regression errors $\varepsilon_{t}$, which are provided by the following assumption.

assumptionThe regression errors $\bs \varepsilon$ satisfies $\mathbb{E}[\bs\varepsilon|\bs F,\bs\Lambda,\bs e]=\bs 0$ and $\mathbb{E}[\bs\varepsilon\bs\varepsilon'|\bs F,\bs\Lambda,\bs e] = \bs\Sigma_{\varepsilon}$ with $\|\bs\Sigma_{\varepsilon}\|\leq M<\infty$ $a.s.$

Finally, several results that we derive can be improved under the following assumption that strengthens (ref) part (i) to ensure that the noise is independent over time.

assumptionIn (ref) part (i), set $\bs R_{T}=\bs I$.

Convergence rates for PCA-based forecasts

Under strong factors, stock2002forecasting show consistency of the PCA forecast (ref) for the conditional mean $\bar{y}_{T+h|T}=\bs f_{T}'\bs\gamma$ when $(N,T)\rightarrow\infty$. bai2006confidence further show that $\hat{y}_{T+h|T}^{pca}-\bar{y}_{T+h|T} = O_{p}(\delta_{NT}^{-1})$.

Under weaker loadings, we have the following result for the PCA forecast.

theoremUnder \crefrange{assA1}{ass:varepsilon}, the PCA-based forecast (ref) satisfies \begin{equation} \begin{split} \hat{y}_{T+h|T}^{pca} - \bar{y}_{T+h|T} &=O_{p}\left(\frac{1}{N^{\alpha/2}}\right)+ O_{p}\left(\frac{1}{T^{1/2}}\right)+ O_{p}\left(\frac{N^{3/2}}{N^{3\alpha/2}T}\right)\\ &\quad +O_{p}\left(\frac{N^{1/2}}{N^{3\alpha/2}}\right) + O_{p}\left(\frac{N^2}{N^{2\alpha}T^{3/2}}\right)+O_{p}\left(\frac{N}{N^{2\alpha}T^{1/2}}\right). \end{split} \end{equation}

The first two terms on the right-hand side reduce to the known results when $\alpha=1$. However, a number of additional terms appear. To interpret these, we analyze three cases of interest. In the first, the loadings are strong, so $\alpha=1$. In the second, the loadings are weaker $\alpha<1$, and $N\asymp T$. In the third setting, again $\alpha<1$, but now the time dimension is of lower order relative to the cross-section dimension in the sense that $T = O(N^{\gamma})$ with $\gamma<1$. Specialized to these cases, (ref) yields the following.

corollaryCase 1. Strong loadings: $\alpha=1$. \begin{equation} \hat{y}_{T+h|T}^{pca}- \bar{y}_{T+h|T} = O_{p}(\delta_{NT}^{-1}). \end{equation} Case 2. Weaker loadings: $\alpha<1$. $N\asymp T$. \begin{equation} \hat{y}_{T+h|T}^{pca} - \bar{y}_{T+h|T} = O_{p}\left(\frac{1}{N^{\alpha/2}}\right) + O_{p}\left(\frac{1}{N^{(3\alpha-1)/2}}\right) . \end{equation} Case 3. Weaker loadings: $\alpha<1$. $T=O(N^{\gamma})$, $0<\gamma<1$. \begin{equation} \begin{split} \hat{y}_{T+h|T}^{pca} - \bar{y}_{T+h|T}&= O_{p}\left(\frac{1}{N^{\alpha/2}}\right) + O_{p}\left(\frac{1}{N^{\gamma/2}}\right) +O_{p}\left(\frac{1}{N^{(3\alpha-3+2\gamma)/2}}\right). \end{split} \end{equation}

Case 1 shows that under strong loadings, we obtain the known result on the convergence rate of PCA-based forecast. Case 2 shows that under weaker loadings, an additional term enters that requires that $\alpha>1/3$ for the forecast to remain consistent for the conditional mean. This result is not surprising, as bai2023approximate already show that the same condition is required for (a rotation of) the factors $\bs F$ to be consistently estimated. We will later show that this restriction can be lifted under (ref). Under Case 3, we see an interaction between the loading strength parameter $\alpha$ and the parameter $\gamma$ that governs how much smaller $T$ is relative to $N$. For consistency of the forecasts for the conditional mean, a sufficient condition is now that $\alpha>1-\frac{2}{3}\gamma$. This indicates that as the cross-sectional dimension increases relative to the number of time periods, we require stronger loadings for the forecast to be consistent.

Convergence rates under regularization

The forecast accuracy of both the ridge-regularized forecast (ref) and the random projection forecast (ref), depends crucially on the choice of the regularization parameter $k$. The following theorem stipulates the conditions imposed on $k$ needed for consistency, as well as the resulting convergence rate.

theoremLet \crefrange{assA1}{ass:varepsilon} hold. If $k$ is such that as $(N,T)\rightarrow\infty$, $\delta_{NT}^{-2}k\rightarrow 0$ and $\sqrt{\frac{N}{N^{\alpha}}}k^{-1}\rightarrow 0$, then \begin{equation} \begin{split} \hat{y}_{T+h|T}^{reg} - \bar{y}_{T+h|T} &=O_{p}\left(\frac{1}{N^{\alpha/2}}\right)+ O_{p}\left(\frac{1}{T^{1/2}}\right)+ O_{p}\left(\frac{N^{3/2}}{N^{3\alpha/2}T}\right)\\ & +O_{p}\left(\frac{N^{1/2}}{N^{3\alpha/2}}\right) + O_{p}\left(\frac{N^2}{N^{2\alpha}T^{3/2}}\right)+O_{p}\left(\frac{N}{N^{2\alpha}T^{1/2}}\right)\\ & + O_{p}\left(\frac{N}{N^{\alpha}}k^{-1}\right)+ O_{p}\left(k \delta_{NT}^{-1}\left[\sqrt{\frac{N}{N^{\alpha}}\frac{1}{T}}+\sqrt{\frac{1}{N^{\alpha}}} + \frac{1}{N^{\alpha/4}T^{1/4}}\right]\right). \end{split} \end{equation}

The first two lines are identical to those obtained for the PCA-based forecast. The second line shows how the regularization parameter balances capturing the signal from the true factors while suppressing the noise coming from the idiosyncratic errors. Again, we specialize this result to three special cases considered before.

corollaryCase 1. Strong loadings: $\alpha=1$. Set the regularization parameter as $k=O((N\delta_{NT}^2)^{1/4})$. \begin{equation} \hat{y}_{T+h|T}^{reg}- \bar{y}_{T+h|T}= O_{p}(\delta_{NT}^{-1}). \end{equation} Case 2. Weaker loadings: $\alpha<1$. $N\asymp T$. Set the regularization parameter as $k=O(N^{(3-\alpha)/4})$. \begin{equation} \hat{y}_{T+h|T}^{reg}- \bar{y}_{T+h|T}= O_{p}\left(\frac{1}{N^{\alpha/2}}\right) + O_{p}\left(\frac{1}{N^{(3\alpha-1)/4}}\right) . \end{equation} Case 3. Weaker loadings: $\alpha<1$. $T=O(N^{\gamma})$, $0<\gamma<1$. Set the regularization parameter as $k=O(N^{(1+2\gamma-\alpha)/4})$. Then, \begin{equation} \begin{split} \hat{y}_{T+h|T}^{reg} - \bar{y}_{T+h|T} &= O_{p}\left(\frac{1}{N^{\alpha/2}}\right) + O_{p}\left(\frac{1}{N^{\gamma/2}}\right) +O_{p}\left(\frac{1}{N^{(3\alpha-3+2\gamma)/4}}\right). \end{split} \end{equation}
figure[figure omitted — 846 chars of source]

In Case 1, we obtain the convergence rate $\delta_{NT}^{-2}$. This is identical to PCA and improves the convergence rate presented in de2008forecasting for ridge regression. In Case 2 with weaker loadings and $N\asymp T$, (ref) shows that for consistency, it is again sufficient that $\alpha>1/3$. However, the second term under the optimal choice of the regularization parameter $k$ decays slower to zero compared to what is the case for PCA-based forecasts. (ref) plots the different convergence rates for the PCA-based forecast and the regularization-based forecasts as a function of the strength of the loadings.

In Case 3, when $T = N^{\gamma}$ for $\gamma<1$, we reach a similar conclusion as for Case 2: the conditions under which the regularization-based forecasts are consistent are the same as for PCA, but regularization results in a lower convergence rate. The right panel of (ref) graphically shows the difference in the convergence rate for $\gamma=\{0.1,0.5,0.9\}$.

Faster convergence rates under serially uncorrelated idiosyncratic errors

The results of the previous section show that in general, we expect the regularization methods to suffer under weaker loadings. However, we now show that a substantial improvement over the convergence rates established in (ref) and (ref) can be obtained by introducing (ref), under which the rows of the idiosyncratic error term are independent. In this case, we have the following.

theoremUnder \crefrange{assA1}{assA5}, the PCA-based forecast (ref) satisfies \begin{equation} \begin{split} \hat{y}_{T+h|T}^{pca} -\bar{y}_{T+h|T}&=O_{p}\left(\frac{1}{N^{\alpha/2}}\right)+ O_{p}\left(\frac{1}{T^{1/2}}\right). \end{split} \end{equation}

We see that the condition that $\alpha>1/3$ is no longer necessary. What is still required is (ref). Case 3, where $T = O(N^{\gamma})$, therefore requires $\alpha+\gamma>1$.

The results for ridge regression and random projection forecasts can be similarly improved to the following.

theoremUnder \crefrange{assA1}{assA5}, and setting $k$ such that as $(N,T)\rightarrow\infty$, $\delta_{NT}^{-2}k\rightarrow 0$ and $\sqrt{\frac{N}{N^{\alpha}}}k^{-1}\rightarrow 0$, the regularization-based forecast (ref) satisfies \begin{equation} \hat{y}_{T+h|T}^{reg} - \bar{y}_{T+h|T}=O_{p}\left(\frac{1}{N^{\alpha/2}}\right)+ O_{p}\left(\frac{1}{T^{1/2}}\right)+ O_{p}\left(\frac{N}{N^{\alpha}}k^{-1}\right)+ O_{p}\left(k\cdot \delta_{NT}^{-2}\right). \end{equation}

For Case 1 with $k = O(\delta_{NT})$ and Case 2 with $k=O(N^{1-\alpha/2})$, we now obtain the same result as presented for PCA in (ref),

equation[equation omitted — 173 chars of source]

However, in Case 3, when $T = N^{\gamma}$ for $0<\gamma<1$, and $k=O(N^{(1-\alpha+\gamma)/2})$, we obtain

equation[equation omitted — 239 chars of source]

The requirement for consistency is $\alpha+\gamma>1$ as it is for PCA. However, in some parts of the parameter space, especially if both $\alpha$ and $\gamma$ are substantially below 1, the convergence rate is lower compared to PCA. The convergence rates are shown in (ref) for different values of $\gamma$ and $\alpha$.

figure[figure omitted — 480 chars of source]

Simulations

We consider the model (ref) and (ref) with $r=2$ factors. The elements of the factors $\bs F$ are independent standard normal random variables. To tune the strength of the loadings, we generate $[\bs\Lambda]_{ij} = \sqrt{\frac{N^{\alpha}}{N}}Z_{ij}$ with $Z_{ij}$ standard normal independent over $i=1,\ldots,N$ and $j=1,\ldots,r$. The parameter $\alpha$ will be varied over $\{0.5,0.75,1\}$ where $\alpha=1$ corresponds to a setting with strong factors. We first consider a setting with the idiosyncratic noise generated as $\bs e_{ti}\sim N(0,r)$, independently across $t=1,\ldots,T$ and $i=1,\ldots,N$. We subsequently consider a setting with serial correlation in the idiosyncratic noise.

In the forecast equation (ref), the parameter $\bs \gamma=(1,1)'$ and $\varepsilon_{t}\sim N(0,r)$ independently over $t$. We consider one-step ahead forecasts and set $\varepsilon_{T+1}=0$ as this is a noise term that we cannot forecast. We generate PCA-based forecasts assuming that the number of factors is known. For ridge regression and random projections, we pick $k$ ex post to minimize the mean squared forecast error. We compare the forecasts in terms of their Mean Squared Forecast Error (MSFE) defined as $\mathbb{E}[(\hat{y}_{T+h|T}-\bar{y}_{T+h|T})^2]$. Finally, we vary $N$ and $T$ over the grid $\{100, 200, 300, 400, 500\}$. Results are averaged over 5,000 realizations of the data generating process.

figure[figure omitted — 935 chars of source]

Independent idiosyncratic noise

The left panel of (ref) shows the one-step ahead MSFE for PCA-based forecasts as well as ridge and random projection forecasts when we set $T=N$. The right panel of the figure displays the optimal choice of the regularization parameter for the latter two methods. We observe the following: regardless of the value of $\alpha$ and the values of $N$ and $T$, PCA-based forecasts offer the lowest MSFE. There is no observable difference between ridge and random projection forecasts in terms of the MSFE. We also see that as the loadings become weaker, the difference between PCA and the regularization methods appears to increase. With regard to the regularization parameter, we see that it increases with increasing $(N,T)$ and with weaker loadings.

figure[figure omitted — 1,021 chars of source]

(ref) shows the MSFE of ridge regression relative to PCA-based forecasts. We omit the random projection forecasts, because the accuracy is identical to the ridge regression forecasts. In the left panel of (ref), we fix $T$ and vary $N$ on the $x$-axis. In the right panel, we fix $N$ and vary $T$ on the $x$-axis. The upper panel has weak loadings with $\alpha=0.75$, while the lower panel has strong loadings with $\alpha=1$. The most important finding is that for strong loadings, the relative performance between regularization-based methods and PCA stabilizes as either $N$ or $T$ increases. This confirms that the convergence rates for these methods are identical under strong loadings. Under weaker loadings (upper panel), we see that the relative MSFE stabilizes with increasing $T$ (right panel), yet it continues to increase with larger $N$ (left panel). This confirms the theoretical results from (ref) that state that under weaker loadings, the regularization methods may attain a lower convergence rate when $N$ is increasing at a faster rate compared to $T$.

Serially correlated idiosyncratic errors

We change the setup by introducing serial correlation in the idiosyncratic component of the predictors. To be precise, we consider $E_{ti} = \rho E_{t-1,i} + v_{t-1,i}$ with $v_{ti}\sim N(0,1-\rho^2)$ independent over $t$ and $i$, $E_{1i}=v_{1i}\sim N(0,1)$ and $\rho=0.7$.

figure[figure omitted — 472 chars of source]

(ref) shows again the MSFE and the optimal choice for the regularization parameter, akin to (ref). There are two important observations. First, we see a noticeable deterioration in the performance of PCA-based forecasts under weak factors and small sample sizes. There is hardly any change in the MSFE for the ridge or random projection forecasts, and we only find a slightly smaller optimal value of $k$. The deterioration of PCA can be explained intuitively by the fact that the eigenvalues are no longer clearly separated from the largest eigenvalue of the idiosyncratic errors. Suppose that the three largest eigenvalues are of the same order of magnitude. Because the autocorrelation increases the maximum eigenvalue of the covariance matrix of the idiosyncratic noise, it may be the case that the second largest eigenvalue is associated with a singular vector that is not related to the unobserved factors, while the third largest eigenvalue is associated with the unobserved factors. In this case, PCA would miss the relevant signal as it hard thresholds the eigenvalues at the supposed true number of factors $r$. The regularization-based methods, on the other hand, will apply comparable soft-thresholding to the largest three eigenvalues, maintaining the majority of the signal and thereby improving over PCA.

(ref) also adds PCA(k), which selects the optimal number of factors ex post. We see that indeed for weak loadings, the optimal number of factors is much larger than the nominal number of factors. Moreover, choosing the optimal number of factors largely mitigates the deterioration of the PCA-based forecasts. For small sample sizes, and weak factors, the regularization-based methods continue to offer more accurate forecasts. The second important observation is however that when $\alpha=0.5$ and we use PCA(k), we see that the MSFE for PCA converges faster to zero relative to random projections and ridge. This confirms the theoretical results in (ref) Case 2 versus (ref) Case 2.

Application: forecast accuracy on monthly and quarterly macroeconomic data sets

To confirm the theoretical results in practice, we setup a forecasting exercise to assess the impact of the ratio $T/N$ on forecasts from PCA, ridge and random projections in an empirical setting. The theory suggests that as $T/N$ decreases, PCA becomes relatively more accurate. However, the simulations also show that in small samples, the regularization methods may yield more accurate forecasts. Our goal is therefore analyse prediction accuracy on variables from the FRED-MD mccracken2015fred and FRED-QD mccracken2020fred. These data sets are large cross-sections of macroeconomic and financial variables and have been studied using PCA and random projections boot2019forecasting. We select vintage 2024-01 for both data sets.

For each variable, forecasts are generated from the model

equation[equation omitted — 99 chars of source]

where the $N_w \times 1$ vector $\bs w_t$ are the first four lags of $y_t$ and an intercept, which are always included. The $N_x \times 1$ vector $\bs x_t$ are the fifth and sixth lag of the target variables, and all the other variables. The $N_x \times k$ matrix $\bs R$ with $k < N_x$ reduces the dimension of the regressors.

The sample runs from 1959M1 (or 1959Q1) up to the end of 2023. Forecasting starts in the first period of 1970, i.e. 1970M1 for the monthly set and 1970Q1 for quarterly data. Out-of-sample evaluation begins in 1980M1 or 1980Q1, after a 10-year burn-in to select the hyperparameters that minimize the squared forecast error. So, the number of factors, the subspace dimension, and the ridge penalty are optimized over an expanding window that starts in 1970. The maximum number of factors in PCA and the subspace dimension in RP are set to 50. Random projection forecasts are averaged over 1,000 draws. For ridge, we set up a grid for the log of the penalty parameter that runs from -14.7 to 15, with increments of 0.3, following boot2019forecasting.

We conduct two experiments to assess the impact of the ratio $T/N_x$, changing the length or frequency of the estimation sample while keeping the number of regressors, $N_x$, fixed.

figure[figure omitted — 635 chars of source]

Moving window size

First, we vary the size of the moving window, $T$, from $\tfrac{1}{3}N_x$ to $\tfrac{4}{3}N_x$. The advantage of comparing across moving windows is that it controls the ratio $T/N_x$ as it is fixed over the evaluation sample. Monthly data allows for a wider range of moving windows and a sufficient number of forecasts for evaluation.

(ref) presents the percentage of times that PCA-based forecasts achieve a smaller MSFE compared to ridge or random projection forecasts. The percentages are below 50 percent for all window sizes, indicating that ridge and random projections yield on average more accurate forecasts than PCA in this dataset. For shorter moving window lengths, the forecasts obtained by PCA improve in terms of MSFE relative to ridge and random projections. With a window length of $T= \tfrac{7}{6} N_x$, PCA yields better forecasts in about 10 percent of the cases, while PCA-based forecasts are more accurate than ridge and random projections for close to 40 percent of the variables when $T=\tfrac{1}{3}N_x$. The increase in relative performance between longer and shorter windows is not monotonic, but there is a clear pattern indicating that ridge and random projections yield relatively less accurate forecasts than PCA at shorter windows.

(ref) presents the distribution of Diebold-Mariano (DM) test statistics diebold1995comparing computed on the difference between the squared error loss from regularization- and PCA-based forecasts. Positive values correspond to more accurate forecasts by ridge or random projections. The distribution of DM statistics is shifted downwards for shorter windows relative to the distribution for longer windows. DM statistics are more often negative with shorter windows, see also (ref). When forecasting CPI: Durables, PCA even produces significantly more accurate forecasts than ridge or random projections with windows shorter than $T=\tfrac{5}{6}N_x$, while the DM statistics are positive with windows larger than $T=N_x$.

figure[figure omitted — 916 chars of source]

(ref) presents the difference in DM test statistics between the largest window ($T=\tfrac{4}{3}N_x$) and the shortest window ($T=\tfrac{1}{3}N_x$), for each of the variables in the FRED-MD. Positive values indicate a shift of relative improvement for PCA compared to ridge or random projections with the shortest window compared to the larger window. The shift in the distribution of DM test statistics is not driven by a particular set of variables because the change is positive for most variables. The shift is particularly large for variables in categories 1 (Output and income) and 7 (Prices) and results in a change in the preferred model. For example, the DM statistic from comparing ridge (random projections) with PCA in the case of IP: Final Products and Nonindustrial Supplies is 1.30 (1.31) with $T=\tfrac{4}{3}N_x$, while it is -1.40 (-1.63) with $T=\tfrac{4}{3}N_x$. As an exception, PCA perform relatively worse for variables in category 3 (Housing).

figure[figure omitted — 1,211 chars of source]

Frequency

As a second experiment, we compare the FRED-MD and FRED-QD data sets using an estimation window covering the same calendar period but at different frequencies, thus varying $T/N_x$. We employ a moving window of 10 years, such that the number of observations is $T=120$ for monthly and $T=40$ for quarterly data. We consider the subset of 102 variables that are in both the FRED-MD and the FRED-QD, see Appendix (ref).

Monthly observations are aggregated to the quarterly level and monthly forecasts are at a $h=3$ month horizon. The quarterly forecasts are at a $h=1$ quarter horizon. Hence, the observations and forecasts measure the same entity across the data sets, only at different frequencies.

figure[figure omitted — 962 chars of source]

In 19.6 percent of the cases, PCA has a smaller MSFE than ridge in the quarterly data set and only in 12.7 percent of the cases when using monthly data. Similarly, PCA leads to more accurate forecasts than random projections with 18.6 percent of the variables for quarterly and 5.9 percent for monthly data.

This corresponds to a shift in DM statistics between using data at the quarterly or monthly frequency, for both ridge and random projections, see (ref). Using quarterly data results in DM statistics that are smaller than when using monthly data, implying that forecasts are preferred more often when using quarterly data. The change of the DM statistic can be substantial and lead to a different conclusion on which model yields more accurate forecasts. For example, the DM statistic for ridge (random projections) versus PCA is 3.32 (3.24) based on monthly data and -2.23 (-2.10) on quarterly data when forecasting the Switzerland / U.S. Foreign Exchange Rate.

Hence, the empirical findings are consistent between both experiments. Ridge and random projections generally yield better forecasts in the FRED-MD and FRED-QD than PCA. However, when the estimation sample becomes smaller compared to the cross-sectional size, i.e. the fraction $T/N_x$ decreases, the accuracy of PCA-based forecasts improves relative to the accuracy of regularization-based forecasts.

Conclusion

We analyze the convergence rates of PCA-based forecasts in the diffusion index forecast model with weaker loadings. We compare these rates with forecasts based on ridge regression and random projections. We find that under strong loadings, all forecasts attain the same convergence rate. Under weaker loadings all forecasts are consistent under the same assumptions on the loading strength, but ridge regression and random projection forecasts converge more slowly relative to PCA when the cross-section dimension is large relative to the number of time periods. We confirm these results in simulations and an empirical study based on the FRED-MD and FRED-QD data sets. These numerical results also show that in small samples, regularization-based methods may be more robust to serial correlation in the idiosyncratic errors.

\ifthenelse{\boolean{BLIND}} { \paragraph{Funding details} Tom Boot acknowledges financial support by the Dutch Research Council (NWO) as part of grant VI.Veni.201E.11. }

\paragraph{Disclosure statement} There are no competing interests to declare.

\ifthenelse{\boolean{APPON}}{