Estimation and Inference for CP Tensor Factor Models

Estimation and Inference for CP Tensor Factor Models

\if00 {

\affil[1]{University of Rochester} \affil[2]{University of Notre Dame}

} \fi

\if10 \fi

Introduction

Factor models have become one of the most popular tools for summarizing and extracting information from high-dimensional data in economics and finance (fan2021review, Bai2016review, Stock2016review). Traditional factor models are designed to manage large panel data, where both cross-sectional and time series dimensions increase. These models admit a low-rank structure and have a common-idiosyncratic decomposition, allowing for the identification of significant variations within the panel of economic data.

In modern economics, researchers increasingly encounter vast, multi-dimensional datasets, or tensor. For example, monthly import-export volume time series spanning various product categories among countries can be represented as a three-dimensional tensor, with unavailable diagonal elements for each product category. Similarly, in portfolio selection, data often involve stock prices and various firm characteristics over time across different firms, forming a two-dimensional tensor. Additionally, macroeconomic studies on growth and productivity analyze multiple macro variables at the country-industry level, enabling cross-country comparative analyses, which are challenging with traditional panel data.

Statistical methods and economic applications for high-dimensional tensor factor analysis are still in their early stages of development. As in the classical panel setting, tensor factor models typically assume low-rank structures, with Canonical Polyadic (CP) and Tucker structures being the most common choices (see, e.g., KB2009review). Recent studies have explored various estimation approaches and extensions. For example, working with Tucker decomposition, chen2022factor consider two estimators based on the autocovariance matrices, while Han2021iterative extend these methods using an iterative procedure with the matrix unfolding mechanism. chen2023statistical propose an $\alpha$-PCA method that preserves the matrix structure and aggregates mean and contemporary covariance through a hyper-parameter $\alpha$. chen2024rank introduce a pre-averaging technique for the Tucker tensor factor model that significantly enhances the model's inherent signal strength under certain conditions. In parallel, lettau20243d provides a pioneering empirical application of Tucker decomposition to a three-dimensional panel of double-sorted portfolios, showing substantial improvements over Fama–French factors and traditional PCA. In the context of CP decomposition, han2024cp propose an iterative simultaneous orthogonalization algorithm with warm-start initialization, while babii2022tensor employ tensor principal component analysis (TPCA), assuming orthogonal factor loadings. chang2023modelling develop an estimation procedure based on a generalized eigenanalysis (GE) constructed from the serial dependence structure of the underlying process.

In this paper, we focus on a tensor factor model with a CP low-rank structure due to its economic relevance and parsimonious features. We propose an iterative projection estimation based on the contemporary covariance and develop new inferential theories. Our contributions advance the existing literature in several ways:

First, unlike han2024cp, which rely on lagged autocovariances, we use the contemporary covariance, which captures the full variance structure of the data. This is crucial in applications where serial dependence is weak, such as asset returns, where the efficient market hypothesis implies little autocorrelation. Our method significantly broadens the applicability of CP tensor factor models.

Second, inspired by techniques in anandkumar2014guaranteed and sun2017provable with noisy CP decomposition, we introduce an initialization strategy via randomized composite PCA (RC-PCA), which accommodates closely spaced eigenvalues. While those works consider CP decomposition of a single tensor, we establish performance bounds in the presence of noise and latent factor randomness---substantially complicating the theoretical analysis.

Third, we provide the first asymptotic normality results for estimated loading vectors in CP tensor factor models. This result enables valid statistical inference, such as confidence interval construction and hypothesis testing, tools that are essential for empirical applications\footnote{The importance of developing inferential theory in factor models is well established in the econometric literature; for instance, bai2003 laid the foundation for inference in high-dimensional factor models and has had a lasting impact. Our work brings similar inferential ideas to the tensor setting, filling an important gap in the existing literature.}. Our inferential framework is based on the asymptotic distribution of singular vectors and differs fundamentally from the technical tools used in both vector factor models (bai2003) and Tucker tensor models (chen2023statistical, yu2022projected). Our approach can be extended to inference in many other matrix and tensor problems, such as the inference of factor loading vectors using lagged autocovariance matrices, as considered in lam2012,chang2023modelling.

Fourth, we extend the eigenvalue ratio criterion---previously used in Tucker models (e.g., han2022rank)---to the CP tensor factor setting and establish its consistency. To our knowledge, this is the first work to do so.

Finally, our empirical analysis using characteristic-sorted portfolios demonstrates the practical value of our method. The proposed CP factor model outperforms benchmark approaches in out-of-sample forecasts. The factors extracted from the CP low-rank structure yield smaller cross-sectional pricing errors than classical factor models in the literature. Moreover, the estimated loadings and factors are interpretable without post-hoc rotation, revealing economically meaningful factors: a market factor, a long-short factor, and a volatility factor.

The remaining sections of this paper are organized as follows. Section (ref) introduces the high-dimensional tensor factor model with a CP low rank structure allowing for non-orthogonal loading vectors. In Section (ref), we present an iterative projection estimation procedure and two generalized eigenvalue ratio-based estimators for the number of latent factors. Section (ref) establishes inferential theory. We assess the finite sample performance through simulation in Section (ref) and apply our method to characteristic portfolios in Section (ref). Finally, Section (ref) concludes the paper with all mathematical proofs and additional simulations included in the Appendix.

Notations and preliminaries

In this subsection, we introduce essential notations and basic tensor operations. For an in-depth review, readers may refer to KB2009review.

Let $[n]$ denote the set $\{1, 2, \ldots, n\}$. Let $\|x\|_q = (x_1^q+...+x_p^q)^{1/q}$, $q\ge 1$, for any vector $x=(x_1,...,x_p)^\top$. We employ the following matrix norms: matrix spectral norm $\|M\|_{2} = \underset{\|x\|_2=1,\|y\|_2=1}{\max} \|x^\top M y\|_2 = \sigma_1 (M)$, where $\sigma_1(M)$ is the largest singular value of $M$. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$, we write $a_n\asymp b_n$ if there exists a constant $C$ such that $|a_n|\leq C |b_n|$ and $|a_n|\geq C |b_n|$ hold for all sufficiently large $n$, and $a_n\lesssim b_n$ if there exists a constant $C$ such that $a_n\le Cb_n$.

