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.
80,660 characters · 18 sections · 81 citation commands
Tensor PCA for Factor Models
\thispagestyle{empty}
\setcounter{page}{0}
Modern empirical analysis often relies on high-dimensional panel datasets with non-negligible dependencies across one or several dimensions. Various economic reasons exist for such dependencies, including common shocks in macroeconomics and finance, interactions in a network, or spatial distribution of economic activity. Since their original introduction in the psychology literature by spearmen1904general, factor models have become a ubiquitous tool in economics for modeling dependencies in two-dimensional panel datasets.
Among the many applications of factor analysis in macroeconomics and finance, we may cite asset pricing, see ross1976arbitrage and chamberlain1982arbitrage; business cycle analysis sargent1977business; industrial production analysis, see foerster2011sectoral and andreou2019inference; forecasting with big data; see stock2002forecasting. Factor models are also widely used for causal inference in panel data; see pesaran2006estimation, bai2009panel, abadie2010synthetic, gobillon2016regional, athey2021matrix, or more generally to model unobserved heterogeneity in microeconometrics; see cunha2010estimating and bonhomme2010generalized. Lastly, structural economic models often naturally lead to a factor structure as seen in liu2024dynamic for the input-output production networks or lewbel1991rank for the demand systems.
Traditional factor models are designed for two-dimensional panel datasets consisting of cross-sectional units tracked over time, where both dimensions can potentially be large; see stock2002 and bai2003inferential.\footnote{We will use the term `traditional' for $2$-dimensional factor models applicable to 2-dimensional panel datasets.} However, many economic data feature more than two dimensions. For example, the input-output production or trade network usually evolves over time, which is not captured by the traditional factor model.\footnote{See, e.g.\ graham2020network for a recent review of static networks.} A panel of macroeconomic indicators often involves regional aggregation of state- or county-level observations, introducing a geographical dimension in addition to the traditional cross-section and time series. Asset pricing models for the cross-section of equities typically involve characteristic-based portfolio sorts, introducing a third dimension. The sorting into deciles is common, but only the lowest and highest deciles are used and combined in a high minus low return spread. Going beyond national borders introduces an international dimension to macro and financial datasets, resulting in a four-dimensional data structure. Each of these examples illustrates that to obtain matrix representations of the data, we often aggregate the multidimensional panel datasets, suppressing more granular information.
Principal component analysis (PCA) is a commonly used method for identifying and estimating traditional factor models for two-dimensional panel data sets; see pearson1901liii. PCA extracts latent factors and their loadings using either the singular value decomposition (SVD) of the original panel dataset collected in a matrix or, equivalently, the eigendecomposition of the associated sample covariance matrices; see jolliffe2002principal for a review of PCA and factor analysis.
A $d$-way tensor is a $d$-dimensional array generalizing vectors and matrices introduced in ricci1900methodes.\footnote{Tensor tools have already been used by economists, e.g., for the identification of finite mixture models; see bonhomme2016estimating.} In this paper, we consider an extension of traditional factor models to multidimensional datasets, called tensor factor models. Similarly to their $2$-way counterpart, the $d$-way tensor factor model can be used to identify the latent factors driving correlations in tensor datasets. In traditional factor models, correlations are captured with a matrix of factors in the time dimension and a matrix of factor loadings in the cross-sectional dimension. Likewise, one can decompose a tensor into a collection of matrices with respect to each dimension. Hence, the matrix component in the `time' dimension correspond to factors that vary over time, and the matrix component in the other dimensions correspond to loadings determining the heterogeneous exposure of each dimension to the latent factors.
Traditional factor models can also be viewed as decomposing a matrix representing a panel dataset as a sum of a low-rank matrix product of loading and factor matrices and a matrix of idiosyncratic shocks. The low-rank component is then approximated using the truncated SVD decomposition of the observed data. The $d$-way factor model also describes a $d$-way tensor dataset as a sum of a low-rank component and an idiosyncratic shocks tensor. However, there several different ways to describe the rank of a tensor. The two most widely known definitions of tensor ranks are called the Canonical Polyadic (CP) rank and the multilinear rank; see carroll1970analysis, harshman1970foundations, and tucker1966some.
The multilinear rank of a $d$-dimensional tensor is described by the $d$-tuple $(R_1,R_2\dots,R_d)$ corresponding to ranks of each of its mode-$j$ matricizations and is more general than than the CP rank.\footnote{The concepts of multilinear rank and CP rank for tensors are implicitly used in hitchcock1927expression. In fact, hitchcock1927expression considered even more general concept of a tensor rank based on arbitrary matricizations, not just along one of its $d$ modes.} It leads to the so-called Tucker model that has been previously considered in statistics literature; see han2020tensor, wang2022high, chen2022factor, han2022rank, and zhang2018tensor among others. Our paper introduces the factor dynamics differently in contrast to these papers. More importantly, we show that under the strong factor model assumption, commonly used in economics and finance, the simple PCA-type estimators for tensors have optimal convergence rates and derive the corresponding asymptotic distributions.
The simple tensor PCA (TPCA) estimators are easy to compute and consists of reshaping (or matricizing) a $d$-way tensor into $d$ matrices along each of its $d$ dimensions and applying the standard PCA these matrices. For this simple TPCA estimators, we (a) demonstrate that they can identify and consistently estimate the factors and loadings; (b) describe the associated convergence rates and asymptotic distribution; (c) show that tensor dimensions can improve the estimation accuracy for factors/loadings and the estimator is rate-optimal for strong factor models; and (d) develop a formal test for the number of factors in the tensor factor model. For the moderately weak factor model, the TPCA procedure is not optimal and we also consider additional results for the optimal iterative alternating least-squares (ALS) algorithm, where the TPCA estimators are used as a starting value. Our theoretical results rely on the powerful perturbation theory results for singular subspaces and singular values which are not commonly used in econometrics, cf. bai2023approximate and references therein.
Monte Carlo simulations support our asymptotic results in finite samples. We find that the $d$-way tensor factor model reduces dimensions more efficiently than the naively pooled $2$-way factor model. We also verify the convergence rates and the distribution theory.
In the empirical application we study the issue of missing firm characteristics in widely used data sources which are the cornerstone of research in finance. bryzgalova2022missing provide a comprehensive analysis of missing data in firm characteristics, and propose a statistical model for imputing missing values and investigate the impact of missingness on asset returns. Firm characteristic data is a 3-dimension tensor, with (a) firm, (b) characteristic and (c) time as the three dimensions. The statistical methods proposed by bryzgalova2022missing do not exploit the tensor data structure, while in this empirical application we do and we show that using our proposed tensor data models improve on their time series of cross-sectional models.
The paper is organized as follows. Section (ref) introduces a new class of tensor factor models. The convergence rates and large sample distributions of the TPCA estimator are covered in Section (ref) which also covers the improved ALS estimator. Section (ref) presents a novel testing procedure for the number of factors. Small sample simulation evidence is reported in Section (ref). An illustrative empirical example appears in Section (ref). Section (ref) concludes. Lastly, in the Appendix, we provide proofs for all the main and auxiliary results.
\paragraph{Notation:} For two matrices $A\in\mathbb{R}^{N_1\times R_1}$ and $B\in\mathbb{R}^{N_2\times R_2}$, we use $A\otimes B\in\mathbb{R}^{N_1N_2\times R_1R_2}$ to denote their Kronecker product. In addition, for a collection of matrices $(\Lambda_j)_{j=1}^d$, we define $\bigotimes_{k\ne j}\Lambda_k=\Lambda_d\otimes\dots\otimes \Lambda_{j+1}\otimes \Lambda_{j-1}\otimes\dots\otimes \Lambda_1$. For a positive integer $p$, we use $I_p$ to denote the $p\times p$ identity matrix. For a matrix $A$, we use $P_A = A(A^\top A)^\dagger A^\top$ to denote the the projection matrix on the column space of $A$, where $^\dagger$ denotes the generalized inverse. For two sequences $(a_n)_{n\in\ensuremath{\mathbb{N}}}$ and $(b_n)_{n\in\ensuremath{\mathbb{N}}}$, we write $a_n\lesssim b_n$ if and only if there exists $C<\infty$ such that $a_n\leq Cb_n$ for all $n\in\ensuremath{\mathbb{N}}$. The operator norm of a matrix $A$ is defined as $\|A\|_{\rm op}=\sup_{\|x\|=1}\|Ax\|$, where $\|.\|$ is the Euclidean norm. More generally, we use $\|.\|_p$ to denote the $\ell_p$ norm. For two tensors $A,B\in\ensuremath{\mathbb{R}}^{N_1\times\dots\times N_d}$, the Frobenius inner product is defined as $\langle A,B\rangle_F = \sum_{i_1,\dots,i_d}A_{i_1,\dots,i_d}B_{i_1,\dots,i_d}$. Let $\|A\|_{\rm F} = \sqrt{\langle A,A\rangle_{\rm F}}$ be the Frobenius norm of a tensor $A$ induced by the inner product. For a matrix $A$ with columns $(a_1,\dots,a_n)$, the $\ell_{2,1}$ matrix norm defined as $\|A\|_{2,1}=\sum_{j=1}^n\|a_j\|$. We also use $\mathbb{O}_{N}$ to denote the set of $N\times N$ orthogonal matrices, i.e. matrices $A$ such that $A^\top A = I_{N}$, where $I_N$ is $N\times N$ identity matrix. Lastly, for $a,b\in\ensuremath{\mathbb{R}}$, put $a\wedge b=\min(a,b)$ and $a\vee b = \max(a,b)$.
Traditional factor models apply to $2$-dimensional panel data represented by a matrix $\mathbf{Y}\in \ensuremath{\mathbb{R}}^{N\times T}$. The factor model with $R$ factors can be expressed as a sum of a low-rank matrix and a matrix of idiosyncratic shocks:
where $F\in\ensuremath{\mathbb{R}}^{T\times R}$ is a matrix of $R$ factors, $\Lambda\in\ensuremath{\mathbb{R}}^{N\times R}$ is a matrix of corresponding factor loadings, and $\mathbf{U}\in\mathbb{R}^{N\times T}$ are idiosyncratic random shocks.\footnote{The factors and loadings can be random, but all results are stated conditional on their realization.} The estimation of factors and their loadings can be done via PCA which can be computed with the singular-value decomposition (SVD) of $\mathbf{Y}$. More precisely, the fators/loadings can be estimated as the leading $R$ right/left singular vectors of $\mathbf{Y}$. Equivalently, the loadings can be estimated as the leading $R$ eigenvectors of $\mathbf{Y}\mathbf{Y}^\top$ while and the factors as the leading $R$ eigenvectors of $\mathbf{Y}^\top\mathbf{Y}$.
A $d$-dimensional panel dataset can be represented by a $d$-way tensor, denoted $\mathbf{Y}\in\ensuremath{\mathbb{R}}^{N_1\times\ldots\times N_d}$. We can describe $\mathbf{Y}$ by enumerating all its elements along the $d$ ways (or modes):
see Figure (ref) for a graphical illustration. Similarly to a 2-way factor model, we can define a $d$-way factor model for a $d$-way tensor $\mathbf{Y}\in\mathbb{R}^{N_1\times\dots\times N_d}$ as follows
Unlike for matrices, the rank of a tensor can be defined in several different ways. To describe the notion of rank used in this paper, we need to introduce the notion of tensor matricization.
A tensor matricization operation can be described in terms of tensor fibers---a generalization of matrix rows and columns. A fiber is defined by fixing all but one of its dimensions, e.g., a matrix column is a mode-$1$ fiber and a matrix row is a mode-$2$ fiber. The mode-$j$ fibers of a higher-order tensor are defined similarly; see Figure (ref) for an illustration in the case of a $3$-way tensor.
Recall that a matrix $\mathbf{Y}\in\mathbb{R}^{N_1\times N_2}$ can be vectorized by stacking all its columns to obtain an $N_1N_2\times 1$ vector. There is also a second way to vectorize a matrix by stacking all its rows into an $1\times N_1N_2$ vector. Similarly, a $d$-way tensor $\mathbf{Y}\in\mathbb{R}^{N_1\times\dots\times N_d}$ can be matricized in $d$ different ways across each of its $d$ modes.\footnote{Sometimes the term “flattening” of a tensor is used. We prefer the term “matricization” as it makes clear we are creating matrices, i.e., tensors could be flattened to lower dimensions that are not necessarily matrices.} The mode-$j$ matricizations of a tensor are obtained by stacking its mode-$j$ fibers as columns of a $N_j\times\prod_{l\ne j}N_l$ matrix, denoted $\mathbf{Y}_{(j)}$. Formally, the elements of a mode-$j$ matricization of a tensor $\mathbf{Y}\in\mathbb{R}^{N_1\times\dots\times N_d}$ is obtained by the following mapping:
where $\dot y_{i_j,k}$ is $(i_j,k)$ element of $\mathbf{Y}_{(j)}\in\mathbb{R}^{N_j\times\prod_{l\ne j}N_l}$; see Appendix (ref) for a numerical example of how a $3$-way tensor is matricized along each of its three ways.
As we have already mentioned there are several different definitions of tensor rank, each leading to a different $d$-way factor model. The two most widely used definitions of tensor rank go back at least to hitchcock1927expression and lead to the so-called Canonical Polyadic (CP) and Tucker decompositions.\footnote{The CP decomposition is also known as CANDECOMP or PARAFAC; see carroll1970analysis, harshman1970foundations. The name Tucker decomposition comes from tucker1966some.} The Tucker decomposition is more general than the CP decomposition and is based on a notion of a multilinear rank which corresponds to the $d$-tuple of ranks of each of its mode-$j$ matricizations:
Note that the matricizations of a tensor can have in general different ranks which is allowed in the Tucker decomposition and is not allowed in the CP decomposition. We can decompose a tensor $\mathbf{X}\in\mathbb{R}^{N_1\times\dots\times N_d}$ with multilinear rank $(R_1,\dots,R_d)$ as
or more concisely as
where $\Lambda_j\in\mathbb{R}^{N_j\times R_j}$ is a matrix of factors/loadings with entries $\lambda_{i_j,r_j}^{(j)}$; $\mathbf{G}\in\mathbb{R}^{R_1\times\dots\times R_d}$ is the so-called core tensor with elements $g_{r_1,\dots,r_d}$; and $\times_j$ is the mode-$j$ product. The mode-$j$ product is defined as a multiplication along the mode-$j$, e.g., the mode-$1$ product is $\times_1:\mathbb{R}^{R_1\times R_2\dots\times R_d}\times \mathbb{R}^{N_1\times R_1}\to \mathbb{R}^{N_1\times R_2\dots\times R_d}$ is
Following the idea of defining a tensor factor model as a sum of a low-rank tensor and a tensor of idiosyncratic shocks, we define the Tucker tensor factor model as
It is known that the Tucker decomposition is in general not unique, however, it is possible to identify the core tensor and all factors/loadings in the Tucker tensor model under additional orthogonality restrictions; see Section (ref). The orthogonality restrictions are also commonly used to identify 2-way factor models.
In the remaining part of this section, we consider several examples of tensor data in economics and finance, showing that such type of data appear in many applications. Using 3-dimensional examples with the notation $\mathbf{Y}$ = $y_{i,j,k}$ or $y_{i,j,t},$ they include:
As well as 4-dimensional examples such as:
In this section, we consider two tensor PCA (TPCA) algorithms to estimate the Tucker tensor factor model in equation ((ref)). We begin by discussing the sufficient identifying conditions and present a simple TPCA algorithm. We argue that the simple TPCA is optimal when factors are strong in the second subsection. The third subsection describes an improved iterative alternating least-squares algorithm and discusses the improvements for the weak factor model. The next subsection provides the large sample distributions for loadings/factors estimated with simple TPCA.
It is known that the Tucker decomposition of a tensor is not in general unique. In this section, we argue that the loadings/factors in the Tucker tensor factor model in equation ((ref)) can be identified under the following assumption:
By Appendix Lemma (ref), the mode-$j$ matricization of equation ((ref)) is
where $\bigotimes_{l\ne j}\Lambda_l=\Lambda_d\otimes\dots\otimes \Lambda_{j+1}\otimes \Lambda_{j-1}\otimes\dots\otimes \Lambda_{1}$. The latter satisfies the following property:
Then the tensor factor model is identified from the PCA applied to $\mathbf{Y}_{(j)}$. Indeed, if $\mathbf{U}=0$, then by Proposition (ref)
where under Assumption (ref), (ii) $D_j=\mathbf{G}_{(j)}\mathbf{G}_{(j)}^\top$ is a diagonal matrix with distinct elements and $\Lambda_j^\top\Lambda_j=I_{R_j}$. Then $\Lambda_j$ is identified as eigenvectors of $\mathbf{Y}_{(j)}\mathbf{Y}_{(j)}^\top$ or equivalently as the left singular vectors of $\mathbf{Y}_{(j)}$. Moreover, we have
since $\Lambda_j^\top \Lambda_j=I_{R_j}$ for all $1\leq j\leq d$. If $\mathbf{U}\ne 0$, then assuming that the factor part dominates the noise asymptotically, cf. Assumption (ref) below, identification is achieved in large samples.
This leads us to the TPCA estimation Algorithm (ref) often credited to tucker1966some with some aspects already in hitchcock1927expression.\footnote{It is also known as higher-order singular values decomposition in the numerical analysis; see de2000multilinear.}
The following assumption imposes mild restrictions on the data generating process.
Assumption (ref) does not impose any restrictions on the dependence structure in factors and loadings. The factors may be generated by a non-stationary process and the loadings may fail to be i.i.d. It is plausible that after controlling for cross-sectional and time-series dependence through the factor structure, the idiosyncratic errors are i.i.d. Nonetheless, Assumption (ref) can also be relaxed to heterogeneous and dependent arrays at costs of heavier notation and proofs and modified TPCA algorithm, see zhang2022heteroskedastic, which is left for future work.
Next, we make an assumption about the factor strength. Let $\sigma_{j,R_j}$ be the smallest non-zero singular value of the matricized tensor $\mathbf{X}_{(j)}$ and let $\delta=\min_{1\leq j\leq d}\sigma_{j,R_j}$ be a measure of factor strength in the tensor factor model. We assume next that factors are not extremely weak in the sense that the smallest singular values of matricizations are well-separated from zero:
The strong factors model in the tensor setting can be described as $\delta^2\sim \prod_{j=1}^dN_j$.\footnote{In the 2-way case, this corresponds to assuming that the squared singular values of $\mathbf{Y}$ scale at the $NT$-rate; see also onatski2012asymptotics,onatski2022uniform for the weak factors in the 2-way case.} In this case, Assumption (ref) is automatically satisfied. Let $O_j$ be the optimal rotation matrix solving $\min_{O\in\mathbb{O}_{R_j}}\|\hat \Lambda_jO - \Lambda_j\|_{\rm F}^2$.\footnote{It is known that the optimal rotation matrix for the Frobenius norm has a closed-form expression $O_j =\mathrm{sign}(\hat\Lambda_j^\top\Lambda_j)$, where the sign of a matrix $A$ with SVD $A=\Lambda DV^\top$ is defined as $\mathrm{sign}(A)=\Lambda V^\top$.} The following result holds:
The proof appears in the Appendix. Theorem (ref) is a non-asymptotic result valid for any values of tensor dimensions $(N_1,\dots,N_d)$, provided that Assumption (ref) is satisfied. It also applies to the CP-factor model, when matricizations have the same ranks and factors/loadings are assumed to be orthogonal. The orthogonality restriction for the CP factor model can also be relaxed; see chang2023modelling and chen2024estimation for recent contributions.
For weak factors/loadings, the estimation accuracy of factors/loadings in Theorem (ref) does not always improve when more data is available in all tensor dimensions because the second term with $\prod_{l=1}^dN_l$ can dominate, creating a bottleneck. This comes from the fact that the matricized tensor $\mathbf{Y}_{(j)}$ has a very large number of columns, namely $\prod_{l\ne j}N_l$. Under Assumption (ref), using the mode-$l$ multiplication of $\mathbf{Y}$ by loadings $\Lambda_l,l\ne j$, we obtain
cf. equation ((ref)). The mode-$j$ matricization of equation ((ref)),
has a much smaller number of columns, namely $\prod_{l\ne j}R_l$ instead of $\prod_{l\ne j}N_l$. Therefore, applying PCA to the mode-$j$ matricization of the left-hand side in the equation ((ref)), may allow us to estimate $\Lambda_j$ more precisely, provided that the noise is negligible; see also equation ((ref)). Iterating this procedure and updating the factors/loadings leads to the alternating least-squares Algorithm (ref) for Tucker decomposition; see kroonenberg1980principal.\footnote{It is also known as the higher-order orthogonal iterations; see de2000best for more details.} The following result establishes the convergence rate of estimators computed using this algorithm:
The proof of this result appears in the Appendix.
Theorem (ref) makes a simplifying assumption that the tensor dimensions increase proportionally to each other which in the 2-way factor case requires that $N/T\to c>0$. Then for the strong factor model, we have $\delta^2\sim N^d$ and the condition $N^d/\delta^4=o(1)$ simply requires that $N\to\infty$ which is trivially satisfied. Consequently, we obtain the rate of order $O_P(1/N^{d-1})$ which is the same as the optimal rate of the naive TPCA Algorithm (ref). However, when the factors are weak, the rate of the ALS Algorithm (ref) can be faster than that of the naive TPCA. Therefore, we can expect that the ALS with naive TPCA used as a starting point may estimate factors/loadings more accurately, provided that a sufficiently large number of iterations is made.
In this section, we consider inference on loadings and factors in the general Tucker model, covering the orthogonal CP model as a special case. The loadings/factors are estimated with the TPCA Algorithm (ref). Let $\hat\Lambda_{j,i:}$ and $\Lambda_{j,i:}$ be $R_j\times 1$ vectors corresponding to the transposed $i^{\rm th}$ rows of $\hat\Lambda_jO_j$ and $\Lambda_j$ respectively with $O_j=\hat\Lambda_j^\top\Lambda_j$. Let $\|.\|_{2,\infty}$ be the entry-wise $\ell_{2,\infty}$ matrix norm.
We need the following set of assumptions for inference:
The coherence condition roughly requires that the values of factors/loadings are spread out instead of being concentrated in a finite number of entries; see candes2012exact, Definition 1.8.
The following result holds for factors/loadings, depending on the dimension $j=1,\dots,d$:
Note that for $d=3$, conditions of Theorem (ref) are satisfied when the tensor dimensions grow proportionally, i.e., $N_j\sim N\uparrow\infty$ for $j=1,2,3$. Note also that the results of bai2003inferential, Theorem 2, for 2-way factor models are not applicable to tensors because the condition $\sqrt{T}/N=o(1)$ does not hold for $d=3$ when the tensor dimensions grow proportionally since $T\sim N^2$ in this case. Lastly, eliminating the bias and/or relaxing the i.i.d. homoskedastic errors may require modifying the TPCA estimator, see zhang2022heteroskedastic, and is left for future research.
We now turn to the estimation of scale components. Let $(\hat \sigma_{j,r}^2)_{1\leq r\leq R}$ be the eigenvalues of $\mathbf{Y}_{(j)}\mathbf{Y}_{(j)}^\top$. The following result holds:
It immediately follows from Theorem (ref) that $\hat\sigma_{j,r}^2/\prod_{j=1}^dN_j$ is a consistent estimator of the asymptotic factor strength constant $d_{j,r}$.
In this section, we develop a novel test for the number of factors in the tensor factor model. Specifically, we consider the following hypotheses:
The quantity $K$ is selected by the user and reflects the upper bound on the total possible number of factors.
Let $\hat\sigma^2_{j,1}\geq \hat\sigma^2_{j,2}\geq \dots\geq \hat\sigma^2_{j,N_j}\geq 0$ be the eigenvalues of $\mathbf{Y}_{(j)}\mathbf{Y}_{(j)}^\top$. Under Assumption (ref), the first $R_j$ eigenvalues diverge from zero, at least in population. Consider the following statistics
and put
where $(\xi_1,\dots,\xi_{K-k+2})$ follow the joint type-1 Tracy-Widom distribution; see karoui2003largest and soshnikov2002note.
Then, consider the sequence: {
}
The following result holds provided that $N_j\lesssim \prod_{k\ne j}N_k$:
Note that the rate condition for $\tau$ is satisfied in the 3-dimensional case when $N_1\sim N_2\sim N_3$. Theorem (ref) leads to the following testing procedure:
The p-values would correspond to the null hypothesis that there are at most $k$ factors for the mode-$j$ matricization. The joint tests for the number of factors across all matricizations can be obtained using the Bonferroni's correction. The $\alpha$-level test would reject $H_0$ whenever
In practice, one could run the test for several different values of $k$ to determine the number of factors and control the size with Bonferroni correction. Note also that the test statistic $S_j$ is decreasing with $k$, which implies that p-values are increasing with $k$.
The objective of this section is to assess the finite sample properties of our estimation procedure.
We consider the Tucker model for the 3-way tensor for $N\times J\times T$ tensor
where the core tensor is $\mathbf{G}\in\mathbb{R}^{R_1\times R_2\times R_3}$, the loadings matrices are $\Lambda\in\mathbb{R}^{N\times R_1}$ and $M\in\mathbb{R}^{J\times R_2}$, and the factor matrix is $F\in\mathbb{R}^{T\times R_3}$.
We set $R_1=1$ and $R_2=R_3=2$ and generate the $1\times 2\times 2$ core tensor $\mathbf{G}$ with $1\times 2$ slices
where $\sigma_1>\sigma_2>0$. The corresponding matricizations are
Clearly, we have $\mathbf{G}_{(j)}\mathbf{G}_{(j)}^\top$ for $j=1,2,3$ begin diagonal, hence, the core tensor satisfies Assumption (ref). The smallest singular values of $\mathbf{G}_{(j)}$ are $\sqrt{\sigma_1^2+\sigma_2^2}$, $\sigma_2$, and $\sigma_2$ respectively for $j=1,2,3$. The factor strength is $\delta=\sigma_2$. We first set $\sigma_1=d_1\sqrt{NJT}$ and $\sigma_2=d_2\sqrt{NJT}$ which corresponds to the strong tensor factor model.
Therefore, we simulate the $\mathbf{Y}$ tensor with $(i,j,t)^{\rm th}$ observation
Note that this rank-$(1,2,2)$ Tucker tensor factor model features fewer parameters than the rank-2 CP model, where the first loading vector would be different in the two terms.
The rest of the DGP is as follows. We generate $T\times 2$ matrix of factors from two independent AR(1) processes
For the factor series $\dot f_r\in\mathbb{R}^T$, we compute $f_r = \dot f_r/\|\dot f_r\|$ to ensure that the factor matrix satisfies Assumption (ref). All loading vectors are generated by taking norm-1 eigenvectors of a symmetric positive semi-definite matrix $A^\top A,$ where each element is $A_{i,j}\sim_{\rm i.i.d.} U(0,1).$
In this subsection, we first assess how changing sample sizes affects the estimation accuracy of factors and loadings, and show that the estimation improvements are perfectly aligned with the convergence rates given in Theorem (ref). We also show that the finite sample distribution of the estimation error is aligned with the asymptotic distribution shown in Theorem (ref).
For convergence rate, we simulate the 3-way tensor as described in Section (ref) where the parameters are set to be: (1) the AR(1) process $\dot{f}_r$ takes $\rho$ = 0.5 and $s_\varepsilon$ = 0.1; (2) the signal strength $d_1=2, d_2=1$ and the noise strength $s_u$ = 1; (4) the sample size of the baseline model is $(N,J,T)$=(30,30,30), and we compare the estimation error of the baseline model with that of the modified model. We consider three different cases of modification: (a) doubling sample sizes of all dimensions, $(N,J,T)$ = (60,60,60), (b) doubling the sample sizes of two dimensions, $(N,J,T)$ = (60,60,30), (c) doubling the sample size of only one dimension, $(N,J,T)$ = (60,30,30). We evaluate the estimates using the $\ell_2$ norm. Although the rank for the $M,F$ is 2, the convergence rates are the same for the first and the second rank and it is enough to examine the first rank. As the signs of $\hat\lambda_r, \hat\mu_r$, and $\hat f_r$ are undetermined, we calculate the errors as follows:
where $\mathrm{sign}(a)=\ensuremath{\mathds{1}}_{a>0} - \ensuremath{\mathds{1}}_{a<0}$.
Figure (ref) plots the histograms of the $\ell_2$ losses of the baseline versus modified DGPs. In panel (a) - (c), as we double the sizes of all dimensions, the estimation of the factor $\hat f_1$ and loadings $\hat \lambda_1,$ $\hat\mu_1$ all improve. The average error is reduced roughly by a half for the factor and two loadings vectors. This is aligned with the convergence rate in Theorem (ref), since doubling all $3$ dimensions of a tensor reduces the $\ell_2$ error by $1/2$. In panels (d) - (f), when we only double $N$ and $J,$ the improvement for the average error of $\hat\lambda_1$ is 0.01/0.015, while the improvement for $\hat\mu_1$ is 0.012/0.017. Both are roughly aligned with the reduction in the $\ell_2$ error by $1/\sqrt{2}$. On the other hand, the improvement for $\hat f_1$ is 0.0083/0.017, which is aligned with the reduction of the $\ell_2$ error by $1/2$. In panels (g) - (i), when we only double $N,$ there is no improvement for $\hat\lambda_1$ because $J$ and $T$ are unchanged in the $O_P(1/\sqrt{JT})$ rate; the improvement for both $\hat\mu_1$ and $\hat f_1$ are 0.012/0.017, which is aligned with the $1/\sqrt{2}$ improvement factor.
\setcounter{subfigure}{0}
For the asymptotic distribution, due to symmetry, we focus on inference for the loading vector $\mu_{j,r}\in\mathbb{R}$. We simulate the 3-way tensor using the same parameters as previous simulation except that we fix the sample size at $(N,J,T)=(30,30,30)$. Theorem (ref) shows that
which is what we aim to verify with simulations. Figure (ref) reports the histograms of the scaled estimation error $d_r s_u^{-1}\sqrt{NJT}(\hat{\mu}_{j,r} - \mu_{j,r})$ in 5000 MC simulations. For convenience, we only report the results of the first two elements of the loadings $\mu_{j,r}$ with $r=1,2$ and $j=1,2$. It can be seen that the finite sample distribution closely matches the asymptotic distribution, the QQ plot of the errors against the asymptotic distribution again confirms the finding. Interestingly we find that the bias appears to be negligible. Overall, these results show that the predictions of the asymptotic theory are valid in finite samples.
\setcounter{subfigure}{0}
\setcounter{subfigure}{0}
In this subsection, we show that ALS (Algorithm (ref)) has faster convergence rate than TPCA (Algorithm (ref)) under weak factor scenarios. According to Theorem (ref) and Theorem (ref), under weak factors, the TPCA rate is $O_P\left(\frac{1}{N^{0.6}} + \frac{1}{N^{0.2}}\right)$ while the ALS rate is $O_P\left(\frac{1}{N^{0.6}}\right)$ for a sufficiently large $N$, provided that the number of iterations is $\bar k$ is large enough. We follow the same DGP as Section (ref) except we set the signal strength $\sigma_1 = \sqrt{(NJT)^{1.6/3}}$ for weak factor. The sizes of all dimensions are set to be the same so that the finite sample properties are the same for all dimensions, and it is enough to examine only one dimension.
Figure (ref) (a) plots the comparison between TPCA and ALS when all dimensions are doubled from $(30,30,30)$ to $(60,60,60).$ The dotted line marks the mean of the $\ell_2$ errors. It can be noted that ALS shows more improvement in error reduction, from the mean of 0.16 to 0.13, while TPCA error is reduced from 0.21 to 0.18. The error reduction rate of ALS - 0.13/0.16 - is roughly in line with the theoretical value of $1/(2^{0.6})$, and the error reduction of TPCA - 0.18/0.21 - is also roughly aligned with the theoretical value of $1/(2^{0.2}),$ which is slightly slower under weak factors.
Figure (ref) (b) shows the improvement of ALS over TPCA under weak factors. It can be noted that ALS with 2 iterations already improves TPCA, and ALS with 10 iterations further improves the estimation, though by slightly less.
\setcounter{subfigure}{0}
We conclude with power properties of our maximum eigenvalue ratio test for the number of factors discussed in Section (ref). The DGP is generated as described in Section (ref). We generate rank $(1,2,2)$ model, and test the null hypothesis that the rank is 1 against the alternative that there are more than 1 but less than $K$ factors for each matricization.
The parameters are designed as follows: (1) the scale component is $\sigma_r$ = $d_r \times \sqrt{NJT}$ with $d_1=2$ and we gradually increase $d_2$ to study the power properties of the test, so when $d_2=0$, the empirical rejection probability corresponds to empirical size of the test, and none zero $d_2$ corresponds to empirical power for rank 2 dimensions, (2) we study cases where the idiosyncratic errors are generated with Gaussian distribution or Student's t distribution (the degrees of freedom are 5), and in both cases the variance of the errors is normalized to be $s_u=1,$ (3) we study finite sample properties of the test using 3-way tensor, and the sizes of the tensors is $(N,J,T) = 30\times40\times 50$, with the first dimension being the smallest.
We perform the test for the number of factors as discussed in Section (ref), where we approximate the asymptotic distribution of the statistics by randomly generating Gaussian matrices 5,000 times, and we also replicate the simulations for 5,000 times to calculate the empirical rejection probability for each scenario.
\setcounter{subfigure}{0}
\setcounter{subfigure}{0}
Figure (ref) reports the empirical rejections probabilities of the test for the null hypothesis of 1 factor against alternative of more than 1 factor but less than $K$ factors, with $K=3,5,7$, on a $3^{\text{rd}}$ order tensor. The test is performed on each matricization individually and plotted separately. The power curves show that the empirical size of test is very close to the nominal level of $5\%$,\footnote{Since the rank is 1 for the first dimension, the power stays at the nominal level.} and the empirical power of the test reaches $1$ as the strength of the second factor $d_2$ increases. The comparison among different $K$ values implies that, when the true number of factors is within range of the alternative hypothesis, the tighter the range of the alternative is, the more likely we reject the null when the alternative is true. This finding is true no matter which matricization we perform the test on.
We also make comparisons of the test for different types of errors. The two plots of Figure (ref) correspond to tensor with: (a) Gaussian errors, (b) Student's t distributed errors. Among different matricizations, the dimension with the larger size (`mat3' ) tend to climb faster to one than smaller size (`mat2'), meaning the matricization with the larger dimension gives the higher empirical power.
Finally, while our simulation results rely on the critical values that assume Gaussian errors, we learn from Figure (ref) (b) that the test performs well for non-Gaussian distributions albeit with some loss of power.
bryzgalova2022missing drew attention to the phenomenon of missing data in firm characteristics, which are the cornerstone of academic research in finance. They provide a comprehensive analysis of missing data in firm characteristics, and show that patterns of missing characteristics vary substantially across characteristics. As a remedy they propose a statistical factor model for imputing missing values and investigate the impact of missingness on asset returns.\footnote{Factor models can also be used to impute the missing counterfactual potential outcomes; see bai2021matrix.}
Firm characteristics data are of the tensor type. Indeed, for each firm, one has across time a set of observed or missing characteristics. The statistical methods proposed by bryzgalova2022missing do not exploit the tensor data structure, while in this empirical application we do and we show that using our proposed tensor data models improve on their time series of cross-sectional models.
We conduct the empirical study using the dataset created by freyberger2020dissecting. The data cover 35 firm characteristics, which are described in Table (ref). The dataset is monthly and ranges from January 1966 to December 2020, and there are a total of 13588 firms in the entire sample. This number is smaller than 22630 in bryzgalova2022missing, because we only have access to the pre-cleaned dataset that has no missingness in the cross-section of characteristics, i.e., all characteristics exist for all firms at any time period. However, this difference does not affect the empirical results because the RMSE are calculated based on simulated (masked) missing characteristics both in this paper and in bryzgalova2022missing. The results calculated based on cross-sectional factor model are also very close to those in bryzgalova2022missing. This tensor dataset is highly unbalanced where only $15.78\%$ of the data are not missing, and the average length of non-missing characteristics of all firms is 104.16 months.
Prior to estimating the tensor factor and the cross-sectional factor models, the raw characteristics are converted into rank quantiles within range $\left[ -0.5,0.5\right].$ This is done because: 1) different characteristics have different scales, and PCA has best performance when applied to data of similar scales, 2) centered rank-normalized characteristics are stationary in the cross-section and over time, and 3) it handles outliers. In practice, the predicted missing rank quantiles can be converted back to the original scale of the raw data.
Instead of fitting cross-sectional two-way factor models to each slice of the tensor (i.e., each matrix at period $t$), we fit three-way factor models on non-overlapping 5-year horizons. There is a considerable amount of parameter proliferation when estimating time series of cross-sections, thus we would expect a better performance from the parsimonious tensor-based approach. We compare our model to the traditional two-way approach in terms of imputation accuracy, and several other alternatives are also considered. See Section (ref) for detailed description of the model and other approaches being compared.
Similar to bryzgalova2022missing, we consider the standard evaluation metrics for model prediction errors - RMSE (root-mean-squared errors). The aggregate RMSE is averaged over all firms, characteristics, and time periods, which is calculated as follows:
We report both in-sample (IS) and out-of-sample (OOS) RMSE. For measuring out-of-sample RMSE, we randomly mask observed characteristics and calculate RMSE comparing predicted characteristics and the true masked ones. We take the same appraoch as bryzgalova2022missing to randomly mask $10\%$ of the characteristics using two schemes:
We also report the $R^2$ that measures the explained variation of the model relative to the total variation, which is a transformation of the RMSE.
Estimating both the tensor and cross-sectional factor models involves determining the number of factors. As noted in bryzgalova2022missing, the optimal number of factors for the cross-sectional models is between 6 and 8, and they select 6 factors as their final model. In comparison, we use the test introduced in Section (ref) to determine the number of factors for each dimension. The test is conducted with the following hypothesis:
where $K = k+1$ and we increase $k$ from $1$ to $20$ ($20$ is the upper bound of the number of factors we deem reasonable). Since the test is conducted multiple times on the 5-year sub-samples, we use Bonferroni correction to reduce the probability of type I errors.
Figure (ref) plots the p-values of test where the x-axis is test parameter $k$, and for each $k$ there are 11 dots representing the p-value of each 5-year sub-sample. The test is conducted on each dimension and for each missing scheme - missing-completely-at-random (MCAR) and missing in blocks of one year (Block). The red dash line is the Bonferroni adjusted significance level. The test result suggests that there is large number of factors in the characteristic dataset (i.e., there is no sharp change in the eigenvalues detected). Therefore, we select 20 factors as the baseline model for optimal results\footnote{The test result is aligned with the empirical finding that, when the number of factors is smaller than 20, the more factors the better the imputation performance. We also considered several alternative testing methods leading to the same conclusion.}, and we also present the results with 6 factors for a comparison with the cross-sectional model.
In this subsection, we compare all the methods listed in Table (ref)\footnote{BF-TPCA is not reported because it shows slight worse performance than B-TPCA, and B-TPCA is the considered the best model among all.}. Some methods require backward and/or forward information, and they become inapplicable when this information is not available. Hence, we use the natural fallback approach, where we replace a method that requires time-series information by a pure lower-level method. For example, B-TPCA is replaced by TPCA and B-XS-ridge is replaced by XS-ridge when the previous characteristic is not observed.
Table (ref) summarizes the main imputation results of different methods. We report the in-sample, out-of-sample missing completely at random (OOS MCAR), out-of-sample block missing (OOS Block) for all characteristics using the full sample period. The cross-sectional (XS-ridge) based methods and the Tensor PCA (TPCA) based method are estimated with different number of factors $R$ with $R\in \{6, 20\}$. The first thing to note is that the cross-sectional median method produces imputation errors that are more than twice as large compared to previous value and autoregression methods, both in-sample and OOS missing completely at random. This demonstrates the importance of considering time-series dependency among firm characteristics. The previous value and autoregression methods are relatively reliable when the characteristics are not missing in blocks. However, these two method are no longer applicable when previous values are not available under the case of block missing, which is very common in firm characteristics. Therefore, the results of previous value and autoregression are very close to that of cross-sectional median. Without the combination of backward/forward information, the previous value and autoregression methods are better than XS-ridge in both in-sample and OOS MCAR, but are a lot worse when characteristics are missing in blocks. However, TPCA is consistently better than the standard approaches in both out-of-sample performances.
The results comparing TPCA and XS-ridge methods indicate that XS-ridge method is overfitting the characteristics data. XS-ridge method has lower RMSE and higher $R^2$ in-sample than TPCA with the same number of factors, while it falls behind TPCA out-of-sample under both missing patterns. The differences are most prominent with 20 factors and under OOS MCAR, where B-XS-ridge is about $40\%$ higher than B-TPCA in terms of RMSE and about $10\%$ lower in terms of $R^2.$ The differences are smaller under block missingness because of lacking in time series information. However, the difference still demonstrates that tensor factor model is superior to cross-sectional factor model by incorporating the time dimension.
TPCA remains superior to cross-sectional factor model when combined with backward and forward information. Since TPCA more accurately predicts the missing characteristics than XS-ridge, the superiority is preserved after being combined with the same information. In addition, ALS further improves the imputation accuracy by providing a faster convergence rate.
Modern datasets are often multidimensional, extending beyond the $2$-dimensional panel data structures used in traditional factor models and PCA. In this paper, we study a class of $d$-way factor models for high-dimensional tensor data, which are a natural generalization of traditional $2$-way factor models. We demonstrate that $d$-way factor models can be estimated with a variation of PCA, which we call TPCA. This simple algorithm is optimal for the strong factor model. We also consider an improved iterative ALS algorithm which is optimal when the factors are moderately weak.
Additionally, we propose the first formal statistical test for the number of factors in a tensor factor model. Our findings indicate that the tensor factor model offers efficient dimensionality reduction compared to naively pooled traditional factor models. Simultaneously, the model is parsimonious with easily identifiable factors and loadings. These conclusions are supported by extensive simulation results. Lastly, we consider an empirical application to imputing missing firm characteristics.
Interesting applications of TPCA and our results could potentially include more refined panel data models with covariates, e.g., see freeman2022multidimensional and beyhum2020factor as well as causal inference and imputations with tensor data; see agarwal2020synthetic.