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.
114,115 characters · 13 sections · 155 citation commands
Robust Tests for Factor-Augmented Regressions with an Application to the Novel EA-MD-QD Dataset
{0.5cm} { \affil[1]{University of Bologna} \affil[2]{BI Norwegian Business School} }
JEL Classification: C12; C13; C33.
Keywords: Forecasting; factor-augmented regression model; common factor model; common correlated effects. \onehalfspacing
Nowadays, forecasting by using so-called diffusion indexes has become increasingly popular thanks to the growing availability of large datasets. In fact, this parsimonious approach allows forecasters to extract the predictive content of many potential predictors into a reduced number of indexes, also referred to as latent factors, when the data follows an approximate factor structure proposed by chamberlain1983funds. Because the factors are unobserved, they must be estimated in a first step before augmenting the smaller forecasting model (possibly of autoregressive nature). The dominant method to estimate the factors is Principal Components (PC) (see e.g. bai2002determining). The simple yet powerful motivation behind this setup is that the estimated factors have been shown to enhance the accuracy of economic forecasts empirically (see stock1999forecasting; forni2003financial; Ludvigson_Ng_2011 or ciccarelli2010global from an extensive literature). Although the rich content of predictive information embedded in the factors represents an opportunity to make more accurate forecasts, it is essential to formally evaluate their predictive power. Fundamentally, the approach entails (pseudo) out-of-sample forecast evaluation between the factor-augmented model and its nested baseline. \\ In this article, we build on recent advances in the forecasting literature: both methodological and data advancements. We develop novel tests of forecast accuracy and encompassing for nested regressions with latent factors. Our contribution is to extend the innovative and extremely useful inferential theory for pseudo out-of-sample comparisons of pitarakis2023direct,pitarakis2025novel to the case of generated regressors, where asymptotic normality of the proposed statistics holds under many empirical regularities: 1) uncertainty about the number of factors, 2) presence of structural breaks, or 3) persistent series. We achieve this by revisiting the methodology by karabiyik2021forecasting and stauskas2022tests who were the first to abandon the usual PC and instead adopted cross-section averages (CAs) of grouped series - as first proposed by pesaran2006estimation - in a forecast comparison setup.\footnote{Note that pesaran2006estimation proposed the Common Correlated Effects (CCE) estimator for the parameters of the panel linear regression model, where CAs constitute only its first step to proxy factors. In our time series setting, we only rely on the CAs.} In contrast to extracting the factors as eigenvectors corresponding to the largest eigenvalues of the covariance matrix of the data, CAs average the noise and settle on the so-called low-rank component of the data, which are the factors identified up to a linear transformation. \\ The notion of CAs comes naturally in practice, as it is becoming progressively standard to organize the predictors into large blocks of variables, such as consumption, money aggregates, prices or exchange rates. In fact, eickmeier2008successful consider 52 studies that use factor-augmented regression models to forecast inflation and/or output, and the datasets used in these studies all have the same structure of blocks of variables. Two notable examples of such structure are the FRED-MD (monthly) database by mccracken2016fred and its quarterly counterpart FRED-QD (see mccracken2020fred). EAdataset recently presented their first European counterpart - namely the EA-MD-QD dataset - which covers the Euro Area as a whole and the main member countries with 100+ series organized in 11 blocks. As such datasets offer a large number of estimated factors and cover a long time horizon, a robust procedure to evaluate their predictive content becomes key. \\ The CAs procedure in relation to the properties of the statistics of pitarakis2023direct,pitarakis2025novel brings advantages relative to PC. First, as averages are taken at each point in time, factors can be integrated with their order being unknown, while such information is integral to the optimization procedure in PC (see bai2004estimating); second, CAs can be combined with structural breaks in factor loadings which can occur in the out-of-sample period, because a single break in loadings is algebraically equivalent to doubling the number of factors.; third, if the number of CAs is greater than the number of factors, controlling for overspecification in the parameter estimation step is simpler and does not require restrictions on the dimensions of the dataset or restrictive normality assumptions on model errors (see moon2015linear). Importantly, we do not require that a single block loads on a single factor only. Indeed, each block can load on all the factors in the system, or many blocks may not load on any factors at all, as long as the average loading matrix has a rank equal to the true number of factors. This relates to and incorporates the notion of local factors, which affect only a subset of series, in the PC setting (see freyaldenhoven2022factor). Finally, missing values are not a problem for CAs as they can still be easily computed in unbalanced panels (see chudik2015large). On the other hand, such case is difficult in PC as it should be coupled with an alternating maximization algorithm (see bai2009panel, or su2025estimation). \\ The studies of stauskas2022tests and gonccalves2017tests (in the PC setting) explore the statistics of clark2001tests and mccracken2007asymptotics, which solve the degeneracy problem faced by the seminal tests of west1996asymptotic in the nested environment (see also a survey in diebold2015comparing). However, the solution comes at the expense of highly non-standard limiting distributions based on functionals of stochastic integrals of Brownian motion and the relative growth rate of in-sample versus out-of-sample observations. Although asymptotically valid critical values can be obtained by simulation (see clark2012reality; or hansen2015equivalence), their practical implementation remains challenging. \\ The encompassing and forecast accuracy tests in pitarakis2023direct,pitarakis2025novel have two key features: 1) they put an end to the joint problem of variance degeneracy and nonstandardness of asymptotics in nested models, and asymptotically follow a standard normal distribution, and 2) the observed predictors can be non-stationary. Therefore, we naturally fill in the gap in the literature by allowing CAs-generated regressors. If combined with the framework of pitarakis2023direct,pitarakis2025novel, the mechanics of CAs solve many issues that arise in stauskas2022tests and gonccalves2017tests. In stark contrast to all existing tests of equal accuracy for factor-augmented regressions, our user-friendly theoretical framework can flexibly accommodate for: 1) different degrees of persistence in the regressors, 2) structural breaks in the factor loadings in the full sample and, most importantly, 3) overestimation of the number of factors. To the best of our knowledge, this is the first result in the CAs literature where an overestimation of the number of factors does not have any asymptotic effect, irrespective of the expansion rates of $N$ and $T$. In spirit, this result resembles moon2015linear in the PC context. Furthermore, extensive Monte Carlo experiments on size and local power exhibit excellent finite sample properties of our tests under a plethora of realistic economic scenarios, thus making our methodology relevant and applicable in many empirical applications. Ultimately, this study is the first to examine the predictive capacity of factors in the EA-MD-QD dataset by EAdataset in a nested setup. We show that they offer powerful predictive insight in the context of macroeconomic forecasting especially in the presence of large shocks such as the Covid pandemic. \\ We also highlight that this work goes in parallel to, and draws some interim results from, the PC-based theory of margaritella2024new for the same statistics. The latter study is the first in the PC setting, and it utilizes the fact that the method can be applied to any dataset with large dimensions. The current study, on the other hand, leverages the block structure that is natural in large macroeconomic datasets. This suggests a simpler estimator with empirically relevant properties, which are not typically possessed by PC. \\ The rest of the paper is organised as follows. Section 2 introduces our forecasting setup, the CCE approach and the tests under the observed predictors together with a list of assumptions. Section 3 dives into the CCE-based statistics and presents the main results. In Section 4, we conduct an extensive Monte Carlo study. The notation adopted in the paper is as follows. Scalars, vectors and matrices are denoted by $a \in \mathbb{R}$, $\*a\in \mathbb{R}^p$ and $\* A\in \mathbb{R}^{c \times p}$, respectively. For any matrix $\* A$, the Frobenius norm is defined as $\left\|\*A\right\| = \sqrt{\textrm{tr}(\* A'\* A)}$ where $\mathrm{tr}(.)$ is the trace operator. Moreover, $\mathrm{rk}(\*A)$ and $\*A^+$ represent rank and Moore-Penrose inverse, respectively, while $\lambda_{\mathrm{min}/\mathrm{max}}(\*A)$ is the minimum/maximum eigenvalue of $\*A$. Next, $\lfloor x \rfloor$ represents the integer part of $x$, $\sup$ denotes the supremum, $k_0$ represents the number of in-sample observations, while $T - k_0 = n$ is the number of the out-of-sample observations. With respect to the notion of convergence, $\to$ and $\to_p$ denote the limit and convergence in probability, respectively.
Consider a vector of predictors $\*z_t=(\*w_t',\*f_t')'\in \mathbb{R}^{q+r}$ such that we have the following forecasting models:
where $\*w_t$ represents a small set of known predictors. While we may choose the observables $\*w_t$ based on economic theory or our experience as forecasters, it is important to remember that nowadays we face an abundance of potential predictors that may or may not improve our forecasts. In order to exploit the information embedded in large dimensional datasets, one could extract the potential predictive content of the many regressors into just a few series of latent factors $\*f_t$. This situation is reflected in ((ref)). Because $\*f_t$ is typically unobserved, we define the infeasible least squares (LS) estimator of $\+\delta$ for each $t=1,\ldots,T$ as
where the out-of-sample forecast error is defined as $\widetilde{u}_{2,t+1}=y_{t+1}-\widetilde{\+\delta}_t'\*z_t$ ($\widetilde{\+\theta}_t$ and $\widetilde{u}_{1,t+1}$ are defined analogously). Throughout, $\*z_t$ can be stationary or non-stationary. We are not interested in specific time series properties of $\*z_t$, but we want to vary its degree of persistence and its effect on the forecast comparison problem. Large databases, such as FRED-MD or EA-QD-MD, offer “recipes” for how to standardise the long series (see e.g. Appendix A in EAdataset). Some persistence may remain and this can be flexibly modeled by moderate integration in the spirit of magdalinos2009limit, which we formalize in our assumptions. \\ To extract $\*f_t$, we assume that our panel of $M$ predictors can be divided into a finite number $m$ of blocks exogenously. The economic rationale underlying this approach lies in the idea that variables with similar economic features are expected to be driven fundamentally by a series of (possibly unobserved) common shocks (see e.g. hallin2011dynamic; moench2013dynamic; ando2016panel) and, as a result, are to be grouped into one block. The FRED-MD database of the Federal Reserve (see Ludvigson_Ng_2011 and mccracken2016fred for more details) is among the primary examples, where the series are grouped in 8 blocks. The new EA-MD-QD dataset by EAdataset contains 11 blocks: (1) National Accounts, (2) Labor Market Indicators, (3) Credit aggregates, (4) Labor Costs, (5) Exchange Rates, (6) Interest Rates (7) Industrial Production and Turnover, (8) Prices, (9) Confidence Indicators, (10) Monetary Aggregates and (11) Others. It is then a reasonable to assume that every series in each block is, due to their very characteristics, sensitive to the same changes in the economic environment. \\ To reflect such data configuration formally, let $I_j = \{1_j, \dots, N_j\}$ be the set of indexes of the series contained in each block for $j = 1, \dots, m$. We denote the predictor panel data variable as $\* x_{i,t} = [x_{i_1,t}, \dots, x_{i_m,t}]^\prime \in \mathbb{R}^m$, where $x_{i_j,t}$ can be seen as the predictor $i$ of the block $j$. That is, we can allow for blocks of different sizes. For simplicity of notation and without loss of generality, let us assume that every block has the same number of predictors $N = M/m$. The data generating process of $\*x_{i,t}$ is assumed to be
where $\* f_t \in \mathbb{R}^r$ is a vector of common factors, $\boldsymbol{\Lambda}_i \in \mathbb{R}^{r \times m}$ is the matrix of factor loadings, and $\* e_{i,t}\in \mathbb{R}^{m}$ are the idiosyncratic components. In this work, we estimate the latent factors using CAs as suggested in pesaran2006estimation:
where $\overline{\* A} = N^{-1} \sum_{i=1}^N \* A_i$ denotes an average for any matrix $\* A_i$. Note that for every fixed $t=1,\ldots, T$, we have $\left\|\overline{\*e}_t\right\|=O_p(N^{-1/2})$ under many empirically relevant assumptions (see e.g. pesaran2011large), and so the factors are identified up to a linear transformation, which is $\overline{\+\Lambda}$. Therefore, given that $\mathrm{rk}(\overline{\+\Lambda})=r$, we obtain
which justifies using averages to approximate the factors up to a linear transformation. Clearly, if $m=r$, then $\overline{\+\Lambda}^+=\overline{\+\Lambda}^{-1}$, which is equivalent to knowing the number of factors $r$. Therefore, the overall set of feasible predictors is $\widehat{\*z}_t=(\*w_t',\widehat{\*f}_t')'\in \mathbb{R}^{q+m}$, meaning that the feasible estimator is
where $\widehat{u}_{2,t+1}=y_{t+1}-\widehat{\+\delta}_t'\widehat{\*z}_t$. Notice how ((ref)) does not restrict one factor to one block, but it is allowed. Similarly, each block can load on all $r$ factors. In fact, we can have $\overline{\+\Lambda}=\left[\overline{\+\Lambda}_r,\+0_{r\times (m-r)} \right]$, so that $m-r$ blocks are driven purely by idiosyncratics. As long as $\mathrm{rk}(\overline{\+\Lambda}_r)=r$, ((ref)) holds and we have $r$ “local” factors that underlie only a subset of blocks. Therefore, CAs seamlessly accommodate the notion of local factors from the PC setting (see freyaldenhoven2022factor). \\ Note that $\left\|\overline{\*e}_t\right\|=O_p(N^{-1/2})$ implies that $\frac{1}{T}\sum_{t=1}^T\widehat{\*f}_t\widehat{\*f}_t'$ is asymptotically singular $\overline{\+\Lambda}'\left(\frac{1}{T}\sum_{t=1}^T\*f_t\*f_t'\right)\overline{\+\Lambda}$ is, unless $m=r$. Consider the case when $m>r$, which is equivalent to overspecifying the number of factors. Indeed, based on the discussion above, EA-MD-QD would provide us with $m=11$ estimated factors ($m=8$ in FRED-MD), while the average number detected in the literature is smaller.\footnote{For instance, in the case of FRED-MD, stauskas2022tests estimate $r=2$ by using the criterion of ahn2013eigenvalue. margaritella2024new estimate 7 on average by using less conservative criteria from bai2002determining.} In order to proceed with the analysis of ((ref)), we need to re-define our target object since we estimate $\*f^0_t=[\*f_t', \*0_{1\times(m-r)}]'$. Hence, similarly to stauskas2022tests, we introduce a so-called rotation matrix $\overline{\*H}\in \mathbb{R}^{m\times m}$, such that
Next, we let $\*D_N=\mathrm{diag}(\*I_r, \sqrt{N}\*I_{m-r})$, such that
The transformed error is $\overline{\mathbf{e}}^0_t=[\overline{\mathbf{e}}^{0\prime}_{r,t}, \overline{\mathbf{e}}^{0\prime}_{-r,t}]'$, where $\left\|\overline{\mathbf{e}}^{0}_{r,t} \right\|=O_p(N^{-1/2})$, but $\left\|\overline{\mathbf{e}}^{0}_{-r,t} \right\|=O_p(1)$. This means that $\left(\frac{1}{T}\sum_{t=1}^T\widehat{\*f}_t^{0}\widehat{\*f}_t^{0\prime}\right)^{-1}$ exists asymptotically. Note that if $m=r$, then $\overline{\*H}=\overline{\+\Lambda}^{-1}$ and $\*D_N=\*I_r$, and so ((ref)) reduces to $\*f_t+\overline{\+\Lambda}^{-1\prime}\overline{\*e}_t$. We leave the definition of $\overline{\*H}$ in the Supplement. \\ We stress that $\overline{\*H}$ serves in ((ref)) as a purely analytical device, which exists as long as $\overline{\+\Lambda}$ does. It provides means to track the contribution of the redundant $m-r$ averages on $\widehat{\+\delta}_t$ and, in turn, $\widehat{u}_{2,t+1}$, which is the key component in the tests. Importantly, we do not need to know or apply $\overline{\*H}$ in practice. Note that by applying the Frisch-Waugh-Lovell Theorem to the second coordinate of ((ref)) (in stacked notation, where we suppress dependence on $t$), we obtain the following:
where $\*M_\*A=\*I-\*A(\*A'\*A)^{-1}\*A'=\*I-\*A\*B(\*B'\*A'\*A\*B)^{-1}\*B'\*A'=\*M_{\*A\*B}$ is a projection matrix onto the orthogonal complement of $\*A$ or $\*A\*B$ for a full rank $\*B$ ($=\overline{\*H}\*D_N$). It is essential that projecting onto $\widehat{\*f}_t$ is equivalent to projecting onto $\widehat{\*f}_t^0$ since $\overline{\*H}\*D_N$ is a full-rank matrix. \\ In this step, it is natural to compare the current study with margaritella2024new in the PC setting. Although PC and CAs have rather distinct theoretical bases, CAs have technical advantages beyond those differences. Firstly, as seen from ((ref)), the overestimation of $r$ is embedded in our analysis. We pursue this route instead of assuming the existence of a consistent estimator of $r$ (the usual step in the PC setting; see e.g. bai2006confidence) to carefully characterize the asymptotic behavior of $m-r$ redundant CAs. Moreover, we effectively work with the predictors of the mixed integration order, which can be handled by CAs (see westerlund2018cce). Although $\*z_t$ is modeled as a mildly integrated process, $\overline{\*e}_{-r,t}^0$ is still stationary. Consequently, redundant $m-r$ CAs can be seen as additional predictors, whose asymptotic contribution needs to be handled very carefully given the overall non-stationary environment. Note that a mixed order of predictors is not allowed in the original studies of pitarakis2023direct and pitarakis2025novel, but we show that this is not necessary in our problem. In the next section, we introduce the tests under $\*z_t$, and then give their representations under $\widehat{\*z}_t$ in Section 3.
The goal of our forecast accuracy tests is to assess which of the two models is more informative in predicting the dependent variable: if $\+\alpha=\*0_{r \times 1}$, then $\*f_t$ does not have any predictive power towards $y_{t+1}$ and it is sufficient to make forecasts with the model in ((ref)). This also means that under the null of $\+\alpha=\*0_{r \times 1}$, the model in ((ref)) nests the ”smaller” model in ((ref)) so that $u_{1,t+1}=u_{2,t+1}$. Such null implies that $\mathbb{E}(u_{1,t+1}(u_{1,t+1}-u_{2,t+1}))=0$. The latter equality can be tested in numerous ways, where one of the most popular tests is given by clark2001tests, which relies on the estimated forecast errors $\widetilde{u}_{j,t+1}$ for $j\in \{1,2\}$. They are based on the recursive estimators $\widetilde{\+\theta}_t$ and $\widetilde{\+\delta}_t$. The test has a highly non-standard asymptotic distribution, which also depends on tuning parameters. pitarakis2023direct, on the other hand, proposed a modification that results in an asymptotically normal distribution free of nuisance/tuning parameters. To compare the forecasting performance of two models in a nested environment, we propose four new tests of forecast nesting and accuracy by extending the work of pitarakis2023direct and pitarakis2025novel to the case of predictive regressions with generated regressors estimated by CAs. Firstly, let us consider the test for forecast nesting from pitarakis2023direct. Assuming the factors are observed, the test statistic for forecast nesting from pitarakis2023direct is given by
where $m_0=\lfloor n \mu_0\rfloor = \lfloor (T-k_0) \mu_0\rfloor$ is a cut-off point to split the average for $\mu_0\in (0,1)$, $\mu_0\neq 1/2$, and $\widetilde{\omega}_1^2$ is the estimated variance of the limiting distribution. As shown in pitarakis2023direct, values of $\mu_0$ closer to 1/2 yield better local power in finite sample, so we set it equal to 0.45 hereafter. Still, the true factors are not observed and must be estimated instead, so the subscript ”$f$” indicates that the test is infeasible. Secondly, pitarakis2025novel proposed two additional tests of forecasting accuracy whose corresponding infeasible statistics are:
where $l^0_j=\lfloor n\lambda_j^0 \rfloor$ for $j=1,2$, with $\lambda _j^0\in (0,1)$ representing two portions of the out-of-sample period. That is, the forecast MSE loss differentials of both models are computed over partially overlapping out-of-sample segments, i.e. $l^0_1>l_2^0$ or vice-versa, and $\widetilde{\omega}_j^2$ for $j=2,3$ are the corresponding variance estimators. Note that if the two segments were fully overlapping then their variance would be asymptotically degenerate as in clark2001tests, so $l^0_1=l_2^0$ is ruled out to avoid incurring into this exact issue. Further, note from the formulation of ((ref)) that, once $l_2^0$ (or $\lambda_2^0$) is fixed, new test statistics can be obtained by averaging ((ref)) over some chosen feasible set of $l_1$ based on the tuning parameter $\nu_0\in (0,1)$. As shown by margaritella2024new, among other averaging possibilities is also a fourth statistic of the form
which is obtained by fixing $l_1^0$ (or $\lambda_1^0$) and averaging ((ref)) over some chosen feasible set of $l_2$. Notice that pitarakis2025novel provide recommendations for $l_1^0$, $l_2^0$ and $\nu_0$. Finally, margaritella2024new also show that $\widetilde{\omega}^2_j$ for $j \in \{1,2,3,4\}$ can be written as
up to a known constant specific to each statistic. In all four cases, we test the null against the one-sided alternative since $s_{f,j}\to_d \mathcal{N}(0,1)+c$ for $c>0$ that depends on the model parameters.
Throughout our analysis, we employ the following set of assumptions. \\
Assumption 1. $\mathbb{E}(u_{t+1}|\mathcal{F}_t)=0$, $\mathbb{E}(u^2_{t+1}|\mathcal{F}_t)=\sigma^2$, and $\mathbb{E}(u^4_{t+1}) < \infty$ for all $t$, where $\mathcal{F}_t$ is the sigma-algebra generated by $\{\mathbf{z}_t,\mathbf{z}_{t-1},...,\*X_t, \*X_{t-1},\ldots, y_t, y_{t-1},\ldots\}$.
Assumption 2.
Assumption 3.
Assumption 4. $\mathbf{z}_t$ and $u_{t}$ are independent of $\mathbf{e}_{i,s}$ for all $t$, $s$ and $i$.
Assumption 5. $\boldsymbol{\Lambda}_i$ is a non-random matrix, such that $\overline{\boldsymbol{\Lambda}}\to \boldsymbol{\Lambda}$ as $N\to \infty$ and $\overline{\boldsymbol{\Lambda}}=[\overline{\boldsymbol{\Lambda}}_r, \overline{\boldsymbol{\Lambda}}_{-r}]$, where $\overline{\boldsymbol{\Lambda}}_{-r}\in \mathbb{R}^{r\times (m-r)}$ and $\overline{\boldsymbol{\Lambda}}_r\in\mathbb{R}^{r \times r}$ is full rank for all $N$, including $N\to \infty$, and $\|\boldsymbol{\Lambda}_i\|<\infty$. If $m=r$, then $\overline{\boldsymbol{\Lambda}}=\overline{\boldsymbol{\Lambda}}_r$.
Assumption 6. $\boldsymbol{\alpha} = T^{-1/4-\tau/2}\boldsymbol{\alpha}^0$, where $\|\boldsymbol{\alpha}^0\|<\infty$. \\
Assumption 1 is standard in the literature and imposes a martingale difference sequence (MDS). It implies one-step-ahead forecasts (see e.g. campbell2006efficient, hjalmarsson2010predicting, or breitung2015instrumental). While some tests in pitarakis2023direct, pitarakis2025novel and allow for $h$-step-ahead forecasts, we focus on the smallest set of conditions to accommodate all of them. However, we explore the theory behind $h$-steps, and explain the extension in Remark 1. The effect of general serial correlation in $u_{t}$ is explored in simulations, where the power remains for low levels of correlation. Assumption 2 (a) provides comparative statics, as $\tau \to 1$ increases the persistence of the process. magdalinos2009limit allow us to manipulate $\tau$ and give theoretical guarantees, such as $\frac{1}{T^{1+\tau}}\sum_{t=1}^T\*f_t\*f_t'\to_p\+\Sigma_{\*F}$ (constant positive definite matrix). The convenient feature of this specification is that the discussion below ((ref)) and ((ref)) still holds. Part (b) ensures that we return to the usual stationarity conditions under $\tau=0$. Part (c) regulates the behavior of a recursive LS estimator under mildly integrated predictors. \\ Assumption 3 (a) reveals that we cannot maintain an arbitrary correlation of idiosyncratics and persistent predictors simultaneously. In the literature of CAs and PC under $\tau=0$, we require that terms of the form of $\frac{1}{T}\sum_{t=1}^T\*z_t\*e_{i,t}'$ vanish for all $i$ in order to demonstrate that the factor estimation error is negligible. When $\tau >0$ and $\{\*e_{i,t}\}_{t=1}^T$ is serially correlated, the rate is too slow, and this problem will remain in the recursive estimation. We have that $\left\|\frac{1}{T^{1+\tau}}\sum_{t=1}^T\*z_t\*e_{i,t}'\right\|=O_p(T^{-\tau/2})$ (see Lemma 3.1 in magdalinos2009limit), but if $\{\*e_{i,t}\}_{t=1}^T$ follows an MDS process or a moving average (MA) process of finite order, then $\left\|\frac{1}{T^{1+\tau}}\sum_{t=1}^T\*z_t\*e_{i,t}'\right\|=O_p(T^{-(1+\tau)/2})$ as needed (see our Auxiliary Results in the Supplementary material). Under Assumption 2 (b), the rate is suitable even if the idiosyncratics exhibit more general time-dependence. Overall, this can be seen as a price for conducting “comparative statics” by altering $\tau$, as this restriction is inherent in the process of magdalinos2009limit.\footnote{In the Supplement, we provide additional Monte Carlo evidence which reveal that correlation in $\*e_{i,t}$ is not really harmful for the statistical power of the tests.}\\ Assumption 3 (a) part (ii) is similar to the one in stauskas2022tests, and it regulates the dependence under stationarity so that $\{ \+\xi_s\}_{s=1}^t$ obeys the invariance principle and is needed to make sure that (the feasible) $\widehat{\+\delta_t}$ does not dominate the asymptotic theory. Assumption 3 (b) allows for weak cross-section dependence together with unconditional time and cross-section heteroskedasticity. For example, a spatial dependence structure employed in stauskas2022tests is a special case. Assumption 4 is standard, but it can be relaxed at the expense of higher moment requirements, while Assumption 5 treats the factor loadings as fixed parameters. They can be made stochastic as in, for instance, pesaran2006estimation. It also ensures that the averages are informative about the factors. Finally, Assumption 6 adapts the parameterization necessary for our statistics in moderately integrated environments, because our asymptotic results will be derived under the local alternative setting.
Notice that the expressions for the feasible versions of the statistics in ((ref))-((ref)) are obtained by simply replacing $\widetilde{u}_{2,t+1}$ with $\widehat{u}_{2,t+1}$. One can reconcile feasible and infeasible statistics using the following asymptotic equivalence decompositions by using
for $j=1,\ldots, 4$, where explicit representations of $\xi_{N,T,j}$ come from margaritella2024new:
It is clear from the expansions of the feasible statistics that $\widetilde{u}_{2,t+1}-\widehat{u}_{2,t+1}$ is at the core of the upcoming analysis. It can be shown that this difference admits the following representation:
where each of the components brings a distinct contribution to the difference between infeasible and feasible forecast. For instance, $III$ takes into account the usual factor estimation error, which stems from the $r$ factors, while $II$ additionally involves an error from estimating $\+\alpha$ under the observed $\*f_t$. The component $I$ is the most subtle one. In particular, here $\widetilde{\+\delta}_t^0=[ \widetilde{\+\alpha}_t',\*0_{1\times (m-r)}, \widetilde{\+\theta}_t' \*]'\in \mathbb{R}^m$, which is $\widetilde{\+\delta}_t$ appended with extra $m-r$ zeros in place of the redundant averages. With $\widehat{\+\delta}_t\in \mathbb{R}^m$ being the LS estimator of the model parameters with $\widehat{\*z}_t$, we have $\*Q_N=\mathrm{diag}(\overline{\*H}\*D_N,\*I_q)$ transforming $\widehat{\*f}_t$ into $\widehat{\*f}^0_t$ inside of ((ref)). In effect, $I$ tracks the asymptotic effect of the factor overestimation error via the $m-r$ redundant averages. Lemma 1 below formalises its behavior uniformly over $t$.
There are significant differences between our Lemma 1 and its counterpart in stauskas2022tests. Firstly, we undertake our analysis under non-stationary predictors $\*z_t$. The subtle detail is that the feasible predictors are $\widehat{\*z}_t$ with the excess averages driven by stationary idiosyncratics. Therefore, we effectively have a mixed integration order of the regressors, which is reflected by the normalization matrix $\*D_T$. Strikingly, the appropriately normalised $\* Q_N^{-1}\boldsymbol{\widehat{\delta}}_t - \boldsymbol{\widetilde{\delta}}_t^0$ is asymptotically negligible even if $m>r$. The intuition behind these results is as follows. Let $\tau=0$, such that $\*D_T=T^{1/4}\*I_{q+m}$, to demonstrate that this difference does not originate solely from the persistence of the regressors. Then, we can show that
where
and $\left\|\*v_t \right\|= O_p(1)$ which is the same martingale difference process as in stauskas2022tests. The representation in ((ref)) immediately implies that $\left\|\sqrt{T}(\* Q_N^{-1}\boldsymbol{\widehat{\delta}}_t - \boldsymbol{\widetilde{\delta}}_t^0)\right\|=O_p(1)$ and so we are back to stauskas2022tests. The reason why we only need to scale up by $T^{1/4}$ and not $\sqrt{T}$ is because the former characterises the local power of the tests in pitarakis2023direct and pitarakis2025novel. In particular, the asymptotic distribution of his infeasible statistics is generated by $\{u_{t+1}^2-\sigma^2 \}$, whereas the components involving (functions of) $\{\*z_tu_{t+1} \}$ - which are the distribution generators in clark2001tests - are negligible. In other words, “the weight” is placed on the former martingale difference instead of the latter. Naturally, the same happens with the components of the feasible statistics, and so the influence of $\*v_t$ becomes negligible, as well. In general, Lemma 1 suggests that the redundant $m-r$ cross-section averages should not interfere with the asymptotic normality of the feasible statistics. In practice, it means that practitioners can use all the available blocks in FRED-MD or EA-MD-QD datasets and stay agnostic about the true number of factors, as long as $m\geq r$. This stands in sharp contrast to the CAs literature, where the redundant $m-r$ averages result in an asymptotic bias, unless $TN^{-1}\to 0$, which is a substantial restriction (see e.g. karabiyik2017role; or de2024cross). To the best of our knowledge, this is the first result in this literature when the redundant $m-r$ CAs stay completely harmless even if $TN^{-1}\to c>0$. \\ In connection to Lemma 1, we also use the facts that uniformly in $t$,
to ultimately demonstrate that uniformly in $t$ and a constant $C>0$
since $\left\|T^{-\tau/2}\widehat{\*z}_t^0 \right\|$ is uniformly bounded by the results in Lemma 3.1 of magdalinos2009limit.\footnote{Expanding $\*z_t$, it is $\sup_{k_0\leq t \leq T-1}\mathbb{E}\left(\left\|\frac{1}{T^{\tau/2}}\sum_{s=1}^t\*R_{z,T}^{t-s}\*u_{z,s} \right\|^2\right)=O(1)$, but $\left\|T^{-\tau/2}\widehat{\*z}_t^0 \right\|\leq \left\| T^{-\tau/2} \*z_t\right\|+\underbrace{\left\|T^{-\tau/2} \overline{\*e}_t^0 \right\|}_{o_p(1)}$ as $\overline{\*e}_t^0$ is stationary.} Note that ((ref)) holds for all $\tau \in (0,1)$, and we can show that no further restrictions are needed when we sum its squares, as required by the expansions of our statistics. These arguments can be used in the expansions ((ref)) - ((ref)) which are formalised by Lemma 2.
Notice that a direct consequence of the last part of Lemma 2 is that the variance estimators of the test statistics ((ref))-((ref)) are all consistent. In fact, it follows from the workings of pitarakis2023direct, pitarakis2025novel and margaritella2024new that $\left|\widetilde{\omega}_j^2 - \widehat{\omega}_j^2\right| = \theta_j\left|\widetilde{\phi}^2 - \widehat{\phi}^2 \right| = o_p(1)$ for each $j$, where $\theta_j$ is a function of tuning parameters specific to the particular statistic. As a result, it follows from the Continuous Mapping Theorem that the respective expressions $\left|\widetilde{\omega}_j/\widehat{\omega}_j-1 \right|$ are all negligible. By using Lemma 2, we can show that $s_{\widehat{f},1}$ - $s_{\widehat{f},4}$ admit the asymptotic representations reported below in Theorem 1:
Since the decompositions in expressions ((ref)) - ((ref)) are functions of the quantities in Lemma (ref), the results in Theorem (ref) follow readily. $\blacksquare$ \\ Three comments are in order. Firstly, the asymptotic theory of our four test statistics follows directly from pitarakis2023direct and pitarakis2025novel together with their major features, such as the ability to handle non-stationary predictors. The asymptotic representation in Theorem 1 holds without any restrictions on the relative expansion rate of $N$ and $T$ so long as they both diverge to infinity, thus making our testing procedure applicable to many large-dimensional macroeconomic and financial datasets. Secondly, note from ((ref)) that the rate of convergence in Lemma 2 is $O_p(T^{(\tau-1)/2})$. This indicates that the contribution of $m-r$ redundant averages may still decay slowly. Under $\mathbb{H}_1:\+\alpha^0\neq \*0_{r\times1}$, power can be reduced due to an exacarbated factor estimation error if $\tau\approx 1$. The fact that the error is still negligible is natural, as a mildly integrated process is weaker than (local-to-) unit root process ($\tau=1$). Lastly, the value of our statistics in empirical settings can be boosted even more due to their robustness to structural breaks in the factor loadings. We formulate this in Proposition 1.\\
Proposition 1. Let the loadings break at some time points $D_1=\lfloor \phi_1 T \rfloor,\ldots, D_d=\lfloor \phi_d T \rfloor$ with $\phi_1<\phi_2<\ldots<\phi_d \in (0,1)$ for a finite $d$, such that
Let $\*Q_i= [\boldsymbol{\Lambda}_{1,i}', \boldsymbol{\Lambda}_{2,i}',\ldots, \+\Lambda_{d+1,i}']'\in \mathbb{R}^{(d+1) r\times m}$. Under Assumptions 1-6 together with $\mathrm{rk}(\overline{\*Q})=(d+1)\times r$, Theorem 1 continues to hold.} \\
Proposition 1 further increases applicability of the tests, because it can be combined with the overspecification and nonstationarity of factors under Theorem 1. It states that under $d$ breaks in the factor loadings, the starting model can be reformulated into a model with $(d+1)\times r$ factors. Letting $d=1$ (a single break), we obtain a simple expression: $\boldsymbol{\Lambda}_{i,t}'\mathbf{f}_t = \mathbb{I}(t<D)\boldsymbol{\Lambda}_{1,i}'\mathbf{f}_t + \mathbb{I}(t\geq D)\boldsymbol{\Lambda}_{2,i}'\mathbf{f}_t$, which was also examined in stauskas2022tests, where we need $\mathrm{rk}(\overline{\*Q})=2r$. However, in their case, the condition of $D\in (1, k_0]$ is necessary. In other words, the effect of the structural break must be subsumed in the initial estimation of the parameters before going into the out-of-sample analysis. The reason is that the asymptotic distributions of statistics in clark2001tests are highly non-standard. If $D\in (k_0, T-1]$, there is a break in the asymptotic variance of the statistics. While it is typically possible to simulate the non-standard distribution, the break in the variance makes the simulation infeasible as the break date is unknown. In our Proposition 2, the break location does not matter. Intuition follows from the discussion below ((ref)). Our tests have the property that normality is generated by $\{u^2_{t+1}-\sigma^2\}$ independently of the predictors. Therefore, the factors and their loadings (even with breaks) do not interfere with asymptotic normality. Clearly, the predictors matter under $\mathbb{H}_1:\+\alpha\neq \*0_{r\times 1}$. However, we demonstrate that the factor estimation error is asymptotically negligible (Lemma 2) and does not affect the power function. It remains negligible even under breaks because we re-write the initial model with $r$ factors as the model with an inflated number of factors. \\ Our last result takes into account the uncertainty around the tuning parameters in statistics $2$-$4$ from pitarakis2025novel and margaritella2024new. Since different choices may lead to varying results in statistical power, Proposition 2 offers feasible power enhancements. \\
Proposition 2. Let $s_{\widehat{f},j}^{P}=s_{\widehat{f},j}+\widehat{\triangle}_j$ for $j=2,3,4$ be a power-enhanced statistic with the enhancement term $\widehat{\triangle}_j$. Then, under Assumptions 1-6 as $(N,T)\to \infty$
} Here, the power-enhancing term $\triangle_j$ is based on $\breve{u}_{2,t+1}^2= \widetilde{u}_{2,t+1}^2-(\widetilde{u}_{1,t+1}-\widetilde{u}_{2,t+1})^2$. Clearly, we only have the feasible $\widehat{\triangle}_j$, which obeys $\widehat{\triangle}_j=\triangle_j+o_p(1)$ according to Proposition 2. Due to the asymptotic equivalence, we will focus on the power-enhanced statistics in our Monte Carlo experiments and in the empirical application. \\
Remark 1 ($h$ steps). We can consider multi-step forecasts. We can update Assumption 1 to $\mathbb{E}(u_{t+h}|\mathcal{F}_t)=0$ ($h>1$) similarly to bai2006confidence. This implies that $\{u_{t+h}\}$ follows an $MA(h-1)$ process, which is sufficient for $\frac{1}{\sqrt{T}}\sum_{s=1}^{t-h}(u_{s+h}^2-\sigma^2)$ to generate Brownian motion, similarly to pitarakis2023direct. In our proofs, many terms that should vanish can be written as $\frac{1}{T}\sum_{t=1}^{T-h}g_tu_{t+h}$. Here, $g_t$ depends on $T^{(1+\tau)/2}(\widetilde{\+\delta}_t-\+\delta), T^{(1+\tau)/2}(\widetilde{\+\alpha}_t-\+\alpha)$ and similar recursive LS deviations. As $MA(h-1)$ induces a very limited serial correlation of forecast errors, we can verify the convergence rates of the latter terms (see Section 2.2 and 5.1 in the Supplementary material). Effectively, $g_t=g_t(\ldots,u_{t-1},u_t)$, and $\mathbb{E}(g_tu_{t+h})=0$ by the law of iterated expectations. It follows that $\mathbb{E}\left(\left|\frac{1}{T}\sum_{t=1}^{T-h}g_tu_{t+h}\right|^2 \right)=\frac{1}{T^2}\sum_{t=1}^{T-h}\sum_{s=1}^{T-h}\mathbb{E}(g_{t+h}g_{s+h})\mathbb{C}ov(u_{t+h},u_{t s+h})=O(T^{-1})$ due to $MA(h-1)$ structure and a truncated sum. In Section 4.5, we examine our statistics in simulations for $h\in \{1,4,12,24 \}$ and find $s_{\widehat{f},1}$ and $s_{\widehat{f},4}$ to show the best power across the horizons. Clearly, ((ref)) is replaced by HAC estimator as recommended in pitarakis2025novel to take serial correlation into account when estimating variances of our statistics. \\
Remark 2 (Number of factors). We can set $m=r$ by estimating the number of factors similarly to the PC literature (see ahn2013eigenvalue), which is relevant to avoid possible small sample distortions of size and/or power when $m>r$ (we will document this for $s_{\widehat{f},1}$ and $s_{\widehat{f},4}$). This can be achieved by using the Information Criterion (IC) by de2024cross (see their Proposition 3), which consistently selects an optimal set of CAs. Although their result holds for $\tau=0$, ditzen2025selection extended it to $\tau>0$. If $\tau$ is not too large (e.g $\tau \leq 0.5$), the IC can be trusted without any modifications. If $\tau \approx 1$, a consistent selection is possible when both $N,T$ are very large. However, if data of such dimensions are available, our Theorem 1 can be trusted to eliminate the effect of $m-r$ redundant CAs as $(N,T)\to \infty$.
In this section, we introduce the data-generating process (DGP) - which is very general and encompasses a plethora of empirically relevant scenarios - of our Monte Carlo simulations. Then present a large number of simulation results in Section (ref).
For the purpose of Monte Carlo simulation, we propose the same DGPs of stauskas2022tests and set, in the factor-augmented regression, $\* w_t = y_t$, $\theta = 0.5$, $\alpha = 0$ in the size experiment and $\alpha \in \{0.1, 0.2, 0.3, 0.4, 0.5, 0.6\}$ in the (graphical) power experiments. Various cross-sectional and time dimensions are investigated depending on the simulation considered, while the choice of $\pi_0$ is 0.25 so that $k_0 = \lfloor\pi_0 T \rfloor$ becomes a quarter of the total sample, and the forecast horizon is $h=1$. We choose the number of true factors to be $r=1$ and and the number of blocks to be $m = 3$ (unless otherwise stated). Consider the following general setup:
We examine a large number of modelling specifications of $\*\Lambda_{i,t}$, $\boldsymbol{\epsilon_{i,t}}$, $u_t$, $\delta$, $\rho_i$, $\beta$ and $\tau$. For conciseness, we summarise our proposed DGPs - labelled (1) to (9) - in Table (ref). To start with Equation ((ref)), we explore different distributional assumptions for $u_t$, which is drawn from $\mathcal{N}(0, 1)$, $t(10)$ or $\mathcal{N}(0, (1-\phi) + \phi u_{t-1}^2)$ for $\phi = \frac{1}{2}$ across the various DGPs. These three modelling choices allow us to investigate the effect differences of a standard normal distribution, fat tails and ARCH effects of first order on our tests. Further, Equation ((ref)) describes the DGP of the factor model. Notice here that the loadings specification encompasses both the time-invariant, nonzero mean
case as well as the time-varying cases in DGP (8), which incorporates a structural shift in the loadings mean from 1 to 2 at time $T/2$ as in breakschendoladogonzalo. Notice that the break is in the out-of-sample period in order to test the location invariance predicted by our theory. On the contrary, DGP (9) has time-invariant loadings with mean zero, which implies that $\overline{\+\Lambda}$ is asymptotically zero and CAs are inconsistent. Finally, we choose in Equation ((ref)) a common autoregressive specification\footnote{Notice that if $\delta = 0$ then the factors can be considered as a 'static' sequence of i.i.d. shocks. Hence, even if $\alpha \neq 0$, the probability of rejecting the null hypothesis when it is false (test power) should converge to the desired significant level (as for the test size) because the factors are expected to be uninformative.} for the time dynamics of the factors (see bai2006confidence, gonccalves2017tests, and StockWatson2002, for instance). In a similar manner, in Equation ((ref)), we allow the panel idiosyncratic components $\* e_{i,t}$ to be weakly dependent in $i$ and $t$ (similarly to bai2002determining, BanerjeeMarcellinoMasten2008, boivin2006more, and breakschendoladogonzalo), meaning that they exhibit serial correlation via the coefficient $\rho_i$ as well as cross-sectional correlation of spatial type by means of the $\beta$ coefficient in Equation $(\ref{eq: networ_dep})$. Note that the shocks $\boldsymbol{\epsilon}_{i,t}$ in this equation are drawn from $\mathcal{N}(\boldsymbol{0}_{m \times 1},\boldsymbol{\Sigma}_{\epsilon \epsilon,i})$, $\mathcal{N}(\boldsymbol{0}_{m \times 1},\boldsymbol{\Sigma}_{\epsilon \epsilon,i,t})$ or $t(10)$, meaning that the first distribution acts as a benchmark, whereas the second and third distributions introduce time-varying volatility and fat tails, respectively. Here, both covariance matrices $\boldsymbol{\Sigma}_{\epsilon \epsilon,i}$ and $\boldsymbol{\Sigma}_{\epsilon \epsilon,i,t}$ are simply diagonal and their non-zero elements are generated from $\mathcal{U}(0.5, 1.5)$ so to accommodate for cross-sectional heteroskedasticity, as permitted under CCE. Lastly, remember that we are working under the theory of moderately integrated systems \`a la magdalinos2009limit, whereby we assume that the parameters $\theta$ and $\delta$ can be written as $\theta = 1-\theta_0/T^\tau$ with $\left\lVert\theta_0\right\rVert < \infty$ and $\delta = 1-\delta_0/T^\tau$ with $\left\lVert\delta_0\right\rVert < \infty$. Hence, we consider $\tau \in \{0, 0.2, 0.4, 0.6, 0.8\}$ in the next section. For $j=2,3,4$, we consider power-enhanced statistics $s_{\widehat{f},j}^P$.
In this section, we report size and local power experiments run on 1,000 MC replications across different DGPs and parameter settings. The significance level is set to 5% and the critical values are obtained from the usual standard normal table. In detail, we present a comparison of test size and local power of the four statistics $s_{\widehat{f},1}$, $s_{\widehat{f},2}^{P}$, $s_{\widehat{f},3}^{P}$ and $s_{\widehat{f},4}^{P}$ in Table (ref) and Table (ref), respectively. In both analyses, we consider a setup under DGP (2) and two levels of persistence, i.e. $\tau = 0.2$ (mild) and $\tau = 0.6$ (moderate). DGP (2) assumes moderate dependence in the idiosyncratics and so it could be regarded as an unsophisticated but realistic representation of economic data under a factor model structure. Hence, it should serve as a benchmark for all other MC settings. Throughout the paper, we set the tuning parameters to $\mu_0 = 0.45$ for $s_{\widehat{f},1}$, $\nu_0 = 0.8$, $\lambda_1^0 = 1$ and $\lambda_2^0 = 0.65$ for $s_{\widehat{f},2}^{P}$ and $s_{\widehat{f},3}^{P}$, $\lambda_1^0 = 0.65$ and $\lambda_2^0 = 1$ for $s_{\widehat{f},4}^{P}$ as in margaritella2024new.
Starting from test size, Table (ref) reports the results in an environment with low persistence on the left-hand side and high persistence on the right-hand side. Observe that the test size converges to the significance level as $T \to \infty$ but it is seems to be insensitive to different levels of $N$. Further, it is evident that size does not deteriorate as $\tau$ grows, as anticipated by the theory of pitarakis2023direct. With respect to the performance of the individual statistics, it is clear that $s_{\widehat{f},2}^{P}$ and $s_{\widehat{f},3}^{P}$ deliver the best test size, while both $s_{\widehat{f},1}$ and $s_{\widehat{f},4}^{P}$ suffer from oversize for smaller values of $T$. However, additional results in the Supplementary material indicate that this oversize experienced by the latter pair is due to our overspecification of the number of blocks relative to the true number of factors, i.e. $m=3$ while $r=1$. By ensuring that $m=r$, it is therefore possible to materially reduce the issue of oversize while maintaining excellent local power. This can be achieved by using the IC by de2024cross - which remains valid under non-stationary factors, as demonstrated by ditzen2025selection - unless persistence is very high, in which case very large $N$ and $T$ combinations are required. In essence, we conclude that reliable test size is obtained with larger $T$ at all levels of persistence and for all statistics.
Next, we turn our attention to the local alternative, where we set $\alpha=0.5$. Given our theoretical results, we should point out that, loosely speaking, our tests should not discern between estimated and observed factors for sufficiently large block sizes, meaning that the theoretical framework of pitarakis2023direct,pitarakis2025novel is virtually restored. This is due to the fact that factor uncertainty inherent in CAs disappears with $N \to \infty$ in each block. The results are reported in Table (ref). Similarly to test size, we see a distinct difference in performance between the four statistics under both persistence environments, with $s_{\widehat{f},1}$ and $s_{\widehat{f},4}^{P}$ exhibiting significantly better local power relative to $s_{\widehat{f},2}^{P}$ and $s_{\widehat{f},3}^{P}$. Such constrast, however, becomes materially unimportant with larger values of $\alpha$. Contrary to our simulations on test size, whereas, we notice that higher persistence in the predictors lowers the local power of all statistics. This is perfectly consistent with our theory in Section (ref) where the convergence rates under the local alternative depend on $\tau \in [0,1)$: as the system of predictors becomes increasingly persistent, the tests devised by Pitarakis - and, by natural extension, our tests under CA-estimated factors - will be diminishingly able to detect enhanced forecasting power in the augmented regression relative to its nested baseline. Overall, we suggest using $s_{\widehat{f},2}^{P}$ and $s_{\widehat{f},3}^{P}$ for better test size and $s_{\widehat{f},1}$ and $s_{\widehat{f},4}^{P}$ for enhanced local power.
To demonstrate the robustness of our tests, we provide in Figure (ref) a comprehensive collection of power curves across DGPs (1) to (9) for all statistics, with the usual values $\alpha \in \{0.1, 0.2, 0.3, 0.4, 0.5,\allowbreak 0.6\}$ and $\tau = 0.5$. All statistics exhibit satisfactory power across the various DGPs but seem to suffer moderately in terms of both size and local power under ARCH effects and regression heavy tail settings. The former can be improved if ((ref)) is substituted by a HAC estimator since $\{u_{t+1}^2-\sigma^2 \}$ is serially correlated under conditional heteroskedasticity.
In general, it is unsurprising that the nondegenerate DGPs, namely (2) to (8), all exhibit improving local power as $\alpha$ moves away from zero. Note that the performance of DGP (8) with a break in loadings in the out-of-sample portion especially stands out. We should remark that the results are not excellent in the case of DGPs (6) and (7). The former performs significantly worse across the group in terms of size and power, whereas the latter exhibits excellent test size but its local power seems to match that of DGP (6) for higher levels of $\alpha$. In fact, these statistics are constructed based on averages of the (squared) predictive regression residuals so it is natural to think that their estimation is affected by ARCH effects and fat tails. Nonetheless, our findings seem to indicate that these issues are only of finite sample nature and vanish asymptotically with increasing time dimension. With respect to the first DGP, whereas, the values are flat around the significant level and should be interpreted as test size. The reason is that the panel data has no factor structure nor any kind of dependence, so that the CAs-estimated factors should not improve upon the baseline specification. Finally, DGP (9) models the rank violation so that the CAs are inconsistent for the factor space. The statistical inference is not valid in this case and, indeed, the local power is visibly impaired.\\ Finally, we apply another battery of simulation scenarios and report them in the Supplement in the interest of space. The interested reader should be informed in advance that, in line with previous findings, a similar analysis almost invariably applies to all four statistics across these additional settings, with excellent local power across the spectrum. The experimental evidence includes: 1) size experiments with a known number of factors ($m=r$), 2) the effect of different redundant blocks, 3) varying levels of persistence (also in the presence of serially correlated idiosyncratics), 4) violations of the MDS assumption of $u_{t}$. Overall, local power deteriorates in all of these settings in small samples, with the notable exception of 1) where the impact almost does not exist (but the size of $s_{\widehat{f},1}$ and $s_{\widehat{f},4}^{P}$ is improves significantly, as discussed above). To summarise, our CAs-based tests perform quite well under a plethora of empirically relevant economic scenarios. More importantly, our MC simulations seem to support two important theoretical contributions of this work to the factor model literature. First, our tests seem to be invariant to uncertainty in the number of factors and breaks in the loadings. Second, they can accommodate for moderately integrated systems of variables with acceptable loss of local power.
Following Remark 1, we explore the evolution of the power for $h=1,4,12,24$. We consider HAC standard errors for all statistics under the empirically appealing DGP(2) and report the results in Figure 2. To our surprise, the statistic $s_{\widehat{f},1,HAC}$ - proposed by pitarakis2023direct to handle multi-step forecasts - suffers from a loss in local power, while the remaining statistics seem to be unaffected. Despite this , $s_{\widehat{f},1,HAC}$ still outperforms $s_{\widehat{f},2,HAC}^P$ and $s_{\widehat{f},3,HAC}^P$ yet underperforms relative to $s_{\widehat{f},4,HAC}^P$. This is consistent to the previous results in our simulations. Indeed, the first and last statistics exhibit the best local power while the second and third statistics offer better statistical size, thus again corroborating the discussion on the trade-off between size and power.
In this section, we provide empirical evidence on the predictive power of CAs-estimated factors in macroeconomic forecasting. For this purpose, we use the novel, publicly available\footnote{Freely available at the page: \url{https://zenodo.org/records/14938961}} EA-MD-QD dataset of EAdataset. The dataset is comprised of quarterly and monthly macroeconomic time series for both the Euro Area (EA) as a whole and its ten primary member countries, and it is updated on a monthly basis and constantly revised. To the best of our knowledge, this is the first article to use this dataset in the context of forecasting and out-of-sample (factor-based) model evaluation. In detail, the EA dataset consists of 118 time series - which include 71 and 47 variables collected at quarterly and monthly frequencies, respectively - that have been recorded starting from 2000:1 to 2024:4. Similarly to the FRED-MD dataset mccracken2016fred, the variables are classified into 11 blocks: (1) National Accounts, (2) Labor Market Indicators, (3) Credit aggregates, (4) Labor Costs, (5) Exchange Rates, (6) Interest Rates (7) Industrial Production and Turnover, (8) Prices, (9) Confidence Indicators, (10) Monetary Aggregates and (11) Others. The country-level datasets have a comparable number of time series and the same classification groups, with the exception of Monetary Aggregates since it only recorded at EA level. We will apply our tests to the EA as a whole as well as a collection of individual countries, specifically Germany, France, Italy and Spain.
We treat the data according to the thorough suggestions of EAdataset. First, we aggregate the monthly data at the quarterly level to address the mixed-frequency nature of the variables. Next, while our testing framework is able to accomodate for different degrees of persistence in the time series, it is necessary to address the issue of missing values and outliers beforehand. To do so, we induce stationarity in the data according to the transformations recommended for each variable in Appendix A of EAdataset. Subsequently, we impute the missing values using the EM algorithm of StockWatson2002 and follow standard screening for outliers as detailed by mccracken2016fred. Completion of these steps results in a first dataset with fully treated stationary data. By reversing all first-differencing transformations, we also recover a second dataset with fully treated persistent data.
In the forecasting exercises, we choose a collection of target variables in our regressions so to display the predictive power and usefulness of the estimated factors. The selected target variables include: real GDP (GDP), industrial production index (IPMN), harmonised index of consumer prices excluding energy and food (HICPNEF), total unemployment (UNETOT), 3-month interest rates (IRT3M), 6-month interest rates (IRT6M), long-term interest rates (\textit{LTIRT}), wages and salaries (\textit{WS}), general government total financial liabilities (\textit{GGFB}), non-financial corporations financial liabilities (\textit{NFCLB}), households total financial liabilities (\textit{HHLB}), real exchange rate against 42 industrial countries (\textit{REER42}), residential property prices (\textit{HPRC}), and share prices (\textit{SHIX}). We then extract the factors from stationary and persistent datasets (after removing the target variable to be forecasted) and use them to augment an AR(1) benchmark. Note that richer specifications can be attained by including, among other possibilities, other exogenous regressors or extra lags in the predictors. While doing so can improve the predictive accuracy, we do not explore it in what follows. In addition, in real life, a forecaster may decide to choose subsets of the blocks so to achieve a more parsimonious forecasting model. For instance, gonccalves2017tests investigates the predictive content of factors for the equity premium associated with the S&P 500 Composite Index by comparing all permutations of models that include a single lag of at least one of the eight estimated macroeconomic factors. However, we are not excessively concerned with this since our tests are invariant under the use of redundant blocks after all.
To guarantee the validity of some of the assumptions in Section ((ref)), we run some diagnostics on the estimated factor model. In the interest of space, we report the results when the EA real GDP is the target variable in the forecasting equation, but similar conclusions for all other variables can be drawn based on additional results in the last section of the Supplement. The first assumption that we examine is that the number of blocks is at least as large as the number of true factors, i.e. $m \geq r$. To estimate the number of factors, we use the eigenvalue-ratio (ER) criterion by ahn2013eigenvalue with $r_{max} = 8$ and find one factor in the stationary data, which is in line with the findings of EAdataset. Following OutputGap, we also borrow this result in the context of the persistent data\footnote{To our knowledge, PilarPoncela2006 is the only paper in the literature that proposes a test for the number of non-stationary factors, although in a dynamic factor model framework.} since, as argued by OnatskiWang2021, the number of factors in levels should coincide with the number of factors in first differences provided that no spurious effects are estimated. EAdataset also implement the criterion of bai2002determining, which select seven factors. Given that ahn2013eigenvalue is often parsimonious in small sample while the information criterion of bai2002determining tends to overselect the number of factors, a researcher should probably accept that the true number of factors in this dataset lies between 1 and 7. As discussed next, this number may also change as the result of structural breaks.
Another assumption that we make is that the factor loadings are constant over time. Since the inception of the euro currency and the start of the dataset, the economies of the European Union (EU) members have experienced systemic shocks such as the Great Financial Crisis, a sovereign debt crisis and the Covid-19 pandemic. Thus, we raise the question about the possibility of structural breaks in the factor model. In fact, the presence of $d$ breaks in the sample requires us to reformulate Equation ((ref)) into a model with $(d+1) \times r$ factors, meaning that the number of CAs now needs to satisfy $m \geq (d+1) \times r$ unequivocably. In order to address our concerns about large shocks, we use the sup-LM statistic-based test of breakschendoladogonzalo to detect big structural breaks in these loadings at unknown dates. In particular, there is evidence against the no break null hypothesis at the timestamps 2020:2 and 2021:4 at the 5% significance level ($\widehat{d}=2$). Unsurprisingly, these dates exactly coincide with the beginning and the end of the national lockdowns imposed by EU countries, which had systemic and material repercussions on the entire Euro area. With this finding in mind, we then confront a consequential concern and inspect if $r$ changes due to these structural breaks in the factor loadings. Indeed, this concern is not new in the literature as it was previously raised, for instance, by breitung2011testing and BaiNg2007 in the context of US macroeconomic data. Further, OutputGap argues that the effect of the Covid shock is pervasive in the time series of most real macroeconomic variables and should therefore be treated as an additional latent common factor. In fact, the Covid shock induced a large shift not only in the levels MarozStockWatson, Ng2021 but also in the volatility CarrClarkMarcMert2022,LenzaPrim2021 of macroeconomic variables. Hence, we apply once more the ER criterion to the subsamples between the structural breaks: 1) 2000:2-2020:1, 2) 2020:2-2021:4, and 3) 2021:4 to the end of the sample. This check follows Section 4.4 of breakschendoladogonzalo to address the possibility of factor heteroskedasticity. We find that $r = 1$ across all three subsamples. It should be noted that even if $r = 3$ over the entire sample, the number of CAs ($m=11$) is still not smaller than $(\widehat{d}+1) \times \widehat{r} = (2+1)\times 3=9$ so we conclude that we have sufficient estimated factors. Having said this, should any of the timestamps 2020:2 and 2021:4 correspond to a break in the factor covariance matrix, the sup-LM statistic will at worst provide an overestimate of the number of true breaks in the factor loadings, meaning that our tests are still robust.
Before turning to our equal predictive ability and encompassing tests, one more diagnostic step remains. Indeed, another assumption in our setup is that the idiosyncratic components in the Equation ((ref)) are at most weakly dependent cross-sectionally and in time with stationary and persistent factors but serially uncorrelated with non-stationary predictors. To illustrate that the assumption of weak cross-sectional correlation in the idiosyncratics holds, we apply the sparse estimation methodology via adaptive thresholding of CaiLiu2011 to the covariance matrix of the panel residuals in the stationary and persistent cases. For the former case, we find that the estimated matrix is almost exactly diagonal with average (resp. average absolute) value of the cross-correlation of 0.002 (resp. 0.007), which suggests that the factors estimated by CAs are able to explain most of the variation in the co-movements of the cross-sectional units. For the latter case, whereas, the estimated matrix is non-diagonal but sparse with average (resp. average absolute) value of the cross-correlation of 0.005 (resp. 0.036). We regard this level of cross-sectional dependence as very weak overall. On the other hand, serial correlation in the idiosyncratics seems to be more assiduous according to standard Ljung-Box Q-testing. In the stationary dataset, we find average (resp. average absolute) serial correlations of order 1 equal to -0.002 (resp. 0.184). Conversely, the corresponding value is 0.558 (resp. 0.559) in the persistent case. While our findings are not of concern in the context of stationary data, they raise apprehension when it comes to persistent data since we do not allow for serial correlation in association with non-stationary regressors in our theory. To alleviate some distress, we report additional Monte Carlo simulations as a benchmark in the Supplement where we explore the effect of strong serial correlation in DGP (3) with $\tau \in \{0, 0.2, 0.4, 0.6, 0.8 \}$ on size and local power. Our results seem to indicate for all four statistics that, in spite of our rates of convergence, test size is fundamentally unaffected while local power suffers a small but acceptable loss. As a result, we suspect that our convergence rates may possibly be conservative in this scenario. Since time and cross-sectional dependence in DGP (3) is stronger than the one observed in the persistent data, we trust that our equal forecasting ability and encompassing tests are still able to convey practical insight into the predictive power of the estimated non-stationary factors. With all these pre-testing analysis in mind, we finally dedicate our attention to testing predictive accuracy.
Upon estimation, we apply our four novel tests of equal predictive accuracy and encompassing to assess whether the predictive content of CAs improves upon the forecasting ability of the autoregressive baseline with the selected target variables. The results are reported in Table (ref).
First and foremost, it is paramount to highlight that the predictive content inherently embedded in our CAs-estimated factors is particularly pronounced with target variables that play a key role in the design of political, fiscal or monetary policies. These variables are economic growth, inflation (excluding energy and food), unemployment and wages. Since monitoring inflation is at the core of the European Central Bank's mandate and policymaking, it is very valuable to observe that the estimated factors have formidable forecasting power for consumer prices under all tests and types of persistence. Focusing on individual tests, we see that $s_{\widehat{f},1}$ exhibits statistical significance across the largest collection of target variables. While this statistic exhibits power, we know from our simulations that it tends to overreject the null. Therefore, the support from other statistics, especially $s_{\widehat{f},2}^{P}$ and $s_{\widehat{f},3}^{P}$, strengthens the evidence for factor predictive power. Unsurprisingly to our theory, our statistics detect better overall forecast accuracy using stationary data, which is still against the popular belief that stationary data incurs a material loss of information content compared to non-stationary data. In fact, this finding was already anticipated by the fact that the convergence rates of our statistics slow down with increasing persistence in the data. Further, while the ultimate objective of our tests is to uncover the predictive specifications where the estimated factors improve upon the (autoregressive) baseline, it is also important to identify those situations where the factors are not informative about the future and should therefore be excluded from our forecasting considerations, as is rightfully and conventionally the case for the changes in, and the level of, the long-term interest rate, the exchange rate and the stock market returns. In fact, these three variables are notoriously difficult, if not impossible, to predict because - assuming the efficient market hypothesis holds - the expectations of rational, informed traders must already be priced into financial markets.
Nevertheless, we acknowledge that our results may inadvertedly be affected by structural breaks in the coefficients of the factor-augmented regression. Although we have dedicated considerable attention to breaks in the factor loadings, there exist no statistical procedure to detect structural instability à la CorradiSwanson1014 under the CAs framework. In fact, these authors propose a test for the joint hypothesis of structural stability of both factor loadings and regression coefficients in the factor-augmented forecasting model in a PC setting. As previously argued, the period associated with the Covid-19 pandemic should be treated by including an additional factor in the specification, which will by construction capture a portion of the shift in the level of the target variables. Hence, overestimating the number of CAs makes us robust to: 1) the possibility of breaks in the factor loadings, and consequently 2) an increase in the number of factors resulting from a large shock, but possibly also 3) a large break in the factor-augmented regression coefficients. However, it is not unreasonable to think that some concerns about the last point have not been dispelled once and for all. At the same time, our empirical results would not be interpretable fairly under the tests of clark2001tests and mccracken2007asymptotics due to their strong foundational assumptions about factor properties. \\ We repeat the same forecasting exercise with data on Germany, France, Italy and Spain. The target variables are the same, with the exception of the 3-month and 6-month interest rates since they are not available at the country level. The results of our tests are reported in Table (ref).
Although the institutions of these countries are structurally distinct, the results will naturally be heterogeneous across countries. That is, the factors extracted from each country dataset are found to be informative about the future evolution of different variables. In sum, it seems from our results that the macroeconomic factors are still able to offer consistent forecasting insight into economic growth, wages and (possibly) consumer prices across the four countries, with the notable exceptions of Spain and Italy. With respect to Spain, the factors appear to show quite weak predictive ability while, for Italy, they anticipate the bulk of economic activity across most variables, excluding the stock market. As before, $s_{\widehat{f},1}$ detects the greatest number of specifications where our macroeconomic factors yield forecasting power. All things considered, we conclude that - although the CAs method may appear deceivingly simplistic at first glance - predictive regressions augmented with the estimated factors offer powerful predictive insight in the context of macroeconomic forecasting in both levels and differences, especially in the presence of possibly big structural breaks.
We present four novel tests of equal predictive accuracy and encompassing for factor-augmented regressions under a CAs setup by revisiting stauskas2022tests. While their results in the context of tests by clark2001tests bring important practical implications, merging CAs with a novel class of tests by pitarakis2023direct and pitarakis2025novel opens new possibilities. The overspecification of the number of factors has an asymptotically negligible influence, thereby eliminating the need to know the number of factors (as long as it is bounded from above by the number of averages). This goes against the usual trend in this literature according to which the excess CAs result in bias even asymptotically. In practice, the suggested upper bound is often substantial, as can be seen from our empirical application on the novel EA-MD-QD dataset. Further, observed predictors and factors can be persistent, and the location of structural breaks in the factor loadings need not be known. Ultimately, Monte Carlo simulations reveal an excellent small sample performance across many empirically relevant scenarios.