Consider two tensors ${\cal A}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}, {\cal B}\in \mathbb{R}^{r_1\times r_2\times \cdots \times r_N}$. The tensor product $\otimes$ is defined as ${\cal A}\otimes {\cal B}\in \mathbb{R}^{d_1\times \cdots \times d_K \times r_1\times \cdots \times r_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 $k$-mode product of ${\cal A}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}$ with a matrix $U\in\mathbb{R}^{m_k\times d_k}$ is an order $K$ tensor of dimension $d_1\times \cdots \times d_{k-1} \times m_k\times d_{k+1} \times \cdots \times d_K$, denoted as ${\cal A}\times_k U$, 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}. $$ Similarly, for a matrix $V\in \mathbb{R}^{d_k\times d_\ell}$, define ${\cal A}\times_\ell\times_k V\in\mathbb{R}^{d_1\times\cdots \times d_{\ell-1}\times d_{\ell+1} \times\cdots \times d_{k-1}\times d_{k+1}\times\cdots \times d_K}$ as $$ ({\cal A}\times_\ell\times_k V)_{i_1,...,i_{\ell-1},i_{\ell+1},...,i_{k-1},i_{k+1},...,i_K}=\sum_{i_{\ell}=1}^{d_{\ell}}\sum_{i_k=1}^{d_k} {\cal A}_{i_1,i_2,...,i_K} V_{i_\ell,i_k}.$$ The mode-$k$ matricization of a tensor ${\cal A}\in \mathbb{R}^{d_1\times\cdots\times d_K}$ is denoted as $\hbox{\rm mat}_k({\cal A})\in \mathbb{R}^{d_k\times d_{-k}}$, where $d=\prod_{j=1}^K d_j$ and $d_{-k}=d/d_k=\prod_{j=1,j\neq k}^K d_j$. It is obtained by setting the $k$-th tensor mode as its rows and collapsing all the others into its columns. And the vectorization of the matrix/tensor ${\cal A}$ is denoted as ${\rm{vec}}({\cal A})\in \mathbb{R}^d$. With a slight abuse of notation, we still define $\hbox{\rm mat}_k(\mbox{vec}({\cal A}))=\hbox{\rm mat}_k({\cal A})$. For nonempty $J\subseteq [K]$, $\text{mat}_J({\cal A})$ is the mode $J$ matrix unfolding which maps ${\cal A}$ to $d_J\times d_{-J}$ matrix with $d_J=\prod_{j\in J}d_j$ and $d_{-J}=d/d_J$, e.g. $\text{mat}_{\{1,2\}}({\cal A})=\text{mat}_3^\top({\cal A})$ for $K=3$.

Model

We consider a tensor-valued time series $\mathcal{Y}_t \in \mathbb{R}^{d_1 \times d_2 \times \cdots \times d_K}$, where $ K \ge 2$ and $1 \le t \le T$\footnote{When $K=1$, $\mathcal{Y}_t$ reduces to a vector, and model ((ref)) becomes the classical factor model extensively studied in the literature (bai2002 and stock2002). The identification condition discussed in Remark (ref) does not apply in this case, and hence we assume $K\ge 2$. }. Our focus is on a tensor factor model with a CP low-rank structure:

where $\otimes$ denotes the tensor product, $f_{it}$ is a one-dimensional latent factor, $\Gamma_{ik}$ denotes the $d_k$-dimensional loading vector, which needs not to be orthogonal\footnote{Our framework naturally accommodates orthogonal loadings or factors as a special case, although our identification condition does not require orthogonal loadings or factors KB2009review.}. Without loss of generality and to ensure identifiability, we assume that $\mathbb{E} f_{it}^2 =1$ and normalize the factor loadings $\Gamma_{ik}$ so that $ a_{ik}= \Gamma_{ik} /\|\Gamma_{ik}\|_2$, for all $1\le i\le r$ and $1\le k\le K$. Consequently, all factor strengths are captured by $w_i$ with $w_i=\prod_{k=1}^K \|\Gamma_{ik}\|_2$. In the strong factor model case, $\|\Gamma_{ik}\|_2\asymp d_k^{1/2}$, which implies that $w_i \asymp \sqrt{d_1d_2\cdots d_K}$. Unlike the uncorrelated factors assumed in han2024cp, we allow for moderately strong correlation structures among the individual factors. The noise tensor $\mathcal{E}_t$ is assumed to be uncorrelated with the latent factors but may exhibit weak correlations across different dimensions. The rank $r$ may either be fixed or divergent.

A natural alternative approach to extracting common factors is to vectorize the data:

where $\mbox{vec}({\cal Y}_t) \in \mathbb{R}^{d}$ with $d=d_1d_2\cdots d_K$ and $F_t=\left(f_{1t},f_{2t},\cdots,f_{rt}\right)^\top \in \mathbb{R}^{r}$. However, this method ignores the tensor structure of the data and hence substantially increases the number of parameters in the loading matrices from $(d_1+d_2+\cdots+d_K)r$ in the tensor case to $(d_1d_2\cdots d_K)r$ in the stacked vector version. In contrast, our proposed model preserves the tensor structure by modeling $\mbox{vec}({\cal Y}_t)$ as $A W F_t + \mbox{vec}({\cal E}_t)$, $W=\text{diag}(w_1,...,w_r), A=(a_1,\cdots,a_r)$ and $a_i=\mbox{vec}( a_{i1}\otimes a_{i2}\otimes\cdots\otimes a_{iK})$. This structure not only reduces the parameter dimension but also leads to improved convergence rates due to its parsimonious structure.

Consider an illustrative example of sorted portfolios, detailed in Section (ref). The observed $\mathcal{Y}_t$ is represented as a matrix, where $d_1$ is the number of characteristic-sorted portfolios, $d_2=10$, and $K=2$. Each entry $\mathcal{Y}_{t,jl}$ corresponds to the excess return of the $l^{\text{th}}$-decile of the $j^{\text{th}}$ characteristic at time $t = 1, \ldots, T$. Figure (ref) shows a time series plot of $\mathcal{Y}_t$ for ten characteristics from January 1990 to December 2022. Model ((ref)) identifies $r$ latent factors, which can be interpreted as systematic risk factors. The element $ a_{i1,j}$ of the loading vector $ a_{i1}$, where $i = 1,\cdots,r$ and $j = 1,\cdots,d_1$, captures the heterogeneous exposure of the $j^{\text{th}}$ characteristic to the $i^{\text{th}}$ risk factor. Similarly, the entry $ a_{i2,l}$ of the loading vector $ a_{i2}$, where $i = 1,\cdots,r$ and $l = 1,\cdots,10$, determines the exposure of the $l^{\text{th}}$ decile to the $i^{\text{th}}$ risk factor. We allow the number of risk factors $r$ to increase with the dimensions $d_1$, $d_2$ and the sample size $T$.

In the literature, an alternative tensor factor model based on Tucker decomposition has been explored (see, e.g., Han2021iterative, wang2017tensor, lettau20243d):

