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,848 characters · 18 sections · 112 citation commands
Diffusion Index Forecasting with Tensor Data
\if00 {
\affil[1]{University of Rochester} \affil[2]{University of Notre Dame} } \fi
\if10 \fi
Since the seminal work of StockWatson2002 and BaiNg2006, diffusion index forecasting has been widely adopted by government agencies, policy institutes, and academic researchers around the world (see, e.g., LudvigsonNg2007,LudvigsonNg2009,Kyle2015). The classical diffusion index model predicts the target variable as a linear combination of factors extracted from a large panel of time series data, as well as other important predictors. Its strength lies in the ability to significantly reduce the dimensionality of the predictor space by summarizing it into a small number of factors, enabling effective use of large datasets while keeping the size of the forecasting model small.
However, as the availability and complexity of economic data have expanded, new challenges have emerged for forecasting models. In particular, multidimensional data, panel data with more than two dimensions, have attracted increasing attention in economics due to their ability to capture richer and more intricate relationships. For example, consider predicting U.S. import/export volumes with China using monthly time series data. While the traditional gravity model focuses on bilateral trade flows, it may overlook the influence of trade patterns between the U.S., China, and other countries due to substitution effects. Such data can be structured as a three-dimensional tensor, where the observed time series $\mathcal{X}_t$ is of dimension $N \times N$ for each period $t$, with $N$ denoting the number of countries in the dataset. This type of multidimensional structure poses challenges to classical diffusion index forecasting, which is based on vector factor models.
A natural approach to tackling this challenge is to flatten or vectorize the tensor time series (see, e.g., LudvigsonNg2007) to fit within the framework of vector factor models. However, this process changes the original data structure, potentially diminishing the interpretability of how information from different dimensions interacts. Furthermore, vectorizing tensors often leads to a significant increase in the number of parameters to estimate, which can result in high computational costs.
In this paper, we consider diffusion index forecasting with tensor and non-tensor predictors, where the tensor structure is preserved with a Canonical Polyadic (CP) tensor factor model\footnote{CP and Tucker structures are the two most commonly assumed low-rank structures for tensor factor models (see, e.g. KoldaBader2009). We adopt the CP low-rank structure due to its parsimonious features.}. Common factors are extracted from tensor data using the contemporary covariance-based iterative simultaneous orthogonalization (CC-ISO) procedure proposed in infCP2024. When the number of potential non-tensor predictors is small, we estimate the diffusion index model with ordinary least squares (OLS) and establish the consistency and asymptotic normality of the estimator. Unlike BaiNg2006, we allow factors to exhibit different strengths. The convergence rate of the conditional mean prediction for the target variable depends on both the strength of the weakest factor and the sample size. To conduct valid inferences, we propose a thresholding-based covariance matrix estimator that is robust to cross-sectional correlation in the idiosyncratic component and demonstrate its consistency.
When the number of potential non-tensor predictors is large, potentially comparable to or exceeding the sample size, we propose a two-step penalized regression approach, applying least absolute shrinkage and selection operator (LASSO) to select important non-tensor predictors. The combination of factor models and sparse regression has been explored in the literature. In a panel data context, Fanfactorandsparse consider the factor-augmented sparse linear regression model, which includes a vector latent factor model and sparse regression as special cases. chenfanzhu2024 extend this framework to matrix-variate data and propose two new algorithms for estimation. However, both papers focus on a single type of predictor, either panel or matrix data. In economic forecasting, researchers often have access to mixed types of data. Consider the trade example again. While trade flows among various countries provide valuable information for predicting U.S. import/export volumes, other economic variables such as GDP, unemployment rates, exchange rates, and interest rates also play a critical role. These different types of data may reflect distinct sources of predictability. The tensor data on trade flows captures global factors while macroeconomic variables act as proxies for local predictability. Our model offers a novel framework for integrating these diverse data sources to improve forecast accuracy.
Our work also relates to several recent developments in econometrics. Within tensor and matrix factor models, recent contributions include BabiiGhyselsPan2025, BeyhumGautier2022, bolivar2025threshold, chen2021statistical, infCP2024, chen2024time, han2022Rank, han2024tensor, yu2024dynamic, among many others. Our framework differs by focusing on diffusion-index forecasting and integrating both tensor and non-tensor predictors within a unified structure. From the perspective of factor-augmented regressions, classical results often find factor estimation to be first-order neutral for OLS inference (StockWatson2002, BaiNg2006 and CaiKongWuZhao2025), though it can matter in some important cases (GONCALVES2014156). We also contribute to the growing literature on high-dimensional covariance estimation (Bickel2008,Rothman2009,Fan2013) and on LASSO methods for dependent data (KOCK2015325, MEDEIROS2016255, ChernozhukovHardleHuangWang2021, BabiiGhyselsStriaukas2022, BabiiGhyselsStriaukas2024 and Beyhum2024). Finally, our empirical application on forecasting international trade follows the standard macro-forecasting tradition, where autoregressive (AR), vector autoregressive (VAR), and diffusion-index models (StockWatson2002) serve as common benchmarks for evaluating forecast performance.
The rest of the paper is organized as follows. In Section 2, we introduce the diffusion index model based on the CP tensor factor model and develop the estimator when the number of non-tensor predictors is small. Section 3 derives the inferential theories for the diffusion index model and proposes a robust covariance matrix estimator. Section 4 introduces multi-source factor-augmented sparse regression to combine information from different sources and discusses the consistency of the proposed estimator. In Section 5, a simulation study is conducted to assess the reliability of the low- and high-dimensional estimators in finite samples. In Section 6, an empirical example on U.S. export/import forecasting highlights the merits of our approach in comparison with some popular methods in the literature. All mathematical proofs and additional simulation results are contained in the Appendix.
In this subsection, we introduce essential notations and basic tensor operations. For an in-depth review, readers may refer to KoldaBader2009.
Let $\|x\|_q = (x_1^q+\cdots+x_p^q)^{1/q}$, $q\ge 1$, for any vector $x=(x_1,\cdots,x_p)^\top$. In particular, $\| x \|_\infty = \max_{1 \leq j \leq p} | x_j |$. We employ the following matrix norms: matrix spectral norm $\|M\|_{2} = \underset{\|x\|_2=1,\|y\|_2=1}{\max} |x^\top M y| = \sigma_1 (M)$, where $\sigma_1(M)$ is the largest singular value of $M$; max entry norm: $\left\| M \right\|_{\max} = \max_{1 \leq i \leq p, 1 \leq j \leq q} | M_{ij} |$ for $M \in \mathbb{R}^{p \times q}$, where $M_{ij}$ denotes the $(i,j)$ entry of $M$. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$, we write $a_n\lesssim b_n$ (respectively, $a_n\gtrsim b_n$) if there exists a constant $C$ such that $|a_n|\leq C |b_n|$ (respectively, $|a_n|\geq C |b_n|$) holds for all sufficiently large $n$, and $a_n\asymp b_n$ if both $a_n\lesssim b_n$ and $a_n\gtrsim b_n$ hold.
Consider two tensors ${\cal A}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}, {\cal B}\in \mathbb{R}^{p_1\times p_2\times \cdots \times p_N}$. The outer product $\otimes$ is defined as ${\cal A}\otimes {\cal B}\in \mathbb{R}^{d_1\times \cdots \times d_K \times p_1\times \cdots \times p_N}$, where $$({\cal A}\otimes{\cal B})_{i_1,...,i_K,j_1,...,j_N}=({\cal A})_{i_1,...,i_K}({\cal B})_{j_1,...,j_N} .$$ The mode-$k$ product of ${\cal A}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}$ with a matrix $U\in\mathbb{R}^{d_k\times r_k}$ is an order $K$-tensor of size $d_1\times \cdots \times d_{k-1} \times r_k\times d_{k+1} \times \cdots \times d_K$, denoted as ${\cal A}\times_k U^\top$, where $$ ({\cal A}\times_k U)_{i_1,...,i_{k-1},j,i_{k+1},...,i_K}=\sum_{i_k=1}^{d_k} {\cal A}_{i_1,i_2,...,i_K} U_{j,i_k}. $$
Given ${\cal A}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}$ and a sequence of $\{U_k\}_{k=1}^K$, where $U_k \in \mathbb{R}^{d_k \times r_k}$, the notation ${\cal A} \times_{k=1}^K U_k^\top$ denotes a sequence of mode-$k$ product: $$ {\cal A} \times_{k=1}^K U_k^\top = {\cal A} \times_1 U_1^\top \times_2 U_2^\top \times \cdots \times_K U_K^\top \in \mathbb{R}^{r_1 \times r_2 \times \cdots \times r_K}. $$
The Khatri-Rao (or column-wise Kronecker) product of two matrices $A=(a_1,a_2,\cdots,a_r)$ and $B=(b_1,b_2,\cdots,b_r)$ is defined as $A * B = (a_1 \odot b_1, \cdots, a_r \odot b_r)$, where $\odot$ denotes the Kronecker product. Denote $d=d_1\times d_2\times \cdots \times d_K$, $d_{\min}=\min_{k\leq K} d_k$ and $d_{\max}=\max_{k\leq K} d_k$.
Assume that a decision maker is interested in predicting some univariate series $y_{t+h}$, conditional on $I_t$, the information available at time $t$, which consists of a tensor-variate predictor ${\cal X}_t\in \mathbb{R}^{d_1 \times d_2 \times \cdots \times d_K}$ and a set of other observable variables $w_t\in \mathbb{R}^{p}$, such as lags of $y_t$. We consider a diffusion index model as
where $h \geq 0$ is the lead time between information available and the target variable. The vector $f_t = (f_{1t}, \ldots, f_{rt})^\top$ consists of $r$ latent factors extracted from the observed tensor data ${\cal X}_t$. Specifically, we model ${\cal X}_t$ as a tensor factor model with a CP low-rank structure:
where $r$ denotes the fixed number of factors and $\widetilde a_{ik}$ denotes the $d_k$-dimensional loading vector, which need not be orthogonal. Without loss of generality and to ensure identifiability, we assume that $\mathbb{E} f_{it}^2 =1$ and normalize the factor loadings $\widetilde a_{ik}$ so that $\| a_{ik}\|_2=1$, for all $1\le i\le r$ and $1\le k\le K$. Consequently, all factor strengths are captured by $s_i$. In the strong factor model case, $\|\widetilde a_{ik}\|_2 \asymp \sqrt{d_k}$, which implies that $s_i \asymp \sqrt{d_1d_2\cdots d_K}$. The construction of $s_i$ is a matter of parametrization, which ensures the order of the estimated factor $\widehat f_t$ to be $O_p(1)$ by convention\footnote{Incorporating $s_i$ into the loadings or factors does not improve the convergence speed for the asymptotic normality of the estimated factors discussed in Section 3.}. The noise tensor $\cal{E}_t$ is assumed to be uncorrelated with the latent factors but may exhibit weak correlations across different dimensions. Unlike classical vector factor models, which suffer from rotation ambiguity, the CP tensor factors are uniquely identified up to permutations and sign changes Kruskal1977,Kruskal1989,Bro2000. Throughout this paper, we assume the sign of factors is known without loss of generality.
To construct forecasts for $y_{T+h}$, the CP factor model (ref) needs to be estimated first. We adopt the CC-ISO method proposed by infCP2024 in our context. Specifically, we estimate $f_t$ via the following algorithm.
Step 1. Obtain the initial value $\widehat A_k^{(0)}=(\widehat a_{1k}^{(0)},\ldots,\widehat a_{rk}^{(0)})\in\mathbb{R}^{d_k\times r}$ via randomized composite PCA infCP2024 or tensor PCA BabiiGhyselsPan2025 and compute $\widehat B_k^{(0)} = \widehat A_k^{(0)}(\widehat A_k^{(0)\top}\widehat A_k^{(0)})^{-1} = (\widehat b_{1k}^{(0)},\cdots,\widehat b_{rk}^{(0)})$, where $1\le k\le K$.
Step 2. Given the previous estimates $\widehat a_{ik}^{(m-1)}$, where $m$ is the iteration number, calculate
for $t=1,\cdots,T$. Then the updated loading vectors $\widehat a_{ik}^{(m)}$ are obtained as the top eigenvector of the contemporary covariance $\widehat\Sigma( {\cal Z}_{1:T,ik}^{(m)} )=\frac{1}{T}\sum_{t=1}^T {\cal Z}_{t,ik}^{(m)} {\cal Z}_{t,ik}^{(m)\top}$, where $1\le i\le r$ and $1\le k\le K$.
Step 3. Update $\widehat B_k^{(m)} = \widehat A_k^{(m)}(\widehat A_k^{(m)\top}\widehat A_k^{(m)})^{-1} = (\widehat b_{1k}^{(m)},...,\widehat b_{rk}^{(m)})$ with $\widehat A_k^{(m)}=(\widehat a_{1k}^{(m)},\ldots,\widehat a_{rk}^{(m)})$.
Step 4. Repeat Steps 2 and 3 until the maximum number of iterations $M$\footnote{In our simulation, we set the maximum number of iterations to $M=100$, but convergence is typically achieved in fewer than 5 iterations.} is reached or $\max_{1\le i\le r}\max_{1\le k\le K}\| \widehat a_{ik}^{(m)} \widehat a_{ik}^{(m)\top} - \widehat a_{ik}^{(m-1)} \widehat a_{ik}^{(m-1)\top} \|_{2}\le \epsilon$, where the default accuracy is set to $\epsilon = 10^{-5}$.
Step 5. Obtain the estimated signal as $\widehat s_i = \sqrt{ \frac{\sum_{t=1}^T \left({\cal X}_t\times_{k=1}^K \widehat b_{ik}^\top \right)^2}{T}}$ and the estimated factors as $\widehat{f}_{it}=\widehat{s}_i^{-1}\left({\cal X}_t\times_{k=1}^K \widehat b_{ik}^\top\right)$, for $i=1,\cdots,r$ and $t=1,\cdots,T$.
The estimated factors, $\widehat f_t$, along with $w_t$, are then used to estimate the coefficients in Equation (ref). When the dimension of the non-tensor predictors $w_t$ is small, we estimate (ref) with OLS and the forecast for $y_{T+h}$ is obtained as
$$ \widehat y_{T+h} = \widehat \beta_0^\top w_T + \widehat \beta_1^\top \widehat f_T, $$ where $\widehat \beta_0^\top$ and $\widehat \beta_1^\top$ are OLS estimates.
The above forecasting procedure assumes that the rank of ${\cal X}_t$ is known. However, we need to estimate it in practice. We adopt the contemporary covariance-based unfolded eigenvalue ratio estimator considered in infCP2024. Other estimators, such as the inner-product-based eigenvalue ratio estimator and autocovariance-based eigenvalue ratio estimator, work as well. More details can be found in han2024cp and infCP2024.
In this section, we consider the asymptotic properties of our estimation when the number of non-tensor predictors is relatively small ($p\ll T$) and thus no regularization is required. infCP2024 propose the CC-ISO method and focus on the estimation and inference of loadings while the asymptotic properties of latent factors are unknown. Hence, we first fill in the gap by presenting the consistency and asymptotic normality of the estimated latent factors in Section (ref). Then we derive the inferential theories for the diffusion index model in Section (ref). Section (ref) introduces a robust covariance matrix estimator of the factor process for conducting inference on the conditional mean forecasts.
We start with some assumptions that are necessary for our theoretical development.
Assumption (ref) (i) assumes that the tails of the error and factor processes exhibit exponential decay, which includes a sub-Gaussian distribution as an important example. This assumption could be extended to account for polynomial-type tails with bounded moment conditions. Unlike lam2012, han2024cp and infCP2024, Assumption (ref)(ii) and (iii) allow both weak cross-sectional and serial correlations in the error term. Assumption (ref)(ii) assumes the $\alpha$-mixing property on the factor process, a standard assumption assumed in the tensor factor literature to capture temporal dependence (e.g., chen2021statistical and han2024cp). We acknowledge that the $\alpha$-mixing condition might not be flexible enough to accommodate certain time series models (Andrews1984). Nevertheless, to maintain focus on the essential theoretical developments and ensure analytical tractability, we adopt the $\alpha$-mixing framework in the main analysis. Possible relaxations of this assumption are discussed in Appendix E.
A sufficient condition for Assumption (ref) is $\max_{j} \sum_{l = 1}^d |{\mathbb{E}} \left[ e_{jt} e_{lt} \right]|<\infty$, which ensures that the aggregate dependence across all pairs of cross-sectional units remains bounded as $d$ increases. This condition is mild and commonly used in large-dimensional factor, panel, and matrix-valued time series models (e.g., Bai2003, chen2021statistical). If the cross-sectional dimension has some natural ordering (e.g., spatial or social network data), ${e_{jt}}$ may be assumed to be $\alpha$-mixing in the cross-sectional dimension as well. Namely, for each $t=1,\cdots,T$, ${e_{jt}}$ is $\alpha$-mixing with mixing coefficients $\alpha_t(m)$ such that $\sup_t\alpha_t(m)\leq\alpha(m)$. Then it is straightforward to verify that Assumption (ref) (iii) holds by the mixing inequality. Alternatively, if we take into account the tensor structure, we can consider an example as in Appendix C, which allows for exponentially decaying error correlation along both tensor modes. If there is no natural ordering for cross-sectional indices, one can follow CHEN201271 by introducing a "distance function" between cross-sectional units to define a weak dependence structure that also satisfies Assumption (ref) (iii).
Assumption (ref) (iv) imposes independence between factors and errors, which simplifies the analysis of both the ISO algorithm and the forecast model. A more general assumption allowing for limited dependence, as suggested by Bai2003, could also be considered, though it would introduce significantly greater theoretical complexity.
Unlike Bai2003 and Fanfactorandsparse, Assumption (ref) allows for varying factor strengths by incorporating a mix of strong and weak factors, with certain conditions on the weakest signal strength, the dimensions of tensor data, and the sample size. Specifically, it ensures that, as $d, T \to \infty$, $\max_{i \leq r, k \leq K} \| \widehat a_{ik} \widehat a_{ik}^\top - a_{ik} a_{ik}^\top \|_2 \to 0$, thereby guaranteeing the consistency of the factor estimation.
Let $\psi_0$ denote the estimation error of the warm-start initial estimates for the factor loading vectors. Define
For ease of notation, we define
which represents the final estimation error for the factor loading vectors. We first present the performance bounds of $\widehat{f}_t$ below.
Theorem (ref) shows that $\widehat{f}_t$ is a consistent estimator of the latent factor $f_t$. However, the convergence rate of $\widehat{f}_t$ to $f_t$ may be slower compared to its convergence to $H f_t$. This discrepancy arises due to the non-negligible estimation error associated with the factor signal $s_i$. Nevertheless, this does not affect the prediction of $y_{t+h}$, as the impact is absorbed by the coefficient $\beta_1$. We will provide further discussion on this point in Section (ref). When all factors are strong, i.e., $\alpha_i = 1$ for all $1 \leq i \leq r$, Theorem \ref*{thm:factor1} implies the following:
If we further assume that the error term is serially uncorrelated and follows a sub-Gaussian distribution, then the rates simplify to:
Denote $B = (b_1, b_2, \ldots, b_r) \in \mathbb{R}^{d \times r} $ where $b_i$ = $b_{iK} \odot b_{iK-1} \odot \cdots \odot b_{i1}$ with $b_{ik}$ defined as $B_k = A_k(A_k^{\top} A_k)^{-1} = (b_{1k},...,b_{rk}) \in\mathbb{R}^{d_k\times r}$, $A_k=(a_{1k},\ldots,a_{rk})\in \mathbb{R}^{d_k\times r}$. And denote $\widehat S = \operatorname{diag}(\widehat s_1,\cdots, \widehat{s}_r)$.
Theorem (ref) establishes the asymptotic normality of the estimated factors, confirming that the normal approximation is valid in this context. This result is consistent with the findings of Bai2003 for vector factor models. Additionally, Theorem (ref) derives the asymptotic variance of $\widehat{f}_T$, which provides a theoretical foundation for inference in the diffusion index model (ref) (or (ref)) discussed below. The scaling matrix $H$ does not affect such inference, as it only involves the inner product $\beta_1^\top f_t$, and $\beta_1^\top f_t = \beta_1^\top H^{-1} H f_t$ for any invertible matrix $H$. Thus, the inference remains valid irrespective of $H$. Theorems (ref) and (ref) complement our earlier results in infCP2024, which focus on estimation and inference of loadings.
We first consider the properties of the OLS estimates when the CC-ISO estimates of the latent factors are used as regressors, and then discuss how to construct a confidence interval for the conditional mean of ((ref)).
To take advantage of the faster convergence rate of $\widehat{f}_t$ to $Hf_t$, we rewrite the diffusion index model (ref) as
where $\widetilde \beta_1 = H^{-1} \beta_1$. The conditional mean of $y_{T+h}$ given the information available at time $T$ is
which is an infeasible predictor since it involves the unknown parameters $\beta_0$, $\widetilde \beta_1$ and latent factors $f_T$.
To obtain a feasible forecast, the factor process $f_t$ is first estimated using the CC-ISO algorithm discussed in Section (ref). Then the coefficients $ \widetilde \beta = (\beta_0^\top, \widetilde \beta_1^\top)^\top$ are estimated via OLS:
where $\widehat z_t = (w_t^\top, \widehat f_t^\top)^\top$ and the feasible prediction of $y_{T+h|T}$ is then given by
Denote $z_t = (w_t^\top, f_t^\top)^\top$. To study the asymptotic normality of the OLS estimator $\widehat \beta$, we impose the following assumptions.
These assumptions are standard in both factor and regression analysis. Assumption (ref) (ii) is weaker than the common assumption that regressors and errors are sub-Gaussian with $\nu_3 = \nu_4 = 2$ (see, for example, Fanfactorandsparse, Huang2022, Gao2024). Given this weaker condition, Assumption (ref) (vi) imposes additional conditions on the dimensionality and strength of the signals to ensure the consistency and asymptotic normality of $\widehat f_t$, $\widehat \beta$, and $\widehat y_{T+h|T}$. In particular, it assumes that $d^{\alpha_r}$ grows faster than $d_{\max}$. In the case where $K=2$ and $d_1 \asymp d_2$ such that $d_{\max} \asymp d^{1/2}$, $\alpha_r$ is assumed to be larger than $1/2$, which is also imposed by Bai2023. In the simulation section, however, we demonstrate that the results in the following theorem are robust to the setting where $\alpha_r < 1/2$ when $T$ is large enough. While Assumption (ref) (ii) could be further relaxed to require only bounded fourth moments for errors and regressors, as in BaiNg2006, doing so would necessitate more complex restrictions on dimension and signal strengths. Assumption (ref) (iv) imposes a martingale difference condition on the errors, following BaiNg2006. This assumption could be relaxed to allow for serial correlation at the cost of estimating the long-run variance. To simplify the analysis and maintain interpretability, we maintain the current assumption framework.
Theorem (ref) shows the asymptotic normality of $\widehat \beta$, centered by $\widetilde \beta$, the scaled true coefficient. This result does not hold for the unscaled true coefficient $\beta= (\beta_0^\top, \beta_1^\top)^\top$ because the estimation error of $\widehat f_t$ with respect to $f_t$ is of order $\sqrt{T}$. Nonetheless, it does not affect the inference for the prediction $\widehat y_{T+h|T}$ as shown below. A consistent estimator of the asymptotic variance of $\widehat{\beta}$ can be obtained by the sample covariance matrix of the residuals:
Under conditional homoskedasticity such that ${\mathbb{E}} \left[ \epsilon_{t+h}^2|z_t \right] = \sigma^2_\epsilon$, Equation ((ref)) can be simplified to
where $\widehat \sigma^2_\epsilon = \frac{1}{T} \sum_{t=1}^{T-h} \widehat \epsilon_{t+h}^2$.
The convergence is understood as conditional on $z_T$, which enters only the forecast evaluation but not the estimation of $\widehat{\beta}$. Specifically, given data $\{y_t, z_t\}_{t=1}^T$, our goal is to forecast $y_{T+h|T}$ for a fixed $h$. The coefficient $\beta_0$ is estimated using $\{ y_{t+h}, z_t \}_{t=1}^{T-h}$, since the future observations ${y_{t}:t>T}$ are unavailable.
The two terms in the asymptotic variance of $\widehat y_{T+h|T}$ decay at different rates, so the convergence rate of $\widehat y_{T+h|T}$ is $d^{-\alpha_r/2}+T^{-1/2}$, which implies the efficiency improves with the increase of both the number of observations $T$ and the dimension of the tensor for factor estimation.
Given consistent estimators of $\operatorname*{A\!var}(\widehat{\beta})$ and $\operatorname*{A\!var}(\widehat f_T)$, the prediction interval for $y_{T+h|T}$ with confidence level $\alpha$ can be constructed as
where $q_{1-\alpha/2}$ is the $1-\alpha/2$ quantile of the standard normal distribution, and $$ \widehat \sigma_{y_{T+h|T}}^2 = \frac{1}{T} \widehat z_T^\top \widehat{\operatorname*{A\!var}(\widehat \beta)} \widehat z_T + \widehat \beta_1^\top \widehat S^{-1} \widehat{\operatorname*{A\!var}(\widehat f_T)} \widehat S^{-1} \widehat \beta_1. $$ With $\widehat{\operatorname*{A\!var}(\widehat \beta)}$ given in Equation ((ref)), a consistent estimator of $\operatorname*{A\!var}(\widehat f_T)$ is still needed. Assuming the components of ${\cal E}_t$ are cross-sectionally independent, such that $\Sigma_e = \operatorname{diag}(\sigma_1^2, \ldots, \sigma_d^2)$, $\operatorname*{A\!var}(\widehat f_T)$ can be consistently estimated by
where $\widehat B$ is the estimated $B$ defined on page 10, with the CC-ISO estimator $\widehat a_{jk}$ replacing the unknown $a_{jk}$, $\widehat e_{t}= \mbox{vec}(\widehat {\cal E}_t)$ and $\widehat {\cal E}_t={\cal X}_t- \sum_{i=1}^r \widehat s_i \widehat f_{it} ( \widehat a_{i1} \otimes \widehat a_{i2} \otimes \cdots \otimes \widehat a_{iK})$. If cross-sectional dependence is allowed, a robust variance estimator will be introduced in the next section.
In the context of vector factor models, BaiNg2006 propose the cross-sectional HAC-type estimator of $\operatorname*{A\!var}(\widehat f_T)$ robust to cross-sectional correlation as $$ \widehat{\operatorname*{A\!var}(\widehat f_T)} = \frac{1}{n}\sum_{j=1}^n \sum_{l=1}^n \widehat \Lambda_{j} \widehat \Lambda_{l}^\top \frac{1}{T} \sum_{t=1}^{T-h} \widehat e_{jt} \widehat e_{lt}, $$ where $n$ diverges at a slower rate than $\min\{d, T\}$, and $\widehat \Lambda_j$ denotes the estimated factor loading. This estimator could be extended to the CP factor model by replacing $\lambda_{j}$ with $\widetilde B_{j.}$. However, it is well documented that HAC-type long-run variance estimators often exhibit poor finite-sample performance, particularly when the cross-sectional dimension is large relative to $T$ (see DENHAAN1997; Kiefer2000). Our simulation study in Appendix C confirms this finding in the tensor setting, where the HAC-type estimator tends to produce unreliable variance estimates.
To obtain a more reliable estimator in high-dimensional settings, we adopt a regularized covariance estimation approach that directly targets the structure of $\Sigma_e$. Specifically, we estimate $\Sigma_e$ via a thresholded sample covariance matrix, which shrinks small off-diagonal elements toward zero and yields a more stable and high-dimensionality-robust estimator. This regularization approach replaces the kernel-based smoothing of HAC estimators with an elementwise shrinkage scheme that adapts to approximate sparsity in the error covariance structure.
Recall from Theorem (ref) that the asymptotic variance of $\widehat{f}_t$ is given by $$ \operatorname*{A\!var}(\widehat f_T)=\Sigma_{Be} = B^\top \Sigma_e B, $$ where $B$ can be consistently estimated using the CC-ISO estimator $\widehat B = (\widehat b_1, \ldots, \widehat b_r)$, as shown in infCP2024. The primary challenge lies in estimating the high-dimensional covariance matrix $\Sigma_e$ in the presence of cross-sectional dependence. To address this, we propose a thresholding estimator $\widehat \Sigma_e^{{\cal T}}$: $$ \widehat \Sigma_e^{{\cal T}} = {\cal T}_\lambda \left(\frac{1}{T} \sum_{t=1}^T \widehat e_t \widehat e_t^\top \right), $$ where $\left(\widehat \Sigma_e^{\cal T}\right)_{(j,l)} = {\cal T}_\lambda \left(\frac{1}{T} \sum_{t=1}^T \widehat e_{jt} \widehat e_{lt} \right)$, ${\cal T}_\lambda (\cdot)$ is a thresholding operator and $\widehat e_t$ is the vectorized estimated error using the CC-ISO algorithm. Following Rothman2009, the thresholding operator ${\cal T}_\lambda(\cdot)$ is defined to satisfy the following conditions:
Examples of generalized thresholding include the LASSO penalty rule: $$ {\cal T}_{\lambda}(z)=\operatorname{sgn}(z)\left(|z|-\lambda\right)_{+} $$ and the SCAD thresholding rule proposed by Fan2001scad: $$ {\cal T}_\lambda(z) =
$$ The bound for the estimation error of $\widehat \Sigma_e^{{\cal T}}$ is established uniformly over a class of covariance matrices, as introduced by Bickel2008 and Rothman2009:
for $0 \leq q < 1$. When $q = 0$, this class represents exact sparse covariance matrices, where the number of non-zero entries per column is bounded by $c_0(d)$. For $q > 0$, this class defines approximately sparse covariance matrices, where most of the entries in each column are small. Additional assumptions are imposed to derive the bound for the estimation error of $\widehat \Sigma_e^{{\cal T}}$.
Let $\widetilde a_{ik,j}$ be the $j^{th}$ entry of $\widetilde a_{ik}$ where $\widetilde a_{ik} = d_k^{\alpha_i/2} a_{ik}$.
Assumption (ref) (i) bounds the maximum entry of the factor loadings in model (ref). Similar conditions are used in the strong factor model literature such as Bai2003 and Fan2013. In strong factor models, Assumption (ref) (i) ensures that the factor loadings for each mode are “dense”, i.e., the number of zero entries in each column of $\widetilde A=(\widetilde a_1,...,\widetilde a_r)$, $\widetilde a_i$ = $\mbox{vec}(\widetilde a_{iK} \odot \widetilde a_{iK-1} \odot \cdots \odot \widetilde a_{i1})$, does not increase with $d$. In weaker factor models, however, this number is allowed to increase in $d$ with the rate depending on the factor strength $s_i$. Assumption (ref) (ii) is imposed to ensure that the bound of $\left|e_{it}e_{jt} - {\mathbb{E}} \left[ e_{it}e_{jt} \right]\right|$ is the same as in Bickel2008 and Rothman2009 to accommodate stationary and ergodic errors. This assumption is also imposed in Fan2011 and Fan2013.
The following theorem provides the rate of convergence for $\widehat \Sigma_e^{{\cal T}}$ over the class ${\cal U}(q,c_0(d),M)$.
Fan2013 show that the thresholding estimator $\widehat \Sigma_e^{\cal T}$ with the adaptive thresholding method developed by Cai2011 achieves the same rate as in Theorem (ref) within the strong vector factor model framework. While these results could, in principle, be extended to the CP factor model, the adaptive threshold method presents significant computational challenges when applied to tensor data. Specifically, it requires estimating $\operatorname{var}(e_{jt}e_{lt})$ for all $j$ and $l$, which substantially increases the computational cost due to the high dimensionality of the tensor data. In addition, the adaptive thresholding approach allows $\Sigma_e$ to have diverging diagonal entries, whereas in the CP factor model, the spectral norm of $\Sigma_e$ is typically assumed to be bounded infCP2024,han2024cp. This boundedness assumption aligns with both the theoretical framework and practical considerations of the CP factor model, making the results in Theorem (ref) sufficient for inference in the diffusion index model.
Define $\widehat \Gamma_2 = \widehat B^\top \widehat \Sigma_e^{\cal T} \widehat B$. Theorem (ref) implies the consistency of $\widehat \Gamma_2$, as summarized below.
Corollary (ref) guarantees a valid prediction interval for $y_{T+h|T}$ that remains robust in the presence of potential cross-sectional error correlations.
While diffusion index forecasting with OLS is effective when the number of predictors is relatively small, some real-world applications might involve a large number of potential predictors, sometimes exceeding the sample size. This high-dimensional setting arises in macroeconomic forecasting, financial modeling and trade analysis, where policymakers and researchers need to integrate information from multiple sources. In such contexts, OLS estimation might become unreliable. Moreover, some predictors may be irrelevant, introducing noise rather than improving forecast accuracy. Therefore, it is important to employ variable selection techniques that identify the most relevant predictors while preserving the predictive power of the model. In this section, we extend diffusion index forecasting to accommodate high-dimensional predictors by incorporating regularization---specifically, Multi-Source Factor-Augmented Sparse Regression (MS-FASR)---to ensure robust estimation and improved out-of-sample performance.
Let $w_t \in \mathbb{R}^p$ denote the set of high-dimensional predictors, alongside the tensor time series ${\cal X}_t$. We consider the diffusion index forecast model:
where $p$ is allowed to diverge with the sample size $T$.
Substituting Equation ((ref)) into Equation ((ref)), we obtain:
where $\beta_1^* = \Lambda^\top \beta_0 + \beta_1$. After estimating the factors $f_t$ and $V_t$ from Equation (ref) and (ref), we obtain the estimators of the unknown parameters $\beta_0$ and $\beta_1^*$ via the following penalized regression:
where $\lambda > 0$ is a tuning parameter. Since $\widehat V_t$ is orthogonal to $\widehat f_t$ by construction, the solution to the penalized regression can be obtained via the following steps:
Step 1. Obtain $\widehat f_t$ using the CC-ISO algorithm described in Section (ref).
Step 2. Estimate $\Lambda$ and $V_t$ via OLS:
Step 3. Obtain the projection residuals $\widetilde y_{t+h}$ by regressing $y_{t+h}$ on $\widehat f_t$:
Step 4. Estimate $\beta_0$ by regressing $\widetilde y_{t+h}$ on $\widehat V_t$ using LASSO: $$ \widehat \beta_0 = \operatorname{argmin}_{\beta_0} \frac{1}{2T} \|\widetilde Y - \widehat V \beta_0\|_2^2 + \lambda \|\beta_0\|_1, $$
where $\widehat V = (\widehat V_1, \ldots, \widehat V_{T-h})^\top \in \mathbb{R}^{(T-h) \times p}$ and $\widetilde Y = (\widetilde y_{1+h}, \ldots, \widetilde y_{T}) \in \mathbb{R}^{T-h}$.
Step 5. Estimate $\beta_1$ by $$ \widehat \beta_1 = \widehat \beta_1^* - \widehat \Lambda \widehat \beta_0, $$ and forecast the conditional mean $y_{T+h | T}$ by $$ \widehat y_{T+h|T} := \widehat \beta_0^{\top} \widehat V_T + \widehat \beta_1^{*\top} \widehat f_T. $$
The algorithm is based on residual-on-residual regression, so $V_t$ in Equation ((ref)) should be interpreted as a projection error, rather than the true error from a structural equation. That is, Equation ((ref)) does not necessarily represent the true data generating process (DGP); $w_t$ may have a nonlinear relationship or no relationship with $f_t$. This formulation simplifies theoretical analysis.
For $\varsigma \geq 0$, define the sparsity index set ${\cal S}_\varsigma := \left\{j: | \beta_{0,j} | > \varsigma \right\}$. Let $p_0 := \left| {\cal S}_0 \right|$ denote the cardinality of the support set of $\beta_0$. The following additional assumptions are imposed.
These assumptions are standard in the analysis of high-dimensional regressions. Assumption (ref)(i) is weaker than the common assumption that regressors and errors are sub-Gaussian, as seen in the high-dimensional regression literature (e.g., Loh2012 and Fanfactorandsparse). Assumption (ref)(iii) imposes a compatibility condition, which is less restrictive than directly assuming the positive definiteness of the sample or population covariance matrix. Since $V_t$ is not directly observable in the data, it is more natural to impose the compatibility condition on the population covariance matrix rather than its sample counterpart, as is often done in the high-dimensional regression literature. This approach is also adopted in Adamek2023.
Theorem (ref) shows that diffusion index forecasting remains consistent even in the presence of a large number of potential predictors. The convergence rate of $\widehat\beta_0$ equals the usual LASSO rate plus an additional component associated with factor estimation, while the rate of $\widehat\beta_1$ depends on the estimation error of $\widehat\beta_0$.\footnote{In a standard linear regression estimated by OLS, the Frisch--Waugh--Lovell (FWL) theorem implies that the estimation of $\beta_1$ is unaffected by the estimation of $\beta_0$. However, under the $\ell_1$-penalized framework, the orthogonality doesn't hold. Because the LASSO penalty applies to $\beta_0$, the shrinkage changes the fitted residuals that determine $\widehat{\beta}_1$, and therefore the numerical value and convergence rate of $\widehat{\beta}_1$ depend on the estimation error of $\widehat{\beta}_0$. This feature has been well documented in the literature (see, e.g., chernozhukov2018; Fanfactorandsparse).} The rate condition on $p$ is imposed to simplify the consistency result. Furthermore, by assuming $d^{\alpha_r} \psi^2 = o(1)$, the result can be improved by eliminating the $\psi$ term in the rates. While selection consistency of the penalized regression could be established with much more involved theoretical derivations and additional assumptions, our primary focus is on prediction. Therefore, we leave this extension for future research to maintain clarity and focus.
In this section, we examine the finite-sample properties of the proposed estimators through a simulation study. We consider the following DGP for ${\cal X}_t$ with $r=3$ and $K=2$:
where $u_{it}$ and entries of $Z_t$ are generated independently from ${\cal N}(0, 1)$. Throughout the section, we let $d_1 = d_2$ and let $\Sigma_{{\cal E},k} = \text{Toeplitz}(0.5, d_k)$, $k=1,2$, such that the $(j,l)^{th}$ entry of $\Sigma_{{\cal E},k}$ is equal to $0.5^{|j-l|}$. Factor loadings $A_k = \left(a_{1k},\ldots,a_{rk}\right)$ are generated as follows: let $\widetilde A_k^{({\cal N})} \in \mathbb{R}^{d_k \times r}$ whose elements are generated independently from ${\cal N}(0,1)$. We first generate $\widetilde A_k$ by orthonormalizing $\widetilde A_k^{({\cal N})}$ through QR decomposition, i.e., $\widetilde A_k = \left(\widetilde a_{1k}, \ldots, \widetilde a_{rk} \right) = \operatorname{QR}(\widetilde A_k^{({\cal N})})$. Then $A_k = (a_{1k}, \ldots, a_{rk})$ is generated by $a_{ik} = \Sigma_{{\cal E},k}^{1/2} \widetilde a_{ik} / \sqrt{\widetilde a_{ik}^\top \Sigma_{{\cal E},k} \widetilde a_{ik}}$. We set the factor strength $s_i = (r-i+1) \sqrt{d^{\alpha}}$ with $\alpha \in \left\{0.6,0.4 \right\}$.
In Section (ref) and (ref), we evaluate the consistency and asymptotic distribution of factor estimators. Section (ref) compares the coverage rates of the prediction intervals by CC-ISO and by PCA. Section (ref) illustrates the convergence rates of the LASSO estimators and associated predictions studied in Section (ref). Additional simulation results, including settings with correlated and persistent factors, stronger error dependence, and heavy-tailed (Student-t) disturbances, are provided in Appendix G. Across all designs, the proposed method maintains strong predictive performance and estimation accuracy, confirming its robustness.
In this section, we evaluate the finite-sample performance of the factor estimator $\widehat f_{t}$. Estimation errors are measured as $\|\widehat f_t - Hf_t \|_2$ at $t = T$ where $H$ is defined in Equation (ref)\footnote{Since our primary interest is in forecasting, we report results for $t=T$. Figures for $t=\frac{T}{2}$ show a similar pattern.}. We vary $d_k$ in $\{20,40,60,80\}$ and $T$ in $\{300,400,500\}$.
Figure (ref) presents boxplots of log estimation errors over 1000 repetitions. In all settings, estimation errors decrease as $d_k$ increases. Additionally, estimation errors decrease as factors are stronger. These findings align with Theorem (ref).
Next, we conduct simulations to assess the asymptotic normality of $\widehat f_{t}$, as stated in Theorem (ref), and to evaluate the proposed covariance matrix estimator in Theorem (ref). We vary $d_k$ in $\{40,60,80\}$ and let $T = 800 + \lceil d^{3/4} \rceil$.
Specifically, we use the SCAD thresholding function developed by Fan2001scad, defined as $$ {\cal T}_\lambda(z) =
$$ where we set $a = 3.7$ as suggested in \cite{Fan2001scad}.\footnote{ The other three thresholding functions (hard thresholding, soft thresholding and adaptive LASSO) considered in \cite{Rothman2009} are also evaluated, yielding similar simulation results.} The threshold $\lambda$ is set to $\sqrt{\log(d) / T} + \sqrt{1/d}$.
Figure (ref) shows the distribution of $\widehat \Sigma_{Be}^{-1/2}\widehat S(\widehat f_T - H f_T)$ over 2000 repetitions under two factor strengths. We note that the distribution of $\widehat f_T$ approximates the standard normal distribution, which validates Theorem (ref). Furthermore, the result remains robust to cross-sectional dependence, supporting the effectiveness of the proposed thresholding covariance matrix estimation.
In this section, we examine the prediction intervals for $y_{T+1|T}$ constructed based on Theorem (ref). The target variable $y_{t+1}$ is generated as $$ y_{t+1} = \beta_0 + \beta_1^\top f_t + \epsilon_{t+1}, $$ where $\beta_0 = 0.5$ and $\beta_1 = \left(0.5,0.5,0.5\right)$. The idiosyncratic error $\epsilon_{t+1}$ is drawn independently from ${\cal N}(0,\nu_t)$ with $\nu_t$ drawn independently from $U[0.5,1.5]$.
Set $d_1 = d_2 \in \{20,40,60,80,120,160\}$ and $T = 800 + \lceil d^{3/4}\rceil$, with a confidence level of 0.95. We assess the finite-sample performance of the prediction interval for $\widehat y_{T+1|T}$ proposed in Equation ((ref)) and compare it with the vector PCA method of BaiNg2006. For this comparison, we apply the classical PCA method to the vectorized tensor $x_t := \Vec({\cal X}_t) \in \mathbb{R}^{d}$ and construct the confidence interval following BaiNg2006 and Bai2023:
where $\widehat \sigma^2_{y_{T+h|T},pca} = \frac{1}{T} \widehat z_T^{(pca)\top} \ \widehat \operatorname*{A\!var}(\widehat \beta^{(pca)}) \widehat z_T^{(pca)} + \frac{1}{d} \widehat \beta_1^{(pca)\top} \widehat \operatorname*{A\!var}(\widehat f_T^{(pca)}) \widehat \beta_1^{(pca)}$. The variance estimator for $\widehat f_T^{(pca)}$ is given by $$ \widehat \operatorname*{A\!var}(\widehat f_T^{(pca)}) = \Tilde V^{-1} \widehat \Gamma_t \Tilde V^{-1}, $$ where $\Tilde V$ is a diagonal matrix with diagonal elements equal to the top $r$ eigenvalues of $\frac{1}{dT} \sum_{t=1}^T x_t x_t^\top$, and $\widehat z_T^{(pca)}$, $\widehat \beta^{(pca)}$, and $\widehat f_T^{(pca)}$ are the corresponding PCA estimators. We consider two types of $\widehat \Gamma_t$. The first one is the $\widehat A^{(PCA)^\top}\widehat \Sigma^{({\cal T})}_{e,pca} \widehat A^{(PCA)}$ where $\widehat A^{(PCA)}$ are factor loadings estimated via PCA and $\widehat \Sigma_{e,pca}^{({\cal T})}$ is the proposed thresholding estimator of the covariance matrix of error terms for PCA. The second one is the HAC-type estimator proposed by BaiNg2006 and Bai2023: $$ \widehat \Gamma_t^{(HAC)} = \frac{1}{n} \sum_{j = 1}^n \sum_{l = 1}^n \widehat A_{j:}^{(\text{PCA})} \widehat A_{l:}^{\top(\text{PCA})} \frac{1}{T} \sum_{t=1}^T \widehat e_{jt}^{(\text{PCA})} \widehat e_{lt}^{(\text{PCA})}, $$
where $\widehat A_{j:}^{(\text{PCA})}$ and $\widehat e_{jt}^{(\text{PCA})}$ are factor loadings and errors estimated via PCA, respectively. The tuning parameter is set as $n = \min\{ \sqrt{d}, \sqrt{T}\}$ as suggested by BaiNg2006. For both CP and PCA approaches, $\operatorname*{A\!var}(\widehat{\beta})$ is estimated using Equation (ref).
Table (ref) shows the coverage rates of three estimated prediction intervals under two different values of $\alpha$, with a confidence level $95\%$. For $\alpha=0.6$, the coverage rates for the CP-based approach are close to the nominal level. For $\alpha = 0.4$, the coverage rate is slightly lower when $d_k = 20$ but converges to the nominal level as $d_k$ increases. In contrast, the PCA-based approach fails to produce reliable prediction intervals: its coverage rates deviate significantly from the nominal level and show no improvement with increasing $d_k$.
In this section, we evaluate the convergence rate of $\widehat \beta_0$ and $\widehat y_{T+1|T}$ in Theorem (ref). Consider the following DGP for $y_{t+h}$ and $w_t \in \mathbb{R}^p$:
where $z_t = \left(1, f_t^\top\right) \in \mathbb{R}^{r+1}$. We set the predictor dimension to $p=200$, with the first three elements of $\beta_0$ equal to $0.5$ and the remaining elements set to 0. Each entry of $\Lambda$ is drawn from the uniform distribution $U[-1,1]$, and the entries of $V_t$ are generated independently from $N(0,1)$. The idiosyncratic errors $\epsilon_{t+h}$ follow the same setting as in Section (ref). We fix $d_1 = d_2 = 40$ and vary $T$.
In this setting, the rate of $\| \widehat \beta_0 - \beta_0 \|_1$ is bounded above by $p_0 \sqrt{\log(p) / T}$, while the forecast error $\left|\widehat y_{T+h|T} - y_{T+h|T}\right|$ is bounded above by $p_0 \log(p) / \sqrt{T}$, given $\nu_5=2$ for Gaussian $V_T$. We choose $T$ such that $p_0 \sqrt{\log(p) / T}$ takes values on a uniform grid in $[0.15,0.5]$, which implies that $p_0 \log(p) / \sqrt{T}$ ranges in $(0.34,1.15)$. The tuning parameter for the LASSO regression is fixed at $\sqrt{\log(d)/T} + 1/\widehat{s_r}$, where $\widehat{s_r}$ is the estimated weakest factor signal, $s_r$.
Figure (ref) reports the estimation and prediction errors. The results provide further support for the theoretical findings established in Section (ref).
Understanding trade flow patterns and forecasting their dynamics are essential for policymaking, firm optimization, and risk management. Trade data inherently form a dynamic sequence of tensor variates, which can capture network-like structures, underlie common dynamics, and reveal intricate interaction patterns. In this section, we consider diffusion index forecasting based on the CP tensor factor model for international trade data, providing a unified framework to estimate global trade factors and predict future variations in U.S. trade.
We analyze monthly bilateral import and export volumes of commodity goods among 24 countries and regions from January 1999 to December 2018, using data from the International Monetary Fund Direction of Trade Statistics (IMF-DOTS). The countries and regions included in the dataset are: Australia (AU), Canada (CA), China Mainland (CN), Denmark (DK), Finland (FI), France (FR), Germany (DE), Hong Kong (HK), Indonesia (ID), Ireland (IE), Italy (IT), Japan (JP), Korea (KR), Malaysia (MY), Mexico (MX), Netherlands (NL), New Zealand (NZ), Singapore (SG), Spain (ES), Sweden (SE), Taiwan (TW), Thailand (TH), United Kingdom (GB), and the United States (U.S.).
In our study, we employ the diffusion index model with a CP low rank structure, as defined in ((ref)) and ((ref)). Specifically, we represent the trade data as a $24 \times 24$ two-dimensional tensor, where each element $x_{i,j,t}$ denotes the monthly variation of exports from country $i$ to country $j$ at month $t$. For simplicity, self-exports are set to zero, i.e., $x_{i,i,t} = 0$ for all $i$ and $t$. The target variables for our analysis are the monthly variation of U.S. aggregate export and import to/from countries in the sample, denoted by $y_t^{ex}$ and $y_t^{im}$, respectively.
The number of common factors is determined using the eigenvalue-ratio-based method proposed by Ahn2013 and infCP2024, which identifies four common factors explaining 51.1% of the total variance. Let $f_t$ denote the common factors extracted from the growth rate of bilateral trade. We then construct one-month-ahead forecasts for monthly variations in U.S. aggregate exports and imports using the following regression:
Table (ref) reports the in-sample forecasting results based on Equation (ref). As a benchmark, regression (a) predicts each target variable using only the first lag of changes in U.S. exports and imports. In contrast, regression (b) demonstrates that incorporating common factors extracted from the tensor data significantly increases predictive power compared to using lagged values alone. Specifically, the common factors explain $29\%$ of the variation in monthly export changes and 21% of the variation in import changes. Regression (c) integrates both the lagged target variables and the common factors, leading to a substantial improvement in explanatory power, with R-squared values increasing to 40% for U.S. exports and 30% for imports. Moreover, all four common factors are statistically significant predictors for both trade flows. These findings highlight the crucial role of common factors derived from the tensor data in enhancing the accuracy of monthly U.S. export and import forecasts.
Since the factors in the CP tensor model are identified only up to sign changes, it is meaningful to explore their economic interpretation. To characterize these factors, we examine their correlations with monthly variations in bilateral trade flows among the selected countries. These correlations are visualized in the heatmap shown in Figure (ref), where stronger correlations between a factor and bilateral trade flows are indicated by deeper blue shades.
The heatmap reveals distinct regional patterns for each factor. Factor 1 is closely associated with exports from Asian countries to the rest of the dataset, with the highest correlation observed in exports from CN. Factor 2 is highly correlated with China's imports from most countries in the dataset and also shows notable correlations with trade flows among key Asian economies, including CN, KR, JP, and SG. Factor 3 is mainly correlated with bilateral trade flows among European countries, while Factor 4 predominantly captures trade flows within North America, specifically among U.S., CA and MX. In summary, Factors 1 and 2 contain information on trade flows within Asia, particularly involving China. Factor 3 relates to trade dynamics within Europe, and Factor 4 captures variations in trade among North American countries. The in-sample analysis in Section (ref) demonstrates that these factors are not only economically interpretable but also provide significant predictive power for monthly variations in U.S. exports and imports.
In this section, we evaluate the out-of-sample performance of the diffusion index model (ref) based on the CP low-rank structure and compare it with alternative methods, in particular the vector factor model studied by BaiNg2006. In addition to Model (ref), we incorporate 126 macroeconomic variables from FRED-MD McCrackenNg2016 and up to 12 lags of monthly variations of U.S. aggregate exports and imports. This allows us to assess the performance of MS-FASR, introduced in Section (ref) and investigate whether including U.S. macroeconomic variables and additional lags of the target variables improves the out-of-sample forecast of U.S. aggregate export and import variations. The MS-FASR model is specified as follows:
where $w_t \in \mathbb{R}^{148}$ includes 126 macroeconomic variables and lagged U.S. aggregate exports and imports from lag 2 to lag 12\footnote{$y_{t}^{(ex)}$ and $y_{t}^{(im)}$ represent lag 1 export and import variations and are already included in the model.}.
The out-of-sample analysis follows an expanding-window approach, where model parameters are re-estimated as new data become available. The process begins with an initial five-year sample from December 1999 to December 2004. Factors and parameters are estimated using data from December 1999 to November 2004,\footnote{The predictor variables span December 1999 to October 2004, while the target variables cover January 2000 to November 2004, forming a five-year training sample.} and the model is then used to forecast monthly variations in U.S. aggregate exports and imports for December 2004. This procedure is repeated iteratively until the end of the sample, resulting in a total of 169 monthly forecasts from December 2004 to December 2018.
The tuning parameter $\lambda$ is selected via an expanding forecast validation scheme following song2011large and han2015direct, which is appropriate for time-series settings. Specifically, we divide the sample into an initial training subsample $t = 1, \ldots, \lceil \gamma T\rceil$ and a validation sample $t = \lceil \gamma T\rceil + 1, \ldots, T$, with $\gamma=0.8$. For each candidate penalty $\lambda_k$, the model is recursively re-estimated and used to generate one-step-ahead forecasts over the validation period. The value of $\lambda_k$ that minimizes the mean squared prediction error is selected.
We compare the performance of Model (ref) and Model (ref) against various alternative methods\footnote{ To our knowledge, there are no well-established forecasting benchmarks for international trade flows. Some research, such as Bussiere2009 and Greenwood2012, employ Global VAR (GVAR) to capture international linkages. However, GVAR relies on pre-specified weighting matrices and requires a consistent set of macroeconomic indicators across countries at the same frequency, which is infeasible for monthly data.}:
Table (ref) presents the mean squared error (MSE) ratios of the one-month-ahead out-of-sample forecast for each model relative to the benchmark. It also presents p-values from the forecast comparison tests of Diebold1995 (DM). These tests are one-sided, with the following alternatives:
Our findings indicate that, for both exports and imports, MS-FASR(CP) achieves the lowest MSE among all the methods considered. First, MS-FASR(CP) significantly outperforms LASSO(w), reinforcing our in-sample results that common factors extracted from tensor data are valuable for predicting U.S. export and import variations. Furthermore, MS-FASR(CP) outperforms both DI(CP) and DI(PCA) with DM test p-values smaller than any conventional significance level, suggesting that macroeconomic variables provide additional predictive power. Additionally, MS-FASR(CP) also outperforms DI(CP) + DI(w) and DI(PCA) + DI(w) models, which attempt to incorporate factors from multiple sources. This result provides strong empirical support for combining the CP tensor model with sparse regression. Notably, between CP and PCA, DI(CP) significantly outperforms DI(PCA); MS-FASR(CP) modestly improves upon MS-FASR(PCA), consistent with the theoretical argument in Remark (ref).
The superior empirical performance of the MS-FASR method can be attributed to its ability to integrate multiple sources of information in a statistically coherent and efficient way. Specifically, the method jointly exploits (i) low-dimensional factors extracted from the tensor predictor, which capture common cross-country dynamics, and (ii) a high-dimensional set of macroeconomic predictors $w_t$, which provide complementary, country-specific signals. Unlike conventional diffusion-index regressions, MS-FASR selectively penalizes only the coefficients on $w_t$ while keeping factor components unpenalized. This structure preserves systematic global information from the tensor factors while preventing overfitting from noisy or redundant local predictors. Moreover, the residual-on-residual estimation step ensures that the penalized regression operates on information orthogonal to the factor space, mitigating multicollinearity and enhancing out-of-sample stability. Competing models either rely solely on factor information or treat all predictors symmetrically, which can reduce forecasting efficiency when predictive sources are heterogeneous. MS-FASR’s hybrid structure thus allows it to combine global coherence with local adaptability, yielding substantial gains in predictive accuracy.
To quantify the relative contributions of global and local information, we conduct a two-component Shapley attribution (Shapley) that decomposes the total gain in forecast accuracy relative to the benchmark into contributions from local predictors ($w_t$) and global factors ($f_t$)\footnote{ The Shapley decomposition provides an order-invariant and symmetric measure of how much each information source contributes to the overall reduction in MSE relative to the benchmark model. Let $\text{MSE}_{\text{Benchmark}}$, $\text{MSE}_{\text{Lags+f}}$, $\text{MSE}_{\text{Lags+w}}$, and $\text{MSE}_{\text{MS-FASR}}$ denote the MSEs of the benchmark, lags and global factor only, lags and local predictors only, and MS-FASR models, respectively. The Shapley contributions are \[ \phi_f = \tfrac{1}{2}\big[\text{MSE}_{\text{Benchmark}} - \text{MSE}_{\text{Lags+f}}\big] + \tfrac{1}{2}\big[\text{MSE}_{\text{Lags+w}} - \text{MSE}_{\text{MS-FASR}}\big], \] \[ \phi_w = \tfrac{1}{2}\big[\text{MSE}_{\text{Benchmark}} - \text{MSE}_{\text{Lags+w}}\big] + \tfrac{1}{2}\big[\text{MSE}_{\text{Lags+f}} - \text{MSE}_{\text{MS-FASR}}\big], \] with $\phi_f+\phi_w = \text{MSE}_{\text{Benchmark}} - \text{MSE}_{\text{MS-FASR}}$. Each $\phi$ represents the average marginal MSE reduction attributable to that component.}. The result shows that both local and global information contribute meaningfully to MS-FASR’s forecasting gains, with local predictors accounting for a slightly larger share. For exports, $54.5\%$ of the total MSE reduction is attributed to local information and $45.5\%$ to global factors; for imports, the shares are $61.3\%$ and $38.7\%$, respectively. This suggests that while local variation remains somewhat more influential, global factors also provide substantial complementary information, particularly for export forecasting.
Factor models are powerful tools for extracting meaningful information from high-dimensional data, which can then be used for prediction. This paper studies the case where the data naturally take the form of a tensor and can be represented by CP decomposition. We develop inferential theories for factor estimation and predictive intervals in the diffusion index forecasting model. We establish that the least squares estimates from predictive regressions are $\sqrt{T}$-consistent and asymptotically normal, even in the presence of weaker factors. Furthermore, we show that the conditional mean remains consistent and asymptotically normal, with its convergence rate determined by $T$ and the strength of the weakest factor. For predictive inference, we propose a consistent estimator for the high-dimensional covariance matrix of cross-sectionally correlated and heteroskedastic errors.
Additionally, we consider settings where multiple data sources with different structures are available and introduce the MS-FASR model, which effectively integrates information across datasets. Simulation studies confirm our theoretical results, and an empirical application demonstrates that leveraging the tensor structure enhances predictive performance. Our findings suggest that incorporating tensor-based factor extraction can lead to substantial improvements over existing forecasting methods.