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.
97,159 characters · 14 sections · 139 citation commands
New Tests of Equal Forecast Accuracy for Factor-Augmented Regressions with Weaker Loadings
{0.7cm}
\etocdepthtag.toc{mtchapter} \etocsettagdepth{mtchapter}{subsection} \etocsettagdepth{mtappendix}{none}
Assessing the out-of-sample forecast performance of different models is fundamental for practitioners deciding which specification to employ. Central banks, investment banks, and government economic planning agencies routinely employ forecasting models to make informed decisions and implement effective policies. Therefore, evaluating the population predictive ability of various alternative models is for them essential.\footnote{We stress the out-of-sample and in-population because “in-sample" mean square error (MSE) based forecasts, and horse races thereof, are void of statistical grounds, see a.o., diebold2015comparing.} Competing forecasts may differ in setup (e.g., forecasting inflation via linear regression or exponential smoothing), but the focus is often on augmenting a linear model with additional predictors. In such cases, models are nested under the null hypothesis of equal forecast accuracy. Nesting comes naturally in the validation of economic theories. For instance, inflation ($\pi_t$) can be forecast using an AR(1) model, or alternatively, using an ARX(1) model as $\pi_t=\gamma_1+\gamma_2 \pi_{t-1}+\gamma_3\iota_{t-1}+\nu_t$, where interest rates ($\iota_t$) serve as exogenous input. Under the null hypothesis then, $\gamma_3=0$. Another example of this is forecast accuracy comparison against random walk models in the field of exchange rate considered by rossi2005testing (see pitarakis2025novel for an overview).
While natural, nowadays nested comparisons go hand in hand with a high-dimensional setting. With the availability of contemporary large datasets, researchers and practitioners face an abundance of potential predictors that may or may not improve forecasts beyond standard autoregressive specifications. By far the most popular approach to 'condense' the predictive information is employing the factor structure (see a survey in eickmeier2008successful). Here, a large number of potential predictors load on a small number of latent series (factors) that drive their co-movement.\footnote{To be precise, factor models are not a dimension-reduction technique per-se; it is only the reduced rank assumption on the common component, combined with white noise-like assumptions on the idiosyncratics, which have the desired reduction effect barigozzi2024dynamic. However, for macroeconomic applications, these are often reasonable assumptions.} Therefore, factor-augmented forecast is the hallmark example of a nested setup. It also constitutes a whole research industry with applications in macroeconomics and finance, where Principal Components (PC) is the predominant method of factor estimation which extracts eigenvectors corresponding to the largest eigenvalues of the sample covariance matrix (see bai2002determining; or bai2006confidence). Hence, it is not surprising to see two parallel strands of econometric literature. The first focuses on the robustness of the PC procedure itself, while the second develops tests to evaluate out-of-sample forecasts with estimated factors. An important example of the former is empirically well-documented and it concerns risks when the loadings through which the factors weight on the high-dimensional set of observables are weak (see examples in stock2002macroeconomic, or ludvigson2009macro). Technically, the eigenvalues diverge at a sub-linear rate (see e.g. uematsu2022inference; bai2023approximate). Practically, it limits the informativeness of a data set and may cloud forecast comparisons. Indeed, the key examples of the second strand are built only on the strong loadings assumption, and they are the seminal work of gonccalves2017tests and stauskas2022tests\footnote{This study, however, exploits the Common Correlated Effects (CCE) estimator by pesaran2006estimation. While this method has elegant properties when the data admit a specific structure (e.g. blocks), we focus on PC due to its versatility.}. They derive conditions under which (highly non-standard) asymptotic distributions of the tests of clark2001tests continue to hold when the competing model is augmented with estimated factors. In the current study, we merge both strands. Firstly, we re-establish standard normal inference by borrowing from the recent developments in the theory of nested environments by pitarakis2023direct,pitarakis2025novel. Next, we introduce the possibility of weaker loadings to quantify their effect on the forecast comparison problem.
Statistical comparison of nested models typically results in a non-standard inference. Particularly, under the null, the population errors of the two competing specifications are identical, leading to zero out-of-sample mean square error (MSE) differentials in the limit, as well as zero asymptotic variances. As a result, the test statistics become asymptotically ill-defined. On the contrary, these problems are bypassed in the non-nested case, and asymptotic normality is relatively easy to establish diebold2015comparing. Such challenges motivated early works by clark2001tests and mccracken2007asymptotics who introduced adaptive normalizations of the MSEs to recover well-defined asymptotic distributions of the test statistics. However, these distributions are highly non-standard and follow stochastic integrals of Brownian motion that depend on the relative growth rate of in-sample versus out-of-sample observations. Simulation-based approaches for estimating asymptotically valid critical values exist (see clark2012reality; or hansen2015equivalence), but their practical implementation is very challenging, which makes the results of gonccalves2017tests elegant, but impractical.
More recently, Jean-Yves Pitarakis (JYP) proposed a set of tests for predictive accuracy and encompassing in nested models, avoiding the non-standardness of the asymptotic distribution of the statistics under the null, as well as the variance degeneracy issue of existing procedures. pitarakis2025novel proposes a forecast accuracy test, still based on the MSE comparison across two nested models. However, in comparing MSEs, it uses partially overlapping out-of-sample segments to compute them, rather than the whole, same out-of-sample span. The intuition goes that as long as the fractions of out-of-sample squared forecast errors associated with the two competing specifications are different, the variance of a suitably normalized test statistic involving the MSEs cannot be degenerate. JYP then proves asymptotic normality of two types of tests under a set of general, nonrestrictive assumptions (more details in Section (ref)).
pitarakis2023direct instead, is an encompassing test. Meaning, it is based on the forecast encompassing principle hendry1982formulation for which if one forecast offers no additional value over another, then an optimal convex combination of the two cannot yield a lower squared error loss. Defining $\widetilde{u}_{1,t+1}, \widetilde{u}_{2,t+1}$ as the one-step-ahead forecast errors associated with two alternative forecasts, this principle boils down to testing the population moment restriction: $\mathbb{E}({u}_{1,t+1}({u}_{1,t+1}-{u}_{2,t+1}))=0$. However, suitably normalized sample statistics to test such restriction are also plagued by the same issues that plague the equal predictive ability tests: variance degeneracy and nonstandardness of asymptotics. To circumvent these, the simple and brilliant idea in pitarakis2023direct is that of considering the linear combination of two subsample means, in place of a unique sample mean, as additive sample counterparts of the population quantities. Simply enough, if $[1:k_0]$ is the set of in-sample observations and $[k_0+1:T]$ the out-of-sample's, this means that the sample counterpart of $\mathbb{E}[{u}_{1,t+1}{u}_{2,t+1}]$ is not going to be $(T-k_0)^{-1} \sum_{t=k_0}^{T-1} \widetilde{u}_{1,t+1}\widetilde{u}_{2,t+1}$ but $1/2 ( m_0^{-1}\sum_{t=k_0}^{m_0+k_0-1}\widetilde{u}_{1,t+1}\widetilde{u}_{2,t+1}+(T-k_0-m_0)^{-1} \sum_{t=m_0+k_0}^{T-1} \widetilde{u}_{1,t+1}\widetilde{u}_{2,t+1})$, for $m_0$ being a split point that should be different from $(T-k_0)/2$. The two subgroups of units will have slightly different variances, and that is the key to circumvent the variance degeneracy issue. As for the previous tests, JYP proves asymptotic normality of his proposed encompassing test statistics under a set of general, high-level assumptions.
JYP's tests, for their simplicity, robustness and generality, should be considered the new standard procedures to test for forecast accuracy and encompassing. Crucially, they have potential to make the out-of-sample evaluation of factor-augmented forecasts much more practical. However, as of now, there is no theory formally justifying their appealing properties in this setting as his high-level assumptions accommodate only observed, but not estimated factors. pitarakis2023direct only considers a factor-augmented data generating process (cf. DGP2) as a robustness check in the simulations. \footnote{Also, upon reviewing of his simulations, it appears that the factors are not re-estimated at every roll of the out-of-sample window, but are just estimated once, prior to the recursion, over the whole out-of-sample span. This neglects the very out-of-sample nature of the setting (pitarakis2023direct simulation scripts are available on the author's GitHub page at \url{https://github.com/jpitarakis/Multi-Step_Encompassing}).} Moreover, given the potential for a wide applicability, it is only natural to provide the justification by reflecting the risks that occur in practice. Indeed, various routes have been taken to accommodate the weakness of factor loadings, which include local factors (see freyaldenhoven2022factor) or sparsity which prescribes a sparse PC as an alternative estimator (see uematsu2022inference). Only recently bai2023approximate provided a comprehensive theory for the usual PC to determine the impact of the loading weakness. It is also more general than in the latter studies. Particularly, if assumptions of uematsu2022inference were applied to the usual PC setting, factor loadings would need to be stronger for the factors to be identified.
As the key contribution, we provide the theoretical framework that justifies the use of the whole set of JYP's tests of forecast accuracy and encompassing, when the competing specification is that of a factor augmented autoregression. The factors are estimated at every roll of the out-of-sample window (recursively), and the factor loadings are allowed to be strong, as well as homogeneously or heterogeneously weak. This leads to a technical contribution, where we re-work the weaker loadings theory in bai2023approximate to the recursive estimation setting, which can be of independent interest. We spell and discuss the assumptions required to show how in all these settings JYP's test statistics maintain Gaussian asymptotic distribution and appealing power properties, as long as the loadings are at most moderately weak. This extra informativeness condition is required to take into account the recursive estimation of factors.
Finally, few words on the notation we are going to use throughout. Firstly, $a$ is a scalar, $\*a$ is a vector and $\*A_t$ is a matrix with $t$ rows. For any generic matrix $\*A$, spectral norm is $\left\|\*A \right\|_{sp}=\sqrt{\lambda_{\mathrm{max}}(\*A'\*A)}$, while $\|\*A \|=\sqrt{\mathrm{tr}(\*A'\*A)}$ is the Frobenius norm with $\mathrm{tr}(.)$ being the trace operator. Vectorization of a matrix $\*A$ is denoted by $\mathrm{vec}(\*A)$, and $\lambda_{\min(\max)}(\*A)$ denote the smallest (largest) eigenvalue. Next, $\lfloor x \rfloor$ represents the integer part of $x$, and $M$ is a positive constant, while $\sup_{a\leq t \leq b}$ ($\inf_{a\leq t \leq b}$) is supremum (infimum). Moreover, $k_0$ represents the in-sample observations, while $T-k_0=n$ are the out-of-sample observations. Convergence in distribution and probability are given by $\to_d$ and $\to_p$, respectively, while weak convergence is given by $\Rightarrow$. Ultimately, $\widetilde{a}$ and $\widehat{a}$ are quantities estimated under the observed and estimated $\*f_t$.
For $t=1,\ldots,T,$ let us consider the following forecasting model
where $\*w_t\in \mathbb{R}^k$, $\*f_t\in \mathbb{R}^r$, which are stacked into $\*z_t\in\mathbb{R}^{k+r}$ with conformable parameter vectors: $\+\theta$, $\+\beta$ and $\+\delta$. We consider $\*w_t$ being the “known factors", this can contain both lags of $y_t$ as well as an intercept, seasonal dummies or time period dummies. Instead, $\*f_t$ is a vector of “unknown factors". In particular, we assume that there exist a panel of $N$ series (which excludes $y_t$) whose components $x_{i,t},\; i=1,\ldots,N,$ can be decomposed into two unobservable and mutually orthogonal components: a common component $\chi_{i,t}$ and an idiosyncratic component $e_{i,t}$. Respectively, they represent the comovements and the individual features of the series. For $\chi_{i,t}$, we assume it being low-rank, i.e., to be driven linearly by an $r$-dimensional vector of common static factors $\*f_t$, such that: $\chi_{i,t}=\+\lambda_i'\*f_t$, with $\+\lambda_i$ being, for every $i$, an $r$ dimensional vector of factor loadings. Thus, the decomposition takes the following form
or in matrix notation
where $\*F=(\*f_1,\ldots,\*f_T)'\in \mathbb{R}^{T\times r}$, $\+\Lambda=(\+\lambda_1,\ldots,\+\lambda_N)'\in\mathbb{R}^{N\times r}$ is the matrix of individual factor loadings and $\*E \in \mathbb{R}^{T\times N}$ is the matrix of idiosyncratic components.
We are interested in understanding whether the factor-augmented model ($\+\beta\neq \*0_r$, “unrestricted") is on average better in terms of out-of-sample forecast accuracy than a simple, possibly autoregressive model ($\+\beta = \*0_r$, “restricted"). For this purpose, we will split the sample $T$ into $T=k_0+n$, where $k_0$ and $n$ are in- and out-of-sample periods, respectively. Conveniently, we let $k_0=\lfloor T\pi_0 \rfloor$ for $\pi_0\in (0,1)$. Then, we produce recursive pseudo out-of-sample forecasts for the restricted and unrestricted models, for $t=k_0,\ldots,T-1$, and compare their errors of the form $\widetilde{u}_{1,t+1}=y_{t+1}-\widetilde{\+\theta}_t'\*w_t$, against $\widetilde{u}_{2,t+1}=y_{t+1}-\widetilde{\+\theta}'\*w_t-\widetilde{\+\beta}'{\*f_t}=y_{t+1}-\widetilde{\+\delta}_t'\*z_t$ for $t=k_0,\ldots,T-1$. Here and henceforth, the notation $\widetilde{a}$ $(\widehat{a})$ indicates an infeasible (feasible) estimator; as such, the least squares (LS) estimator
is indeed infeasible, given that the factors $\*f_t$ are unobserved. However, we are first going to proceed as if $\*f_t$ was given. Reason being that the expansions of the feasible test statistics reveal how the key component in this analysis is indeed the difference between the infeasible and feasible forecast error: $\widetilde{u}_{2,t+1}- \widehat{u}_{2,t+1}$, similarly to gonccalves2017tests and stauskas2022tests. Boundedness of functions of this quantity is indeed the key to the main theoretical results in this paper.
To begin with, we consider in total 3 tests of equal forecast accuracy and encompassing. The test for forecast encompassing is given by pitarakis2023direct and it has the following form:
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)$, with $\mu_0\neq 1/2$; $\widetilde{\omega}_1$ is the estimated standard deviation of the limiting distribution of the test statistic. Two further tests of forecast accuracy come from pitarakis2025novel:
and
Here, $l^0_j=\lfloor n\lambda_j^0 \rfloor$ for $j=1,2$ and $\lambda _j^0\in (0,1)$, controls the two portions of the out-of-sample period over which the forecast errors of both models are compared. One can have $l^0_1>l_2^0$ or vice-versa, which means that both portions are overlapping. However, an equality is ruled out in order to avoid the asymptotic degeneracy of variance. This test has a slight disadvantage when compared to the pitarakis2023direct, because of the overlapping evaluation over the effective sample size, some data, though in principal minimal, is lost in the MSEs comparison. Note also how ((ref)) is simply an average of ((ref)) over some chosen feasible set of $l_1$ for the fixed $l_2^0$ (or, fixed $\lambda_2^0$). The tuning parameter $\tau_0\in (0,1)$ helps to pick that set.\footnote{For a further discussion on how to choose the tuning parameters, we refer to pitarakis2025novel, Section 3.} Intuition behind the averaging is the following: if $l_2^0$ is fixed and $l_1$ changes, then the MSE of the restricted model accumulates. In effect, the uncertainty over the possible choices of $l_1$ is integrated.\\ There are several suggestions for different avenues of averaging in pitarakis2025novel. One that is significant both practically and theoretically is averaging over $l_2$ while $l_1^0$ is fixed, such that the MSE of the unrestricted model accumulates. Because in practice the feasible version of the unrestricted model will contain a factor estimation error, we can track if it interferes with an integration of the uncertainty around $l_2$. This results in a new statistic
This statistic is not present in pitarakis2025novel, therefore, as an additional contribution, we provide its full analysis. Lastly, $\widetilde{\omega}^2_j$ for $j=1,\ldots,4$ represent the variance estimators of the four statistics.\\ Algorithm (ref) below summarizes $g_{f,1}$ - $g_{f,4}$ and their respective variance estimators $\widetilde{\omega}^2_1,\ldots, \widetilde{\omega}_4^2$. It also provides recommendations for the tuning parameter values, which serve to boost statistical power of the tests. The recommendations stem from a theoretical and simulation-based investigation in pitarakis2025novel, while in our simulations we experiment with different values to explore the balance size and power.
Note that for $j=1,\ldots, 4$, we have $\omega_j^2=\phi^2\gamma^2_j$, where $\gamma_j^2$ is known. We only need to estimate $\phi^2$ with
where in practice we use the feasible estimator $\widehat{\phi}^2$. Similar sample variance estimators have been widely used in different exercises for factor-augmented regressions (see e.g. bai2006confidence, gonccalves2017tests, or yan2022factor)\\ For further exposition it is convenient to split ((ref)) - ((ref)) into two components, such that for $j=1,\ldots,4$ we have
where $g_{f,j,1}$ generates the distribution under the null, $g_{f,j,2}$ generates power under the alternative. Then, as $T\to \infty$, we shall have that, under the null, the $\widetilde{\omega}_{j}^{-1}$-scaled first component converges in distribution to a standard Gaussian density function. To see its key difference from the tests in clark2001tests, let $\sigma^2=\mathbb{V}ar(u_{t+1})$ (unconditional variance). Then it is the process $\{u_{t+1}^2-\sigma^2 \}$ that generates the distribution, instead of the predictors. For the local power analysis in $g_{f,j,2}$, we impose $\+\beta=\+\beta^0T^{-1/4}$ as the local alternative similarly to pitarakis2023direct, pitarakis2025novel. Thus, we require the following assumptions on the forecast error (A.1) and on the factors and predictors (A.2). They are slightly more primitive than in the original studies in order to make them compatible with and sufficient for the feasible setup explored later.
A.1 (Forecast Error)
A.2 (Factors and Predictors)
A.1 i) considers the filtration at $t-1$, hence this implicitly means that we are looking at the one-step-ahead forecast. As we provide a unified theory for all JYP's tests, and those in pitarakis2025novel are only for one-step-ahead, this choice allows us to treat them all at once. Martingale difference sequence (MDS) assumption is natural in this setting (see the same assumption in stauskas2022tests or an equivalent one in bai2006confidence, cheng2015forecasting, or karabiyik2021forecasting). However, we shall note that this is merely for convenience, as further steps ahead and serial correlation can be considered. In particular, we will comment on how to relax the MDS requirement in Remark (ref). As it is likely to occur in practice, we can also allow for conditional heteroskedasticity (e.g. $u_t=\sigma_t\varepsilon_t$, were $\varepsilon_t$ is $\mathrm{IID}(0,1)$, and $\sigma_t^2$ represents ARCH/GARCH effect) in the forecast error if ((ref)) is replaced by a consistent HAC. A.1 ii) and iv) are absent in gonccalves2017tests, but they are fairly standard and allow to bring down moment requirements on factors and idiosyncratics. Indeed, in the equivalent assumption of A.2 in gonccalves2017tests, 16th moments of factors are required. In our case, it is sufficient to have 8th moments. Generally, parts i) - iii) of A.2 ensure that the tests in pitarakis2023direct, pitarakis2025novel are asymptotically normal as required, but under lower-level conditions. For example, part iii) is similar to the one in clark2001tests and stauskas2022tests, where $\{ \mathrm{vec}(\*z_t\*z_t'-\+\Sigma_\*z)\}_{t=k_0}^{T-1}$ follows a mixing sequence of specific size (see hansen1992convergence).
Proposition (ref) below gives the said results for $g_{f,j,1}, g_{f,j,2}$ under the above assumptions, whereas in pitarakis2023direct, pitarakis2025novel they are obtained under high-level conditions.
Recall the nulls for the encompassing test: $\mathbb{E}(u_{1,t+1}(u_{1,t+1}-u_{2,t+1}))=0$, and for the forecast accuracy tests: $\mathbb{E}(u_{1,t+1}^2-u_{2,t+1}^2)=0$. Both imply that $\+\beta=\*0_r$, and thus that $g_{f,j,2}\to_p 0$. Under the alternatives instead, $g_{f,j,2}$ converges to an expression that explicitly depends on the difference between the diagonal matrix $\+I_r$ and the quadratic form of the covariances between known and unknown factors with the precision matrix of the known factors: $\+\Sigma_{\*w\*f}'\+\Sigma^{-1}_\*w\+\Sigma_{\*w\*f}$. Note that by denoting limit of $g_{f,j,2}$ by $\psi$, we can give credence to the suggested values in Algorithm (ref). For a standard normal CDF $\Phi(.)$ and a quantile related to size $\alpha$, the power function for each test is $1-\Phi(q_\alpha-\omega_j^{-1}\psi)$. As $\omega_j^{-1}\propto\gamma_j^{-1}$ for $\gamma_j$ that depends on the tuning parameters, we see that suggestions in Algorithm (ref) increase precision and boost the power.
Before we move on to the case of the estimated factors, a word on the further compact notation. We note that the components of ((ref)) - ((ref)) can be expressed as $\sum_{t=\lfloor f_1(T)\rfloor }^{\lfloor f_2(T)\rfloor}d_{u}(\widetilde{u}_{1,t+1}, \widetilde{u}_{2,t+1})$ for a loss-differential $d_{u}(.)$, where $\frac{\lfloor f_j(T)\rfloor}{T}\to q_j$ for $j=1,2$, such that $q_2>q_1$. This helps conduct analysis uniformly over different out-of-sample paths. Therefore, we will compactly formulate the results in terms of $\frac{1}{\sqrt{d_T}}\sum_{t=\lfloor f_1(T)\rfloor }^{\lfloor f_2(T)\rfloor}d_u(.)$, where $d_T=(\lfloor f_2(T)\rfloor-\lfloor f_1(T)\rfloor+1)$. They will then apply to ((ref)) - ((ref)) simultaneously.
Clearly, the estimator $\widetilde{\+\delta}_t$ in (ref) is unavailable as the factors are unobserved. Its feasible counterpart is
where $\widehat{\*z}_t=(\*w_t', \widehat{\*f_t}')'$. Therefore, to use ((ref)) in the “for” loop of Algorithm (ref), we additionally obtain $\widehat{\*F}_t$ by PC in the same loop, whose procedure is defined as
where $\*U_{Nt,r}$ is a $t\times r$ matrix of eigenvectors of $(Nt)^{-1}\*X_t\*X_t'$. For the minimization, the usual required normalizations are that ${t}^{-1}\widehat{\*F}_t'\widehat{\*F}_t=\*I_r$ and $\widehat{\+\Lambda}_t'\widehat{\+\Lambda}_t$ diagonal, as $\+\Lambda$ and $\*F$ are not identified separately. We stress that $\widehat{\+\Lambda}_t$ is indexed by $t$, because we estimate a different loading matrix for each recursion. This implies that the forecast error for the unrestricted model is given by $\widehat{u}_{2,t+1}=y_{t+1}-\widehat{\delta}'_t\widehat{\*z}_t=y_{t+1}-\widehat{\+\theta}_t'\*w_t-\widehat{\+\beta}_t'\widehat{\*f_t}$. This estimated forecast error is used in the feasible versions of ((ref)) - ((ref)).
Our goal is now to show that ((ref)) - ((ref)) have the same null asymptotic distribution and power properties as in the original references of pitarakis2023direct,pitarakis2025novel. This applies both when $\*f_t$ is observed (infeasible setting), and when the factors are estimated (feasible setting) where the loadings are allowed to be strong or weak. This means, the $r$ eigenvalues of the common component covariance matrix can diverge at a sublinear rate in $N$, i.e., at rate $N^{\alpha}$, for $\alpha\in (0,1]$. Accommodation of weaker loadings will require specific assumptions, and we follow bai2023approximate, who require that $N^{-\alpha}\+\Lambda'\+\Lambda$ has a positive definite limit. To provide intuition, $\alpha=1$ gives the usual strong loading case. Within this setting, factors and loadings can be consistently estimated with PC yielding the usual rate of $O_p(\max(1/\sqrt{N}, 1/\sqrt{T}))$. Instead, $\alpha=0$ leads to absolutely uninformative loadings similarly to onatski2012asymptotics. For instance, let $r=1$, then if $N^{-\alpha}\+\Lambda'\+\Lambda =\sum_{i=1}^N\lambda_i^2 <\infty$, we have square-summable loadings, which implies that individual loadings are practically zero for highly indexed individuals. All situations in-between give weaker (or “weakly influential") loadings de2008forecasting, onatski2012asymptotics. This means that as $N\to \infty$ the loadings are too small or too sparse for the corresponding eigenvalues to diverge at rate $N$ barigozzi2024dynamic. Because of this sublinear divergence, consistent estimation of the common component with PC is less straightforward. As we also prove later, in the recursive estimation setup, PC only allows recovering the factors associated with eigenvalues that diverge at least at rate equal to ${N}^{\alpha}$, for $\alpha >1/2$, which is a new result in the PC literature freyaldenhoven2022factor, bai2023approximate. \\ Recall that $\widehat{\*F}_t$ is the matrix of $r$ eigenvectors corresponding to the $r$ largest eigenvalues of $(Nt)^{-1} \*X_t\*X_t'$ and $\*D_{Nt,r}$ ($\*D_{Nt,r}^{2}$) is a matrix with its singular values (eigenvalues) arranged in a decreasing order.\footnote{Let us note that the SVD (and the subsequent decomposition) is here conducted on the $t\times t$ covariance for $\*X'$ i.e., $\frac{1}{Nt}\*X\*X'$, following the work of bai2023approximate. The same decomposition and asymptotic expansions below can be done in the (more traditional) case of taking eigenvectors of the $N\times N$ covariance $\frac{1}{Nt}\*X'\*X$ stock2002forecasting, barigozzi2023fnets. Both SVDs return the same set of singular values and nothing changes as long as the interest is in modeling static principal eigenvectors.} Similar decompositions as in bai2023approximate can therefore be obtained to get an explicit expression of the difference between the estimated factors and a rotated version of the true factors, i.e., $\widehat{\*F}_t-\*F_t\*H_{Nt,r}'$. This is important, as one never estimates the true factors but always a rotated version of them. Note that for the rotation matrix we have
which accommodates the fact that we employ the usual PC procedure when the loadings can be weaker. The next assumption imposes structure on the weaker loadings.
A.3 (Loadings)
We treat the loadings as random and make them independent from $\*F_t$ in order to simplify some arguments. Alternatively, we can impose $\mathbb{E}(\*f_t\*f_t'|\+\Lambda)=\*I_r$. They can also be treated as fixed, similarly to gonccalves2017tests. Moreover, as in the latter study, we impose a convergence rate to facilitate quantification of our analysis (under $\alpha=1$, it coincides with the natural $N^{-1/2}$ rate) in the recursive setup. Naturally, this slightly strengthens the original assumption in bai2023approximate. In general, i) is one of our central assumptions that allows us to conduct an asymptotic analysis uniformly in loading strength. A similar formulation is used by uematsu2022inference, but in the context of the sparse PC estimator, and by he2025huber, who exploit the minimization of Huber loss function to estimate the loadings. However, neither of these studies deals with hypothesis testing in a forecasting setting. Recently, boot2025diffusion used the assumptions of bai2023approximate to compare PC forecasts with ridge regression and random projections. The remaining assumptions are the higher-level conditions from both gonccalves2017tests (e.g. iii) and bai2023approximate (iv and v). The latter are also utilized in boot2025diffusion, but they would be directly implied if $\*f_t$ and $e_{i,s}$ were independent for all $i,s,t$. The equivalent of iii) in gonccalves2017tests requires 16-th moment, whereas we deem it very strong and bring it down due to A.1 iv).
We obtain an expansion of $\widehat{\*F}_t-\*F_t\*H_{Nt,r}'$ by implementing further modifications from both bai2023approximate and gonccalves2017tests. We start with the case of homogeneous loadings and move in the next section to the case of heterogeneous loadings, i.e., when the loadings are allowed to be weaker to different degrees. We define the following scalar quantities:
Note that only the terms that depend on $\+\lambda_i$ are directly scaled by $N^{-\alpha}$. Since ((ref)) and ((ref)) are functions of the idiosyncratics only, $N/N^{\alpha}$ can be seen as a “penalty” term on the overall rate, as we estimate potentially weaker loadings with PC. Therefore, we obtain the following decomposition for a row $s$:
This is an important term appearing when expanding $\widetilde{u}_{2,t+1}- \widehat{u}_{2,t+1}$, that is, the difference between the infeasible and feasible estimated forecast errors. Specifically, let $\+\Phi_{Nt,r}=\mathrm{diag}(\*I_k, \*H_{Nt,r})\in \mathbb{R}^{(k+r)\times (k+r)}$, where $\*H_{Nt,r}$ is the rotation matrix defined in ((ref)). Then, in line with gonccalves2017tests and stauskas2022tests, we obtain
Because it is a scalar, we have that $III=(\widehat{\*f}_t-\*H_{Nt,r}\*f_t)'(\*H_{Nt,r}^{-1})'\+\beta$, and $III$ is therefore absent when $\+\beta=\*0_r$ (under the null). The component $I$, instead, reveals that the asymptotic equivalence is ensured when the factors are estimated consistently and the infeasible forecasting model is well-specified ($\|\widetilde{\+\delta}_t-\+\delta\|=o_p(1)$). Term $II$ additionally requires that the feasible OLS estimator of the parameters is asymptotically equivalent to the infeasible one. Because $\widehat{\+\delta}_t$ employs the factors which are only identified up to a rotation, the former is naturally rotated, as well. \footnote{Indeed, by the Frisch-Waugh-Lovell argument, the second component is $(\*H_{Nt,r}^{-1})'\widetilde{\+\alpha}_t=(\*H_{Nt,r}'\*F_t'\*M_\*W\*F_t\*H_{Nt,r})^{-1}\*H_{Nt,r}'\*F_t'\*M_\*W\*y_t$, where $\*M_\*W$ is the projection matrix onto the orthogonal complement of the observed predictors.}
\fi
The expression ((ref)) is the integral part of the feasible versions of ((ref)) - ((ref)). To demonstrate their asymptotic equivalence, it is useful to introduce the following quantities:
which are expressed in the compact notation in order to be applied to every statistic under consideration. Then, for $j=1,\ldots,4$
where $q_j(.)$ is a function, such that $|q_j(A,B,C,D) |=o_p(1) $ if $A,B$ and $C$ are negligible and $\widehat{\omega}_j=\widetilde{\omega}_j+o_p(1)$. Therefore, the asymptotic equivalence holds for $j=1,\ldots,4$ if ((ref)) - ((ref)) are negligible. Expansion of all the infeasible statistics can be found in the Online Supplementary material.
The structure of ((ref)) - ((ref)) reveals how the behavior of the out-of-sample average of ((ref)) dictates the overall asymptotic analysis. For this, we state more assumptions regarding the idiosyncratic components.
A.4 (Idiosyncratics)
A.4 is similar to conditions applied in gonccalves2017tests, except for the lower moment requirement. It allows the idiosyncratics to be weakly dependent over time and cross-sectionally. Part vii) is the same as in bai2023approximate, but it is required to hold uniformly in $t$. We formulate it in terms of $T$, because $t=\lfloor s T\rfloor$ for some $s\in (0,1)$ due to the recursive setup. It is primarily used to simplify proofs when the loading weakness is heterogeneous, but it also helps to improve convergence rates.
Lemma (ref) below formalizes the asymptotic behavior of the out-of-sample average of ((ref)) by providing the rate of the average square factor approximation error over the recursive samples, when the loadings are weaker. In addition, it provides the uniform rate of the in-sample factor approximation error.
Lemma (ref) can be seen as an extension of Theorem 4.1 in gonccalves2017tests to the case of weaker loadings and different alternative out-of-sample paths. If $\alpha=1$ (strong loadings), we return to the usual rate of $O_p(\max\{1/N, 1/T\})$, which coincides with the result in the latter study. Note that the result immediately implies that we must have $\frac{N}{N^\alpha}\frac{1}{\sqrt{T}}\to 0$ and $\alpha>0.5$ to consistently estimate the factor space in both cases of in- and out-of-sample. In contrast, the theory in bai2023approximate requires a lower bound of $\alpha$ different from 0 only for inference exercises, but not consistency. The difference arises, because in part (i) we consider a recursive setup, where the rotation matrix in ((ref)) changes for every $t=k_0,\ldots,T-1$, and therefore this extra informativeness condition needs to be satisfied. Part (ii.) of Lemma (ref) deals with the average square in-sample factor estimation error, where we average over $s=1,\ldots, t$ as $ \frac{1}{t}\left\|\widehat{\*F}_t-\*F_t\*H'_{Nt,r} \right\|^2= \frac{1}{t}\sum_{s=1}^t\left\| \widehat{\*f}_s - \*H_{Nt,r}\*f_s\right\|^2$. We use $\widehat{\*F}_{t-1}$ when obtaining the feasible $\widehat{\+\delta}_{t}$ for each recursion. Interestingly, $\alpha>0.5$ is sufficient, but not necessary for part (ii) to hold, since we average for a given $\*H_{Nt,r}$. Hence, the approximation rate can be improved with higher-level conditions in both homogeneous and heterogeneous cases, as we point out in Remark (ref). However, part (i) is responsible for the general out-of-sample approximation and will determine the behavior of $g_{\widehat{f},j}$ for $j=1,\ldots, 4$. \\ The following result employs (ii.) of Lemma (ref) to demonstrate the uniform equivalence of feasible and infeasible OLS estimators.
Similarly to Lemma (ref), Lemma (ref) can be seen as a generalization of Lemma 4.1 of gonccalves2017tests to weaker loadings. This follows from the direct application of our Lemma (ref) (ii) when establishing the uniform consistency of $\widehat{\+\delta}_t$ for (the rotated) $\widetilde{\+\delta}_t$. Hence, $\alpha>0.5$ plays a role, as well. The implication of the lemma is that the rate of consistency is $o_p(T^{-1/4})$, unlike in the latter study, where it is $o_p(T^{-1/2})$. The difference arises through two related channels. Firstly, the tests of pitarakis2023direct, pitarakis2025novel have different local power properties. Indeed, we specify $\+\beta =\+\beta^0T^{-1/4}$, while the tests of clark2001tests explored in gonccalves2017tests require $\+\beta =\+\beta^0T^{-1/2}$. Secondly, we chose a different proving technique, since we only needed to demonstrate consistency, but the exact rate is less important due to different asymptotic properties of our tests. \\ In order to move on to the main results, we introduce the last assumption.
A.5 ($N,T$ Expansion Rates)
i)\;$\frac{N}{N^{\alpha}}\frac{1}{T^{1/4}}\to c>0$, as $(N,T)\to \infty$,\quad ii)\;$\sqrt{T}N^{-\alpha}\to 0$, for $\alpha\in (0.5, 1)$.
A.5 i) should be seen as a device to make convergence rates more transparent when the loadings are indeed weaker. Specifically, to prove that $B$ and $C$ in ((ref)) and ((ref)) are negligible, we will need $\frac{N}{N^{\alpha}}\frac{1}{\sqrt{T}}=o(1)$ as before, but at the same time $T^{1/4}\frac{N}{N^{\alpha}}\frac{1}{\sqrt{T}}=O(1)$. Thus, under A.4 i), $\frac{N}{N^{\alpha}}\frac{1}{\sqrt{T}}=O_p(T^{-1/4})$, which is the rate that can be incorporated in the subsequent analysis. Requirements in spirit of A.5 ii) are often met in the PC literature, and it helps to ensure that factor estimation error does not accumulate too fast with the expansion of $N$ and $T$. Here, it is identical to the assumption needed for inference in Lemma 4 of bai2023approximate. Under $\alpha=1$ (strong loadings), this naturally coincides with the requirement of $\sqrt{T}N^{-1}=o(1)$ in gonccalves2014bootstrapping or gonccalves2017tests. Both here and in the latter study this requirement targets sums over an MDS process $\{ (\widehat{\*f}_t-\*H_{Nt,r} \*f_t)u_{t+1}\}_{t=k_0}^{T-1}$ and makes sure that they remain asymptotically negligible. Such terms appear by applying ((ref)) to components $B$ and $C$ in (ref) and (ref), respectively.\\ Lemma (ref) below is the central outcome that utilizes the interim results discussed above.
To our knowledge, this is the first result that controls factor estimation error in the out-of-sample context uniformly in loading strength and alternative out-of-sample paths. Apart from the desired results on $A,B$ and $C$, we can also see that the feasible variance estimator is asymptotically equivalent to the infeasible one, because
under our assumptions. This provides the last missing piece to establish the equivalence result.
The main message of Theorem (ref) is that we are able not only to use the battery of new statistics ((ref)) - ((ref)) in the popular context of factor-augmented forecasts, but that they are also robust to weaker factor loadings, as long as $\alpha>0.5$. This result is a companion to Theorem 4.2 in gonccalves2017tests. While it does not necessarily nest their results due to different statistics and their local power properties, we have the same approximation rates under strong loadings ($\alpha=1$).
Before we generalize our results to the heterogeneously weak loadings, it is important to illustrate how Theorem (ref) can be used to improve the statistics. The tests of equal forecast accuracy in ((ref)) and ((ref)) depend on the tuning parameters $\lambda_1^0$ and $\lambda_2^0$. Naturally, while suggestions on their values are provided in Algorithm (ref), your choice still alters the power properties of the tests, as argued in pitarakis2025novel. To bypass this issue, power-enhanced versions of the statistics are presented, e.g.
where $\widetilde{\zeta}(\lambda_1^0,\lambda_2^0)=\frac{1}{\widetilde{\omega}_2}\frac{1}{\lambda_2^0}\frac{1}{\sqrt{n}}\sum_{t=k_0}^{k_0+l_2^0-1}(\widetilde{u}_{1,t+1}-\widetilde{u}_{2,t+1})^2$ in spirit of fan2015power. Clearly, the power-adjustment term is infeasible and we must replace it with $\widehat{\zeta}(\lambda_1^0,\lambda_2^0)$. Proposition (ref) below demonstrates that the power enhancement procedures remain valid.
The adjustment term $g^{adj}_{\widehat{f},4}$ concerns our new statistic $g_{\widehat{f},4}$ in ((ref)). We relegate its asymptotic analysis together with a broader discussion on how the power adjustment terms are constructed to Section 4.4 in the Online Supplementary material.
\bigbreak
What done in the previous section imposed that the loadings are weaker to the same degree, or “homogeneously". Alternatively, we can allow for “heterogeneously" weak loadings by employing the normalizing matrix $\*B_N=\mathrm{diag}(N^{\alpha_1/2},\ldots, N^{\alpha_r/2})$, where $1\geq \alpha_1>\alpha_2>\ldots >\alpha_r>0$, where the weakest loading cannot still be absolutely uninformative. This means that as $N\to \infty$ some loadings are too small or too sparse for the corresponding eigenvalues to diverge at rate $N$, while others are relatively stronger or right-out strong such that the corresponding eigenvalues diverge at slighlty sub-linear rate or linear rate in $N$. Clearly, the heterogenous case nests the homogeneous one when $\alpha_1=\cdots=\alpha_r.$ Note that $\|\*B_N \|\leq MN^{\alpha_1/2}$ and $\|\*B_N^{-1} \|\leq m N^{-\alpha_r/2}$ for some positive constants $M$ and $m$, which means that the order of (the inverse of) this normalization matrix is dominated by the (weakest) strongest factor loading. Again, recall that $ \widehat{\*F}_t\*D_{Nt,r}^2=\frac{1}{Nt}\*X_t\*X_t'\widehat{\*F}_t$ by the eigenvalue-eigenvector relationship. Then, by using the fact that both $\*B_N$ and $\*D_{Nt,r}^2$ are diagonal, we obtain $\widehat{\*F}_t\*D_{Nt,r}^2\*B_N^{-1}=\widehat{\*F}_t\*B_N^{-1}\*D_{Nt,r}^2 =\widehat{\*F}_t\*B_N(\*B_N^{-2}\*D_{Nt,r}^2)=\frac{1}{Nt}\*X_t\*X_t'\widehat{\*F}_t\*B_N^{-1}.$
As for the above, we are after an expansion of $\widehat{\*F}_t-\*F_t\*H_{Nt,r}'$.\footnote{To be precise, we are after the expansion of $(\widehat{\*F}_t\*B_N-\*F_t\*B_N\overline{\*H}_{Nt,r}')\*B_N^{-1}$, for $\overline{\*H}_{Nt,r}:=(N\*B_N^{-2}\*D_{Nt,r}^2)^{-1}\*B_N^{-1}t^{-1}\widehat{\*F}_t'\*F_t\+\Lambda'\+\Lambda \*B_N^{-1}.$ However, as we explain in the Online Supplement, the component $\*B_N^{-1}t^{-1}\widehat{\*F}_t'\*F_t\+\Lambda'\+\Lambda \*B_N^{-1}$ is bounded in probability for $t=T$ as argued in bai2023approximate, and we also show that it is uniformly bounded in Lemma 1 of the mathematical Online Supplement. Furthermore, using the fact that the product of $\*B_N$ and $\*D_{Nt,r}^2$ commutes we can show how $\*B_N^{-1}\overline{\*H}_{Nt,r}\*B_N=\*H_{Nt,r}.$} Hence, we will use the same definitions of the scalar quantities in ((ref)) - ((ref)). However, we replace $\alpha$ with $\alpha_r$ (in the interim, absence of a subscript $\alpha$ ($\alpha_r$) means that there is no scaling in terms of $N$). Next, to accommodate heterogeneity in loading weakness, we re-define
where $\left\|N^{-1/2}\*B_N \right\|=O(N^{(\alpha_1-1)/2})$ that is $O(1)$ when $\alpha_1=1$. Additionally, we introduce $\*Q_{N,\alpha_r}=\*B_N^{-1}\sqrt{N^{\alpha_r}}$, which is $O(1)$, as well. Then, by again using the commutative product, we get
\bigbreak
To handle heterogeneous weak loadings we require adaptations of the former assumptions A.3, A.5. More specifically:
A.3* (Heterogeneously weak loadings)
A.5* ($N,T$ Expansion Rates for heterogeneously weak loadings)\\ i)\;$\frac{N}{N^{\alpha_r}}\frac{1}{T^{1/4}}\to c>0$, as $(N,T)\to \infty$, \quad ii)\;$\sqrt{T}N^{-\alpha_r}\to 0$, for $\alpha_r\in (0.5,1)$. \\
Finally, Lemma (ref), whose proofs are in the Online Supplement 3.2.2 & 4.2, links the results of Lemma (ref), Lemma (ref) and Lemma (ref) to the heterogeneous loadings context. With this, results in Theorem (ref) follow directly.
We design a DGP similar to DGP2 in pitarakis2025novel but where factors are specified in the same way as DGP2 of bai2023approximate with strong/ weak homogeneous/ heterogeneous loadings. Throughout, we set the number of factors $r=3$.
where: $\*G_i\sim N(0,\*I_3)$, $\*D^2=\operatorname{diag}(3\;2\;1);$ $ \*B_N=\operatorname{diag}(N^{\alpha_1/2}\;N^{\alpha_2/2}\;N^{\alpha_3/2})$, $(\alpha_1,\alpha_2,\alpha_3)=(1, 1, 1)$ for strong homogeneous loadings, $(0.51, 0.51, 0.51)$ for weak homogeneous loadings, $(0.51, 0.7, 1)$ for mixed strong/weak heterogeneous loadings. \footnote{Note that $\+\lambda_i$ is simulated differently than in bai2023approximate. The reason is that the rate of $\|\*B_N^{-1}\+\Lambda'\+\Lambda\*B_N^{-1}- \+\Sigma_{\+\Lambda} \|$ plays an important role in our asymptotic analysis, whereas it did not matter for bai2023approximate. This simulation method mimics A.3 i), because we can show that under such design $\*B_N^{-1}\+\Lambda'\+\Lambda\*B_N^{-1}= \*D^2 + O_p(N^{-\alpha_r/2})$ as desired. } Cross-sectional and time dimensions are: $(N,T)=(800,500)$. For the practical implementations, following pitarakis2025novel, we set $c=1.25$, $\theta_1=0.5$, $\rho_i=0.3+N(0,1)_i\times 0.5$, $\+\beta=(0,0,0)'$, $\pi=24$, for size, and $\+\beta=(j,j,j)'$ for $j\in \{0.1, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6\}$ for power. Before we go to the results, let us mention how we performed the same simulations using bai2002determining criterion “$IC_{p1}$" (see their eq. 9) to actively select the number of factors. We select this number once during the in-sample period to reflect the assumption that $r$ is fixed over time. Since there is no recursion in the selection of the number of factors, this choice is justified by the results in bai2023approximate who find that in order to estimate factors with weakly convergent loadings ($\alpha>0$), the criteria in bai2002determining remain valid. This turned out to be identical to the results presented here as $IC_{p1}$ always correctly estimates the number of factors.\footnote{We shall mention, however, that others of the bai2002determining criteria did not perform as well as $IC_{p1}$ in selecting the factors; we reported $IC_{p1}$ as it is the most remarkable.}
We report here a summary of the Monte Carlo in the form of power curves for the different test statistics considered. The setting is baseline, meaning that $u_{t+1}$ is uncorrelated over time and $e_{i,t}$ is uncorrelated cross-sectionally. In Figure (ref) - Figure (ref), we find both Encompassing and Forecast Accuracy tests to display satisfactory sizes and powers when factors are included in the alternative forecasting model. In terms of choice of the parameters $\mu_0, \tau_0, \lambda_1^0, \lambda_2^0$ we report here the best performing size-wise (i.e., closest to nominal level $5\%$, see Online Supplement Section 5 for the extended results). As expected, weaker loadings have a dampening effect on power, especially, of course, when the signal ($\+\beta$ value) is low. This effect is more pronounced when looking at the Forecast Accuracy tests $g_{\widehat{f},2}, g_{\widehat{f},2}^{adj}, g_{\widehat{f},3}, g_{\widehat{f},3}^{adj}$, $g_{\widehat{f},4}$, $g_{\widehat{f},4}^{adj}$, which are already affected --power-wise-- by the data-loss due to the sample overlapping discussed earlier. The power adjustment for $g_{\widehat{f},2}, g_{\widehat{f},3}, g_{\widehat{f},4}$ is paramount as the unadjusted versions can lead to severely undersized tests, as can be seen in the tables in Online Supplement Section 5. The $g_{\widehat{f},4}^{adj}$, which averages over the second coordinate thus to let the MSE of the unrestricted model accumulating, has an entirely analogous behavior as $g_{\widehat{f},3}^{adj}$ with some slightly higher power for lower $\+\beta$'s. Overall, we show how all the tests have good finite sample performances when factors are included in the alternative model specification, and loadings can be either strong or weak, homogeneous or heterogeneous. We refer to Online Supplement Section 5 for all the results, including the heterogeneous ones.
In Figure (ref) - Figure (ref) we introduce weak cross-section dependence in $e_{i,t}$ similarly to stauskas2022tests. We see that the results are virtually the same as in the baseline scenario. This signals robustness to both temporal and, for example, spatial dependence structures in idiosyncratics.
In Figure (ref) - Figure (ref) we introduce, on top of the cross-section dependence of the idiosyncratics, conditional heteroskedasticity in the forecast errors $u_{t+1}$, in the form of a GARCH$(1,1)$ as specified above. Throughout the simulations, we still use ((ref)) to estimate $\phi^2$. The main effect for all tests is a mild inflation of the size\footnote{As one would expect, if the persistence of the past conditional variance of $u_t$ is high(er), i.e., the $\eta$ of the GARCH$(1,1)$ is large(r), the size would suffer more. In such cases HAC-type corrections of the variance are recommended.}, which is however milder for $g_{\widehat{f},4}^{adj}$ compared to $g_{\widehat{f},1}$, $g_{\widehat{f},2}^{adj}$ and $g_{\widehat{f},3}^{adj}$. The power is instead higher for lower values of the coefficient $\boldsymbol{\beta}$ (e.g., $\boldsymbol{\beta}=[0.10,\ldots,0.30]$) if compared to the previous results, though this might just be a byproduct of the size increase or of the GARCH time varying variance structure, which potentially makes certain periods in the sample more informative. At the same time, even though ((ref)) estimates unconditional variance, it can be sensitive to conditional heteroskedasticity, especially as it utilizes 4th moments. While a theoretical justification for HAC-type corrections in our testing framework is left for future research, we have experimented with the practical implementations proposed by pitarakis2023direct for $g_{\widehat{f},1}$, specifically using Newey-West and Andrews standard errors. These methods demonstrate good empirical performance in our simulations and are thus recommended for practitioners, particularly in settings where conditional heteroskedasticity—such as GARCH-type volatility—is likely present in the forecast errors. The Online Supplement contains heterogeneous weakness results, as well.
{0.7cm} We partially replicate the inflation forecast exercise conducted by pitarakis2023direct. Our focus here is to explore whether global inflation can enhance the accuracy of country-level inflation forecasts. This debate is not new; previous studies have provided evidences that global inflation trends can significantly improve domestic inflation forecasts monacelli2009international, ciccarelli2010global. However, other research suggests that the relevance of global inflation in forecasting domestic rates may stem solely from its ability to capture slow-moving trends in inflation mikolajun2016advanced, or that a global inflation factor improves forecasting accuracy primarily at longer horizons gillitzer2019does. The question is, also, how to measure global inflation. One approach is to calculate a grand average of country-level inflation rates, represented as $\Bar{\pi}_t=N^{-1}\sum_{i=1}^N \pi_{i,t}$. Another method involves treating global inflation as a few latent factors that can be estimated using PC from the pool of country-level inflation rates. The former approach is employed by pitarakis2023direct to illustrate his encompassing test, while we will utilize the latter method to demonstrate how JYP's encompassing tests work in an empirically relevant context where PC factors are considered. There are several reasons why PC global inflation serves as a more accurate measure than a simple grand average of inflation rates. A straightforward sample mean does not consider differences in economic size, inflation volatility, or other factors that may make some countries' inflation rates more indicative of global trends than others. In contrast, PCs can uncover patterns of co-movement in inflation rates that might not be evident from the raw data. Additionally, it is more robust to potential outliers and can adapt to time-varying relationships among countries' inflation rates. However, one could argue that if we assume the existence of only one factor, then all the information contained in that factor is effectively the same as a cross-sectional average of the countries' inflation rates. This would be true, upon essentially three assumptions: (i) the existence of an exact factor model underlying the data, i.e., $\pi_{i,t}=\+\lambda_i'\+f_t+e_{i,t}$ with $Cov(e_{i,t}, e_{j,s})=0$, $t,s \in \mathbb{Z}$, $i,j=1,\ldots,N$, $i\neq j$; (ii) a large cross-sectional dimension $N$ (in principle $N\to \infty$); (iii) all or most of loadings being non-zero (i.e., pervasiveness of factors). Assume that the loadings are fixed. If all these are satisfied, it is clear how $Var(\Bar{\pi}_t)=N^{-2}(\sum_i \lambda_i)^2 Var(f_t)+N^{-2} Var(\sum_i e_{i,t})\to \Bar{\lambda}^2=Var(\Bar{\chi}_t),$ as $N\to \infty,$ meaning how the aggregation of the observed data recovers the same information contained in the factor. Now, (i) is clearly too strong (see also our A.4, v), (ii) is what is referred to as the “blessing of dimensionality" but in practice it clearly depends on the available data, (iii) is precisely what is challenged by the weaker loadings treated in Section 2. Hence, there are good reasons to re-run this exercise using factors and employing the JYP's tests to check if global inflation computed by means of PC factors improves the country-level inflation forecast. We employ the same dataset provided in pitarakis2023direct\footnote{Freely available on JYP's GitHub page: \url{https://github.com/jpitarakis/Multi-Step_Encompassing}}, based on the World Bank global inflation database covering the period 1970-2023 for 23 countries at quarterly frequency. We choose the tuning parameters in line with the simulations to balance the power and size. \footnote{For the exact treatment and transformations of the raw data we refer to pitarakis2023direct, Section 7.}
The base-line model is an $AR(1)$ while the alternative model is a factor-augmented $AR(1)$ where factors are recursively estimated and their number is determined via bai2002determining information criterion ($IC_{p1}$, max number$=10$).\footnote{It is well known how the bai2002determining criteria depend quite substantially on the maximum number of factors as selected by the practitioner, as well as the relative magnitude of $N,T$ forni2009opening. In this case, we find on average that all criteria, including $IC_{p1}$ but excluding $AIC_3$, return the maximum as estimated number of factors. We experimented using $AIC_3$ and the ABC criterion of alessi2010improved too, which both on average estimate 4/5 common factors, but we found a completely similar picture in terms of significance, with only slightly higher p-values.} We find patterns of significance across the four test statistics, though rarely all of them at once (only AUT). $g_{\widehat{f},4}^{adj}$ finds significant better forecast accuracy one quarter ahead when using PC global inflation in $9$ of the $23$ countries. It is followed by the encompassing test $g_{\widehat{f},1}$ with $8$ countries ($9$ with CAN if considering $10\%$ nominal level), $g_{\widehat{f},3}^{adj}$ with $7$ and $g_{\widehat{f},2}^{adj}$ with $4$. Noticeably, for some large economies such as USA, FRA, ESP, at least two test statistics of the four are found significant (JPN too if considering $10\%$ nominal level). Some other large/medium-to-large economies like GBR, DEU, ITA have only one or no significance at all. Overall, this paints a mixed picture with regard to the use of global inflation to better forecast country-level inflation, which seems to be clearly a country-specific issue. These tests can therefore be used to gain a sense as to whether PC global inflation might help beyond the simple AR specification. In Online Supplement Section 6 we repeat the exercise without factors but using a grand average of country level inflations (which excludes the country of reference each time) as global inflation. What comes out is an overly suspicious abundant significance across most countries and all tests. Outliers and multicollinearity are likely to affect these results, as evident from the boxplots and pairwise correlation heatmap in Figure (ref) and Figure (ref), respectively. It is therefore safe to say that the analysis including PC global inflation is much more trustworthy.
We continue the empirical exercise by considering more macroeconomic series, this time from an even larger pool of series. Table (ref) collects the p-values for the tests of encompassing and forecast accuracy described in Section (ref), using the same tunings as in the previous application in Section (ref), for 14 important macroeconomic series taken from the FRED-MD dataset mccracken2016fred (July 2024 vintage, starting January 1960). The FRED-MD dataset includes a wide range of U.S. macroeconomic series such as output (e.g., industrial production), income, labor market indicators (e.g., unemployment rate, payrolls), consumer and producer prices, money supply, interest rates, exchange rates, and financial variables like stock indices and spreads. These series capture economic activity, inflation dynamics, monetary policy, and financial market conditions. All series are cleaned and stationary transformed as prescribed by the Matlab routines provided with FRED-MD. After the necessary cleaning of NAs and outliers, the final dataset contains $117$ series for $772$ data points. Again we choose the tuning parameters in line with the simulations to balance the power and size (see our extensive simulations in the Online Supplement).
We forecast one-month-ahead these $14$ macroeconomic series using the PC factors in the alternative forecasting specification (i.e., $AR(1)$ vs factor augmented $AR(1)$, where factors are estimated via bai2002determining $IC_{p1}$ with max number$=15$).\footnote{We choose to follow a common convention in the factor modeling literature by setting a relatively high maximum number of factors, from which the optimal number is selected using IC. The goal of this is to ensure that the true number of factors is not underestimated due to an arbitrarily restrictive upper bound.}\footnote{We collect the full names of the $14$ FRED-MD macroeconomic time series in the Online Supplement (Section 6, Table 6.2)} With an average of $7$ estimated common factors, the results are overwhelmingly significant across the four tests. Apart from “Real personal consumption expenditures" (DPCERA3M086SBEA), all the other series have at least one significant (at nominal level 5%) test statistic. This shows how even in a high(er) dimensional dataset such as the FRED-MD, these test statistics are very useful in assessing whether PC factors could be employed to improve the forecast of macroeconomic series.
{0.7cm}
We develop the theoretical framework that allows the forecast accuracy and encompassing tests proposed by pitarakis2023direct, pitarakis2025novel to be applied when the alternative forecasting model incorporates estimated PC factors. These factors can have loadings that are either strong or weak, whether homogeneously or heterogeneously distributed. Our theoretical findings are supported by both a Monte Carlo simulation and two empirical applications in macroeconomic forecasting.