where the core tensor ${\cal F}_t\in \mathbb{R}^{r_1\times\cdots \times r_K}$ is the latent factor process in a tensor form, and $ A_i$'s are $d_i\times r_i$ loading matrices. As discussed in babii2022tensor and han2024cp, unlike the CP decomposition, the Tucker decomposition is generally non-unique, leading to identification difficulties: even with the usual identification restrictions, only the column spaces of loading matrices are identified. Consequently, estimation results from model ((ref)) may exhibit ambiguity, undermining meaningful discussions of individual factors\footnote{For classical factor models, stock2002 point out, “because the factors are identified only up to a $k\times k$ matrix, detailed discussion of the individual factors is unwarranted”. The same comment applies to Tucker tensor factor models.}. In contrast, the CP tensor factor model ((ref)) yields $r$ scalar factors, aligning with the conventional wisdom in economic applications regarding common factors. This set of one-dimensional latent factors serves as natural inputs for diffusion index forecasting and factor-augmented regressions\footnote{A prominent example is the modeling of global yield curve dynamics. diebold2008 show that the generalized Nelson-Siegel model accurately captures the dynamics of yield curves and delivers strong predictive performance. Their model assumes that global yields follow a linear factor structure with interpretable loadings. This specification fits naturally within our CP factor model framework, where $f_t$ captures global factors, and $a_{i1}$ and $a_{i2}$ represent country-specific and maturity-specific loadings, respectively.}. We regard the CP tensor factor as a more parsimonious yet flexible and economically relevant alternative. Further comparison of the performance of these two tensor factor models will be presented in Section (ref).

Estimation

We consider a two-step estimation procedure to derive the loading vectors and latent factors. This approach begins with initialization via RC-PCA, followed by an iterative refinement step utilizing an iterative simultaneous orthogonalization (ISO) procedure.

We start by defining the contemporary covariance as the expected value of the outer product of ${\cal Y}_t$:

where $\Theta_{ij}=w_i w_j {\mathbb{E}} \left[ f_{it}f_{jt} \right]$. Its sample analogue, denoted as $\widehat{\Sigma}$, is computed as the average outer product over $T$ observations:

We aim to estimate the loading vectors by minimizing the empirical quadratic loss, formulated as:

where $\|{\cal A}\|_{F}$ denotes the Frobenius norm of a tensor ${\cal A}$. However, this optimization problem is non-convex and prone to multiple local optima. To counter this problem, we employ a two-step approach. The first step focuses on obtaining a suitable initialization close to the global optimum.

The contemporary covariance $\Sigma$ in ((ref)) can be unfolded to a $d\times d$ matrix

where $\Theta= W (\mathbb{E} F_t F_t^\top) W$, $F_t=(f_{1t},\cdots,f_{rt})^\top$, $W=\operatorname{diag}(w_1,...,w_r)$. This unfolding enables classical PCA estimation if the columns of the loading matrix $ A$ are orthogonal. Our framework accommodates general non-orthogonal $a_i$'s and hence the PCA procedure introduces a bias component, which motivates the second stage refinement. The accuracy of the PCA estimator hinges on the maximum correlation among the loading vectors. When the additional orthogonality condition is imposed as in babii2022tensor, the maximum correlation reduces to $0$ and hence bias disappears. The first step, termed initialization via RC-PCA, is detailed in Algorithm (ref).

To further relax the eigengap assumption imposed in babii2022tensor and han2024cp, we incorporate randomized projection (Procedure (ref)) into RC-PCA approach\footnote{Closely spaced eigenvalues do not imply that the covariance matrix is singular. A covariance matrix can remain full rank even when some of its eigenvalues are identical or nearly identical. In our theoretical analysis, “closely spaced eigenvalues” specifically refers to identical or nearly identical eigenvalues of the population covariance matrix of the noiseless data, defined as $\Theta= W (\mathbb{E} F_t F_t^\top) W$, where $F_t=(f_{1t},\cdots,f_{rt})^\top$, $W=\operatorname{diag}(w_1,...,w_r)$. The estimation error for an individual eigenvector depends on the corresponding eigen-gap yu2015useful. When this gap is small, or when the factor number $r$ grows and $\lambda_1 \asymp \lambda_r$, the performance of standard PCA deteriorates, whereas our RC-PCA remains effective. Moreover, the theoretical investigation of RC-PCA under relaxed assumptions is non-trivial and holds independent value.}. Random projection, also known as random slicing anandkumar2014guaranteed,sun2017provable is a well-recognized initialization method in noiseless tensor CP decomposition, which accommodates closely spaced eigenvalues. We extend this approach to the tensor CP factor model. In Algorithm (ref) and Procedure (ref), the tuning parameters $c_0$ and $\nu$ are not highly sensitive to their specific values. In both the simulation and empirical studies presented in the main text, we use the default settings $c_0 = 0.1$ and $\nu = 0.8$, and we find that small or moderate deviations from these values do not affect the results\footnote{The parameter $c_0$ is related to the noise magnitude. Since the randomized projection step guarantees consistency under orthogonal loadings regardless of whether the eigen-gap condition holds, using a larger $c_0$ is a safe choice. The parameter $\nu$ is a pre-specified threshold used to filter out tuples that are too similar to the tuple already selected. Each tuple consists of vectors that estimate factor loadings for an unknown rank $i \in \{ 1, \ldots, r\}$. Therefore, the role of $\nu$ is to ensure that redundant tuples corresponding to the same factor are removed. Further simulation studies regarding $c_0$ and $\nu$ can be found in Appendix (ref).}. Note that the condition $\min\{|\widehat\lambda_i-\widehat\lambda_{i-1}|,|\widehat\lambda_i-\widehat\lambda_{i+1}| \} > c_0 \widehat\lambda_r $ is checked within Algorithm 1 and automatically determines the appropriate regime based on information from the data.

\SetAlgorithmName{Procedure}{procedure}{List of Procedures}

Following initialization, we refine the estimation using an ISO procedure (Algorithm (ref)). This step aims to enhance estimation accuracy and extract latent factors. The procedure is motivated by the vector factor structure of the denoised ${\cal Y}_t$:

where

$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 we have used the fact that $b_{ik}$ is orthogonal to all $a_{jk}, j\neq i$ by construction. The orthogonalization projection, which takes place in all except the $k$th mode simultaneously in each computational iteration, transforms the tensor ${\cal Y}_t$ to a $d_k\times 1$ vector, reducing dimensions and noise substantially. This transformation enables easy and accurate estimation of the classical vector factor model in equation ($\ref{eq:cp-ideal}$).

Consider a simple illustrative example with $K=2$, $r=2$:

The ideal projection with true $b_{12}$ yields

where ${\cal E}_{t,11}^*={\cal E}_t\times_2b_{12}\in \mathbb{R}^{d_1}$. This projection reduces the dimension from $\mathbb{R}^{d_1\times{d_2} }$ to $\mathbb{R}^{d_1}$, while simultaneously increasing the signal-to-noise ratio.

In practice, we do not observe $b_{ik}$ and iterations can be applied to update the estimations. Given the previous estimates $\widehat a_{ik}^{(m-1)}$, where $m$ is the iteration number, ${\cal Y}_t$ can be denoised via

for $t=1,...,T$, and consequently, updated loading vectors $\widehat a_{ik}^{(m)}$ are obtained through eigenanalysis based on 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}$. The iteration continues until convergence or the maximum number of iterations is reached. The default accuracy is set to $\epsilon = 10^{-5}$; a smaller value indicates more precise estimation, which typically results in a larger number of iterations. When our proposed RC-PCA is used as the first-stage (initial) estimator, Algorithm (ref) converges very rapidly. In our simulation and empirical applications, we set the maximum number of iterations to $M=100$, but convergence is typically achieved in fewer than 5 iterations.

\SetAlgorithmName{Algorithm}{algorithm}{List of Algorithms}

The above estimation procedure assumes that the rank $r$ is known. However, we need to estimate $r$ in practice. We consider two estimation procedures based on the eigenvalue ratio method proposed by ahn2013.

For the first procedure, we unfold the sample contemporary covariance $\widehat\Sigma$ in (ref) to a $d\times d$ matrix $\widetilde\Sigma$. Let $\widehat{\lambda}_{1} \geq \widehat{\lambda}_{2} \geq \cdots \geq \widehat{\lambda}_{r} \geq 0$ be the ordered eigenvalues of $\widetilde\Sigma$. The CP tensor factor model (ref) can also be adapted to a vector factor model (ref) with the same number of factors $r$. Thus, the eigenvalue ratio-based estimator derived from the unfolded covariance matrix $\widetilde\Sigma$ can be defined as

and $r_{\max}$ is a selected upper bound.

Alternatively, we can define the mode-$k$ covariance with the inner product:

Let $\widehat{\lambda}_{1k} \geq \widehat{\lambda}_{2k} \geq \cdots \geq \widehat{\lambda}_{rk} \geq 0$ be the ordered eigenvalues of $\widehat\Sigma_k$. The eigenvalue ratio-based estimator using the inner product can be defined as

where $\widehat{r}_k = \operatorname*{arg\,max}_{1 \leq i \leq r_{\max}} \frac{\widehat{\lambda}_{ik}}{\widehat{\lambda}_{i+1,k}}$. We have adopted the setup in the CP tensor factor model (ref) where the number of spiked eigenvalues $r_k$ remains constant across different mode-$k$ covariance. Further details of these two procedures can be found in Algorithms (ref) and (ref).

It is noteworthy that han2024cp explore a similar tensor CP factor model as ((ref)), albeit within a distinct setting where latent factors are assumed uncorrelated and noise follows a white noise process. Methodologically, their approach rely on the autocovariance between ${\cal Y}_{t-h}$ and ${\cal Y}_{t}$, where $h\geq 1$, whereas our method employs contemporaneous covariance. The autocovariance-based method may not be ideal for datasets with low temporal dependence, such as asset return data, which often exhibit minimal serial correlation possibly due to market efficiency.

Another closely related approach is tensor PCA (TPCA) proposed in babii2022tensor. They consider a CP tensor factor model with orthogonal loading vectors. Unlike Tucker factor models, the identification of CP factor models does not necessarily require orthogonality. Applying TPCA to models with non-orthogonal loadings introduces a bias component of higher order than our first-stage RC-PCA. Even when the loadings are orthogonal, our contemporary variance-based iterative estimation exhibits a faster convergence rate than TPCA due to dimension and noise reduction. A comparison of our estimator with the autocovariance-based estimator and TPCA through simulation will be presented in Section (ref).

Theory

In this section, we study the statistical properties of the algorithms introduced previously. Our theoretical framework offers guarantees for consistency and outlines the statistical error rates for estimating the factor loading vectors $ a_{ik}$, where $1\le i\le r, 1\le k\le K$, given certain regularity conditions. Considering that the loading vector $ a_{ik}$ can only be identified with a change in sign, we employ

to quantify the discrepancy between $\widehat a_{ik}$ and $ a_{ik}$. This measure provides a meaningful and practical way to assess the discrepancy between the estimated and true loadings, and we shall apply it in our simulation as well.

To present theoretical properties of the proposed procedures, we impose the following assumptions.

Assumption (ref) allows for cross-sectional dependence in errors and is closely aligned with the noise conditions presented in seminal works such as bai2002, bai2003, lam2011, lam2012, and others within the factor model literature. It requires that $\xi_{it}$ has exponential-type tails, allowing us to apply large deviation theory. This is a standard condition in the literature on ultrahigh-dimensional data analysis (see, e.g., chang2023modelling,fan2013). When $\vartheta=2$, the condition corresponds to a sub-Gaussian tail. Although it might be possible to replace this assumption with finite moment conditions through a significantly more involved theoretical analysis, we have chosen to retain it to focus on the essentials. For simplicity, we assume that the noise tensor remains independent across time $t$, allowing for weak cross-sectional dependence. While incorporating weak temporal correlation among the noise, as suggested by bai2002, is plausible, it substantially complicates our theoretical analysis. Therefore, we defer this exploration to future research. Nonetheless, our simulation studies show that the proposed methods remain robust even when the noise exhibits weak temporal dependence.

Assumption (ref) ensures the unique identification of all factor loading vectors $ a_{ik}$ up to sign changes. Unlike the eigen decomposition of a matrix, if some $\lambda_i$ are equal, the estimation of the loading vectors $ a_{ik}$ isn't subject to rotational ambiguity but only to the signed permutation of loading vectors. Consistent with the classical literature on factor models (e.g., Assumption A in bai2002, Assumption M(d) in stockwatson1998), Assumption (ref) allows factors to be arbitrarily correlated. Furthermore, Assumption (ref) specifies that the tail probability of $f_{it}$ must exhibit exponential decay. Specifically, when $\gamma_1 = 2$, it implies that $f_{it}$ follows a sub-Gaussian distribution. We exclude the case of identical factors. In extreme cases, such as when $f_{1t}=f_{2t}$, our composite PCA method (steps 3 and 4 in Algorithm (ref)) will fail. However, based on Remark (ref), the random projection method (Procedure (ref)) should still yield reasonable initializations. We do not pursue this direction further in our theoretical analysis, as it is beyond the scope of our current paper.

Assumption (ref) is a widely recognized standard condition that accommodates a broad range of time series models, including causal ARMA processes with continuously distributed innovations, as further detailed in works such as tong1990non, bradley2005, tsay2005analysis, fan2008nonlinear, rosenblatt2012markov, tsay2018nonlinear, among others.

While Assumptions (ref) and (ref) currently assume exponential tails for both noise and factor processes, these conditions can be extended to accommodate polynomial-type tails (under bounded moment conditions) when the number of factors $r$ is fixed, albeit at the cost of a more complex theoretical analysis.

Recall $ A_k$ defined in equation ((ref)) with $ a_{ik}$ as its columns. As $ \| a_{ik}\|_2^2=1$, the correlation among columns of $ A_k$ can be measured by

Similarly we use

to measure the correlation of the matrix $ A = ( a_1,\ldots, a_r)\in \mathbb{R}^{d\times r}$ with $ a_i = \mbox{vec}(\otimes_{k=1}^K a_{ik})$ and $d=\prod_{k=1}^K d_k$. Let $\delta_{\max} =\max\{\delta_1,\cdots \delta_K \}$. Using properties of the Kronecker product, we can show that $\delta\leq \prod_{k=1}^K \delta_k < \min_{1\leq k\leq K} \delta_k\le \delta_{\max}$.

Theorem (ref) below presents the performance bounds, which depends on the coherence (the degree of non-orthogonality) of the factor loading vectors.

The first term of the upper limit in (ref) and (ref) arises from the non-orthogonality of the loading vectors $ a_{ik}$, which can be interpreted as bias. Meanwhile, the second term in (ref) stems from a concentration inequality for random noise and thus reflects a form of stochastic error.

When the eigengap condition is not met, we employ randomized projection to determine the statistical convergence rate as shown in (ref), which is slower than the rate in (ref). A broader result than (ref), permitting a more general eigen ratio $\lambda_1/\lambda_r$ for part (ii), is detailed in the appendix. In practice, since the sample covariance tensor includes both the average of signal-by-noise cross-products and the average of noise-by-noise cross-products, it is uncommon to encounter nearly identical sample spiked eigenvalues. Our simulation study demonstrates that while the original composite PCA provides viable initializations when $\lambda_1 = \lambda_r$, its performance is not as good as that of RC-PCA using Procedure (ref).

Let the statistical error bound of the initialization used in Algorithm (ref) be $\psi_0$ (for example, the right hand side of (ref)), and also let

It is important to note that the error bound $\psi_0$ for initialization is intended for each individual loading vector $ a_{ik}$. When applying Algorithm (ref), which requires the inverse of $\widehat A_k^\top \widehat A_k$, the condition $\sqrt{r}\psi_0\lesssim 1$ ensures a reliable initial estimate of the loading matrix $\widehat A_k$. Although the other components in (ref) may seem complex, they are designed to ensure the error contraction effect in each iteration. This ensures that as iterations progress, the error bound will approach the desired statistical upper bound. A more detailed discussion of initial estimates can be found in Appendix (ref).

The following Theorem (ref) specifies the convergence rate for the estimated factors $f_{it}$\footnote{Asymptotic normality results for estimated factors are developed in our companion paper on diffusion index forecasting: \url{https://papers.ssrn.com/sol3/papers.cfm?abstract_id=5213594}, which builds directly on the theoretical framework established here.}.

Note that $w_i$ is a scalar representing the factor strength. In a strong factor model, consistent with bai2003,stock2002,chen2023statistical, it is typically set as $w_i=\sqrt{d}=\sqrt{d_1d_2\cdots d_K}$, and we use $\widehat w_i \widehat f_{it}/w_i$ as the estimated factors for further interpretation and analysis.

We now demonstrate the feasibility of obtaining a more precise bound by closely examining the leading order term. This process allows us to ascertain the asymptotic behavior of the estimator $ a_{ik}$. Specifically, we will establish that

where $P_{ a_{ik},\perp}=I_{d_k}- a_{ik} a_{ik}^\top$ and $\Theta=(\Theta_{ij})_{r\times r}$, with $\Theta$ defined in Assumption (ref). This enables the determination of asymptotic distributions for linear forms of $a_{ik}$.

The following theorem shows the asymptotic distribution of a linear form of the factor loading vector $u^\top a_{ik}$ for some fixed vector $u$. Note that in the strong factor model, we have $\Theta_{ii}\asymp w_i^2 \asymp d$ for all $1\le i\le r$.

Since $\sigma_{u,ik}^2\asymp d_k$, Theorem (ref)(i) implies that $u^\top\left(\widehat a_{ik}^{{\rm iso}} -\text{sign}(\widehat a_{ik}^{{\rm iso}\top} a_{ik})\cdot a_{ik} \right)=O_{\mathbb{P}}(\sqrt{d_k/(\lambda_r T})$. In contrast, Theorem (ref)(ii) yields $u^\top\left(\widehat a_{ik}^{{\rm iso}} -\text{sign}(\widehat a_{ik}^{{\rm iso}\top} a_{ik})\cdot a_{ik} \right)=O_{\mathbb{P}}(1/\lambda_r) $. In Theorem (ref), we focus on vectors $u$ with the property that $\|P_{ a_{ik},\perp} u\|_2>0$ when $d_k$ is large, which effectively assumes $\sin\angle(u, a_{ik})> 0$. Conversely, when $\sin\angle(u, a_{ik}) =0$, the convergence rate of the estimated linear form is faster, and its asymptotic distribution is a mixture of $\chi_1^2$ distributions.\footnote{Theorem 4.5 complements Theorem 4.4 by addressing the case where $\sin\angle(u, a_{ik}) =0$, i.e., when the direction of interest $u$ is aligned with the loading vector $a_{ik}$. In most economic applications, however, we are interested in inference for a particular component of the loading, in which case $u$ is a unit vector aligned with a coordinate axis. Hence, the condition $\sin \angle(u, a_{ik}) > 0$ is typically satisfied, and Theorem 4.4 remains applicable. As a safeguard, one can check whether $\|P_{\widehat a_{ik}, \perp} u\|_2$ is close to zero using the estimated $\widehat a_{ik}$ and the user-specified $u$. } Similarly, we will establish that

Drawing parallels with traditional PCA is insightful; in PCA, a debiasing process is often necessary to achieve asymptotic normality in linear combinations of the principal components, as discussed in koltchinskii2016asymptotics, koltchinskii2017concentration, koltchinskii2020efficient. For the CP tensor factor model, however, merely meeting the signal strength requirement $T/(d_k \Theta_{ii})\to 0$ is enough to render the bias inconsequential. This observation aligns with findings by bai2003 regarding vector factor models.

The estimators are constructed with a specified rank $r$, although in the theoretical analysis, $r$ is allowed to increase. Practically, $\widehat r$ can be estimated using the generalized eigenvalue ratio-based estimators detailed in Algorithms (ref) or (ref). The asymptotic validity of $\widehat r^{\rm uer}$ and $\widehat r^{\rm ip}$ are established in Theorem (ref) below.

Theorem (ref) derives the consistency of rank estimators $\widehat r^{\rm uer}$ and $\widehat r^{\rm ip}$ based on the eigenvalue ratios. This can be viewed as a generalization of Theorem 1 of ahn2013 from vector factor models to CP tensor factor models.

Simulation

In this section, we conduct empirical comparisons among different methods for estimating loading vectors across various simulation scenarios, and verify the limiting distribution of the estimated loading vectors. We assess the performance of the contemporary covariance-based iterative simultaneous orthogonalization procedure (CC-ISO) proposed in this paper, auto-covariance-based iterative simultaneous orthogonalization procedure by han2024cp (AC-ISO), and TPCA by babii2022tensor. The auto-covariance considered by han2024cp is defined by the following lagged-cross product operator: $$ \boldsymbol{\Sigma}_h = \mathbb{E}\left[ \frac{\mathcal{Y}_{t-h} \otimes \mathcal{Y}_t}{T-h} \right] \in \mathbb{R}^{d_1 \times \cdots d_K \times d_1 \times \cdots \times d_K}. $$ In this section, we fix $h = 1$. The estimation error measures the angle between the estimated loading vector and the true loading vector, computed as: $$ \max_{i,k} \| \widehat a_{ik} \widehat a_{ik}^\top - a_{ik} a_{ik}^\top \|_2 = \max_{i,k} \sqrt{1-(\widehat a_{ik}^\top a_{ik})^2} . $$

Throughout our analysis, the observations ${\cal Y}_t$'s are simulated according to model (ref) with $K=2$ and $r=3$. As mentioned in Section (ref), we fix $c_0 = 0.1$, $\nu = 0.8$, $L = 2r^2$, $\epsilon = 10^{-5}$ and $M = 100$. The true loading vectors are generated as follows. The elements of matrices $\widetilde{A}_k=\left(\widetilde{a}_{1 k}, \ldots, \widetilde{a}_{r k}\right) \in$ $\mathbb{R}^{d_k \times r}, 1 \leqslant k \leqslant K$, are drawn from i.i.d. $N(0,1)$ distributions and then orthonormalized via QR decomposition. To conveniently control nonorthogonality, we introduce the parameter $\eta=\max_{1\le i<j\le r}|a_i^\top a_j|$, analogous to $\delta$ defined in (ref). It can be shown that $\delta \leq (r-1)\eta$. If $\eta=0$, we simply set $A_k=\widetilde{A}_k$; otherwise, we set $a_{1 k}=\widetilde{a}_{1 k}$ and $a_{i k}=\left(\widetilde{a}_{1 k}+\theta \widetilde{a}_{i k}\right) /\left\|\widetilde{a}_{1 k}+\theta \widetilde{a}_{i k}\right\|_2$ for all $i \geqslant 2$ and $1 \leqslant k \leqslant K$, with $\theta=\left(\eta^{-2 / K}-1\right)^{1 / 2}$. By varying $\eta$, we control the correlations between loading vectors. It is evident that an increase in $\delta$ (or equivalently $\eta$) leads to a higher degree of linear dependence among the loading vectors.

The factor processes $f_{it}$ exhibit weak temporal dependence and are generated as an independent AR(1) process, with the factor strength $w_i$ being a scalar depending on $d_1, d_2$ and $r$:

where $\phi=0.1$, $\operatorname{Var}(\epsilon) = 1 - \phi^2 = 0.99$. In this case, the tensor factor model in (ref) represents a typical strong factor model, where $\lambda_i = w_i^2 = (r-i+1)^2 d_1 d_2 / 25 = O(d)$ when $r$ is fixed. The results with $\phi=0.5$ are reported in Appendix (ref).

The following three configurations are adapted from babii2022tensor and han2024cp with modifications made for comparative analysis of the empirical performances among TPCA, AC-ISO and CC-ISO. In these configurations, $(d_1, d_2) \in \{(40,40),(40,60),(60,60)\}$, $T \in \{100,300,500\}$ and $r = 3$.

The following configuration aims to assess the robustness of our proposed algorithm under weak factor structures.

For each configuration, we conduct the experiment 500 times and present the box plots of the results. Figure (ref) shows the estimation errors for CC-ISO, AC-ISO and TPCA under configuration I. Notably, CC-ISO consistently outperforms the other two algorithms across various dimensions. The estimation by AC-ISO deviates significantly from the true value due to the weak signal in the auto-covariance matrix resulting from the weak temporal dependence in the factor process.

In Configuration II, we assess the impact of non-orthogonal factor loadings on estimation. Figure (ref) shows the ratio of the estimation errors of CC-ISO to AC-ISO (first panel) and to TPCA (second panel) across different values of $\eta$ and dimensions. In the first panel, the error ratio of CC-ISO to AC-ISO remains around 0.15, indicating the superior accuracy of CC-ISO. The ratio remains relatively stable because the signal part in AC-ISO, albeit small, also increases with dimensions, limiting its relative improvement. In contrast, in the second panel, the error ratio of CC-ISO to TPCA decreases and converges as dimensions grow. This is because TPCA is unable to recover non-orthogonal factor loading vectors, resulting in stable estimation errors regardless of dimension. Meanwhile, CC-ISO effectively identifies non-orthogonal loadings, leading to estimation errors that diminish and eventually converge to zero.

Figure (ref) shows the ratio of estimation errors of CC-ISO to AC-ISO under configuration III, designed to evaluate the robustness of the proposed CC-ISO algorithm against serial correlation in the error term. We observe that CC-ISO's performance improves monotonically as $T$ increases. In contrast, AC-ISO's performance deteriorates as the serial correlations in the error term strengthen. This decline is due to the contamination of signals in the auto-covariance by the serial correlations in the error terms. However, CC-ISO demonstrates robustness against such serial correlations.

In Figure (ref), we show the box plots of the logarithm of the estimation errors of CC-ISO algorithm under a weak factor configuration. It is evident that the estimation errors decrease as $T$ increases. Additionally, the rate of decrease in estimation errors depends on the value of $\alpha$: a higher $\alpha$ leads to a faster decrease. These results validate the robustness of the CC-ISO algorithm against a certain degree of weak factor structure, in line with the conclusions drawn in Theorem (ref).

We also examine the performance of Randomized Projection (RP) from Procedure (ref) and compare it with composite PCA.

Given the close empirical performances of CC-ISO under both initialization methods under configuration V, our focus shifts to the estimation errors of the initial estimations, as illustrated in Figure (ref). RC-PCA outperforms composite PCA in terms of the accuracy of initial estimations, particularly pronounced when $d$ is smaller and $T$ is larger. This occurs because the sample covariance of $\operatorname{Vec}({\cal E}_t)$ approaches the identity matrix when $d$ is small and as $T$ increases. Consequently, $\widetilde\Sigma$ is more likely to have eigenvalues that are closer together.

The subsequent simulation verifies the robustness of the CC-ISO algorithm against weak mis-specification of the model. The data is generated following the Tucker factor model with $K = 2$:

$$ {\cal Y}_t = \bfm A_1 {\cal F}_t \bfm A_2^\top + {\cal E}_t, $$ where ${\cal F}_t \in \mathbb{R}^{r \times r}$ is the factor process in the Tucker factor model. In the CP factor model, ${\cal F}_t$ is diagonal with the $i^{th}$ diagonal element equal to $f_{it}$. In the mis-specification setting, we allow the off-diagonal entries to deviate from 0. Denote the $(i,j)^{th}$ entry of ${\cal F}_t$ by ${\cal F}_{t,ij}$. Let ${\cal F}_{t,ij} = w_{ijt} f_{ijt}$, where $f_{ijt}$ is generated as specified in (ref). The configuration is as follows:

Figure (ref) shows the results under configuration VI. Given $\alpha$, the estimation error decreases in $T$ or in $d$, which illustrates the robustness of CC-ISO against weak mis-specification.

Next simulation is conducted to verify the results in Theorem (ref)(i). The configuration is as follows:

Figure (ref) shows the QQ plots and histograms of $ \sqrt{T \Theta_{ii} }u^\top\left(\widehat a_{ik}^{{\rm iso}} -\text{sign}(\widehat a_{ik}^{{\rm iso}\top} a_{ik})\cdot a_{ik} \right) / \sigma_{u,ik}$ derived from Theorem (ref) under Configuration VII. It is observed that the normalized empirical distribution closely approximates the standard normal distribution. Figure (ref) in Appendix C shows the corresponding histograms when the variance $\sigma_{u,ik}^2$ is replaced by its estimate $\widehat\sigma_{u,ik}^2$, using the thresholding estimator proposed in chy2025. The approximation to the standard normal distribution remains accurate even with estimated variance, demonstrating the robustness of the inference.

Next, we evaluate the performance of two proposed rank estimation algorithms: the unfolded eigenvalue ratio method and the eigenvalue ratio method via inner product, over the following DGP configuration:

The results are presented in Table (ref). The numbers in the table denote the relative frequency of correct rank estimation over 500 replications. Both methods perform very well with accuracy levels close to 1.

Finally we evaluate the factor estimation of the proposed algorithm and compare it with AC-ISO.

Figure (ref) shows the results under Configuration IX. The CC-ISO method outperforms AC-ISO in terms of both the accuracy and variability of the factor estimates. In addition, the estimation errors decrease as $d$ and $T$ increase, which is consistent with Theorem (ref).

Empirical Application: Characteristic decile portfolios

Latent factor models have been extensively employed in the asset pricing literature to identify risk factors and measure risk exposures. In this section, we conduct an empirical analysis using the dataset compiled by ChenZimmermann2021, which consists of over 200 characteristic-sorted portfolios drawn from previous studies on stock market anomalies. This data has also been analyzed by babii2022tensor. We utilize the August 2023 release of the database, which contains monthly portfolio excess returns sorted into 10 deciles based on firm-level characteristics, spanning from January 1990 to December 2022. To ensure a balanced panel of portfolios throughout the sample period, we restrict our analysis to 127 characteristics\footnote{The issue of missing values in tensor factor models is both important and challenging. To the best of our knowledge, the only existing study that addresses this issue is cen2025tensor, which adopts a Tucker tensor factor structure; no prior work has considered missing values under the CP tensor factor framework. Tackling this problem entails substantial additional effort, which we plan to undertake in a separate project.}. Consequently, the dimension of the tensor-valued time series ${\cal Y}_t$ is $127 \times 10$ with a sample size of $T=396$. The risk-free rate, obtained from the Kenneth French data library, is used to compute the excess return of each portfolio.

As noted by babii2022tensor, pelger2020RFS also study a similar dataset, which includes 10 single-sorted decile portfolios for each of 37 characteristics. However, rather than utilizing the natural tensor structure, they flatten the data into a total of 370 portfolios and apply risk-premium PCA based on a vector factor model. In contrast, we exploit the inherent tensor structure by estimating a CP tenor factor model via the proposed CC-ISO method. The CP tensor factor model effectively identifies risk factors and factor exposures without rotation ambiguity. This key property ensures that the estimated factors and loadings are uniquely identified up to sign, independent of rotations or the choice of a varimax algorithm, therefore enhancing the reliability of their economic interpretation.

To implement this, we estimate the CP tensor factor model introduced in Equation (ref). As discussed in Section (ref), the latent factors can be interpreted as systematic risk factors, while the factor loadings represent heterogeneous risk exposures.

To determine the number of factors, we examine the eigenvalues of the sample covariance matrix. The first factor is substantially stronger than the others, with its corresponding eigenvalue accounting for 86% of the total sum of eigenvalues, which might suggest the presence of a dominant market factor. To ensure that our model captures effects beyond this dominant factor, we set the minimum number of factors to two. Guided by the eigenvalue ratio plot and the scree plot (Figure (ref)), we select three factors for our analysis. For the Tucker factor models, the mode-wise eigenvalue ratio and scree plots (Figure (ref) in Appendix D) suggest the rank pair $(2,3)$.

We estimate the CP factor model in Equation $\eqref{eqn:fm}$ using various algorithms: CC-ISO, TPCA, AC-ISO with $h=1$, GE proposed by chang2023modelling with $K = 1$, the CC-based iterative projection methods based on Tucker decomposition (Tucker-CC-IP) proposed by Han2021iterative with rank $(2,3)$, the HOOI algorithm based on Tucker decomposition applied by lettau20243d (3D-PCA), and the observed factor model with Fama-French 3 factors (FF3). We compare the in-sample performance by $R^2$ defined as: $$ R^2 = 1 - \frac{ \| \mathcal{Y} - \widehat{\mathcal{Y}} \|_{F}^2 }{ \| \Tilde{\mathcal{Y}} \|_{F}^2 }, $$ where $\widehat{\mathcal{Y}}$ is the fitted value of $\mathcal{Y}$ with various methods and $\Tilde{\mathcal{Y}}$ is the demeaned tensor of $\mathcal{Y}$ with the $t^{th}$ entry defined as $\Tilde{\mathcal{Y}}_{t} = \mathcal{Y}_{t} - \frac{1}{T} \sum_{t=1}^T \mathcal{Y}_{t}.$

To evaluate out-of-sample performance, we conduct one-step-ahead forecasts over the final two years of the sample, following the scheme proposed by chang2023modelling. Specifically, for each $s \in \{1,2,\ldots,24\}$, we estimate the model using the sample $\{{\cal Y}_t\}_{t=1}^{371+s}$ and fit the estimated latent factors with a VAR(1) model. The one-period-ahead excess return ${\cal Y}_{372+s}$ is then predicted using the estimated loadings and the predicted one-period-ahead factors. The out-of-sample MSE is defined as: $$ \text{Out-of-sample MSE} = \frac{1}{24} \sum_{s=1}^{24}\left\| {\cal Y}_{372+s} - \widehat {\cal Y}_{372+s} \right\|_F^2. $$

Table (ref) displays the $R^2$ and out-of-sample MSE ratios (normalized by Tucker-CC-IP) for six algorithms along with $p$-values from the one-sided Diebold1995 (DM) Forecast comparison tests. In the DM tests, row “DM I” evaluates whether each competing algorithm outperforms Tucker-AC-IP, while row “DM II” tests whether CC-ISO outperforms the respective competing method. Among the CP-based methods, CC-ISO achieves the best in-sample fit, with the highest $R^2$. Within the ISO methods, CC-ISO outperforms AC-ISO. In contrast, The GE algorithm yields the lowest $R^2$ value of 0.768.

Given that Tucker decomposition generally involves a much larger number of factors than CP decomposition ($6$ vs $3$ in our example), it is not surprising that Tucker-CC-IP achieves slightly higher $R^2$ than the CP-based methods. Nonetheless, the improvement is modest, with Tucker-CC-IP increasing $R^2$ by only $2$ to $3$ percentage points compared to CP models. This observation suggests that the characteristic decile portfolio data might exhibit a CP-like factor structure, where a low-rank CP model effectively captures the underlying risk factors.

In terms of forecast performance, smaller MSE ratios indicate better prediction accuracy. CC-ISO significantly outperforms Tucker-CC-IP, 3D-PCA, and GE, while performing comparably to AC-ISO and TPCA. Interestingly, although TPCA achieves the lowest MSE, its improvement over Tucker-CC-IP is not statistically significant. Overall, the Tucker-based methods tend to produce higher MSEs, likely due to overfitting from the increased number of factors. These findings highlight the advantage of CP-based methods in balancing model parsimony with predictive performance.

Nest, we evaluate cross-sectional pricing errors in the spirit of lettau20243d. Specifically, for each decile portfolio, we estimate

$$ r_{ijt} = \alpha_{ij} + \beta_{ij}' f_t + u_{ijt}, \quad i = 1, \ldots, 127, j = 1, \ldots, 10, $$ and compute the mean squared pricing errors (MSPE) of $f_t$ as $$ MSPE_{f} = \frac{1}{127 \times 10} \sum_{i,j} \alpha_{ij}^2. $$

For the Tucker factor model with factor matrix $F_t \in \mathbb{R}^{r_1 \times r_2}$ from Tucker factor model, we vectorize it prior to regression, i.e. $f_t = \Vec(F_t) \in \mathbb{R}^{r_1r_2}$.

Table (ref) reports MSPE and MSPE ratios of different methods relative to Tucker-CC-IP with rank (2,3), as suggested by the mode-wise scree plots. For completeness, we include the 3D-PCA approach of lettau20243d and the Fama-French three-factor model. As shown, the factors obtained from CC-ISO and GE achieve the lowest MSPE, demonstrating that our model performs strongly not only in out-of-sample forecasting but also in cross-sectional pricing.

Finally, we further examine the estimated factors and loadings in our CP factor model. The estimated loadings reveal several notable patterns. Figure (ref) presents the estimated loadings along with their 95% confidence intervals across characteristics. The loadings on the first factor are uniformly positive, resembling “long-only” portfolio weights. The loadings on the second and third factors display more variation, with predominantly positive elements. These loadings are closely linked to the volatilities of the decile portfolios. In particular, Figure (ref) compares the volatilities of decile 10 portfolios with the corresponding characteristic loadings on the third factor. The figure highlights a strong correlation (0.875) between the two: characteristics with higher decile 10 portfolio volatility tend to have larger third-factor loadings, while those with lower volatility have smaller or even negative loadings.

The factor loadings along the decile mode follow a distinct level-slope-curvature pattern, as shown in Figure (ref). The first factor follows a “long-only” structure, with strictly positive and nearly uniform weights across deciles. The second factor shows a slope pattern, with decile-10 and decile-1 loadings having opposite signs but similar magnitudes. This structure resembles a “long-short factor”, where high-return deciles are offset by low-return deciles. The third loading shows a convex pattern, forming the “curvature” component.

The consistent “long-only” pattern of the first loading across both modes, combined with the dominant eigenvalue of the first factor, suggests a strong association with the market factor. Figure (ref) presents the time series of the first factor along with the excess market return from the Kenneth French data library. The two series exhibit a high correlation (above 0.9), reinforcing the interpretation that this factor primarily captures systematic market risk. In summary, the three estimated factors from the CP factor model can be interpreted as the market factor, the long-short factor, and the volatility factor, respectively.

Conclusion

Modeling high-dimensional tensor time series has attracted growing attention, driven by the increasing availability of multidimensional datasets that go beyond the classical panel data structure. This paper considers matrix and tensor factor models with a CP low-rank structure, extending traditional vector factor models to higher-order settings. We develop ISO procedures based on contemporary covariance, preserving the tensor data structure. Theoretical properties such as the rate of convergence and limiting distributions are investigated under settings where each tensor dimension is comparable to or exceeds the number of observations, and the tensor rank may be either fixed or divergent.

Unlike auto-covariance-based estimation methods, our method explores information from contemporary covariance and can consistently estimate both loadings and factors even when observations are uncorrelated, settings where auto-covariance methods may fail. Additionally, we propose two generalized eigenvalue-ratio estimators for rank selection and justify their consistency. Extensive simulation studies highlight the merits of our method over existing methods. Furthermore, the empirical application to sorted portfolios demonstrates the practical relevance of our approach.

\setcounter{page}{1}