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.
104,410 characters · 16 sections · 124 citation commands
Rank Determination in Tensor Factor Model
Factor models have become a popular dimensional reduction tool in economics and statistics, especially for analyzing high dimensional time series. In practice a few common factors can often capture a large amount of variations and dynamics among a large pool of variables and time series. In the finance literature, chamberlain1983 exploited factor analysis to extend classical arbitrage pricing theory. In macroeconomics, bai2002,bai2003,stock2002 considered static factor models for modeling macroeconomic time series. forni2005 studied the identification of economy-wide and global shocks using generalized dynamic factor models. fan2011, fan2013, fan2019 established large covariance matrix estimation based on the static factor model. Factor models are also used to evaluate the impacts of various policies; see, e.g., bai2014, ouyang2015 and li2017estimation. Recently, large matrix or tensor (multi-dimensional array) data has become ubiquitous. wang2019 proposed a matrix factor model and applied it to matrix-valued financial data. chen2021factor analyzed the multi-category import-export network data via tensor factor model.
A critical step in building a factor model is to correctly specify the number of factors used in the model. Estimation and forecasting procedures are all depended on the number of factors. Moreover, in some cases the number of factors may have some crucial economic interpretations and important theoretical consequences. For example, in finance and macroeconomics, it provides the number of sources of nondiversifiable risk or the fundamental shocks driving the macroeconomic dynamics. See, forni1998let, stock2016dynamic, giannone2006, forni2009, among others.
Over the past decades, many methods have been developed to determine the number of common factors needed for modelling high dimensional vector time series. The most widely studied approach is to utilize the behavior of the eigenvalues of the covariance matrix (see, e.g., bai2002), or the singular values of the autocovariance matrix (see, e.g., lam2012). By the definition of factor models, the eigenvalues or the singular values corresponding to the systematic components must increase with the number of cross-sectional units. The rest of the eigenvalues, which represents idiosyncratic components, stay bounded or remain to be zero. In the static factor model, bai2002 proposed to estimate the number of factors by separating diverging eigenvalues from the rest using threshold functions, in the form of an information criterion. Alternative criteria based on random matrix theory have been studied in kapetanios2010 and onatski2010 for the static factor model. Specifically kapetanios2010 developed sequential tests and employed a subsampling method to obtain an approximation of the asymptotic distribution of the estimated eigenvalues. onatski2010 constructed tests based on the empirical distribution of the eigenvalues. In addition, onatski2012 proposed an alternative estimator using {the difference of consecutive eigenvalues}. bai2007 and amengual2007 extended the work of bai2002 to the restricted dynamic factor model. hallin2007 further extended the framework to the generalized dynamic factor model through thresholding eigenvalues of the spectral density matrix. They also proposed a data-dependent method to adjust the multiplicative constant of the penalty function. alessi2010improved introduced a tuning multiplicative constant in the penalty for dealing with approximate factor models. kong2017 employed similar ideas to study continuous time factor model with high frequency data. li2017 modified bai2002's procedure to the case that the number of factors is allowed to increase with the sample size. trapani2018randomized proposed a randomized sequential test procedure to determine the number of factors.
An alternative approach is to study the ratio of each pair of adjacent eigenvalues, with the insight that ratio of the smallest eigenvalue among these corresponding to the system component and the largest eigenvalue among these corresponding to the idiosyncratic component goes to infinity. Under stationary conditions, ahn2013 developed such an estimator based on the sample covariance matrix. lam2012 used such a ratio based estimator based on singular values of the aucovariance matrix, under an alternative definition of factor models proposed in pan2008 and lam2011.
Other than the eigenvalue-based methods, ye2003 developed an eigenvector based order determination procedure. luo2016 proposed a new estimator that combines both the eigenvalues and the bootstrap eigenvector variability. jung2018 suggested to sequentially test skewness of the squared lengths of residual scores that are obtained by removing leading principal components. However, these works assumed that the data are temporally independent, which are unlikely to hold for economic data.
These studies all focus on panel (vector) time series. Recently there is a growing interest in analyzing matrix- or tensor-valued time series, as such time series is encountered more and more frequently in applications, including Fama-French 10 by 10 series wang2019, a set of economic indicator series among a set of countries chen2021autoregressive, multi-category international trading volume series Hoff2011,chen2019matrix, multi-type international action counts among a group of countries Hoff2015, sequence of realized covariance matrices lunde2016econometric, kim2019factor, sequence of gray-scale face recognition images chen2021statistical, dynamic networks barabasi1999emergence,jiang2020autoregressive, dynamic human brain transcriptome data liu2022characterizing, multivariate spatial-temporal climate series chen2020semiparametric, neuroimaging data zhang2019cross,zhou2013. Factor model is again developed as an effective dimension reduction tool wang2019,chen2021factor,han2020. Same as for the vector factor models, it is important to determine the number of factors in these models.
In this paper, we consider the determination of the dimension of the core tensor factor in the tensor factor model in chen2021factor and han2020, which assumes the form \[ {\cal X}_t={\cal M}_t+{\cal E}_t={\cal F}_t\times_1 A_1\times_2\ldots\times_K A_K+{\cal E}_t. \] Similar to lam2012, the noise tensor ${\cal E}_t$ is assumed to be a white tensor process with potentially strong contemporary correlations among the elements of the noise tensor, and all common dynamics is absorbed in the signal process ${\cal M}_t$. This model setting is different from the approximate factor model in bai2002 and the dynamic factor model in hallin2007, in which the noise process is allowed to have weak auto-correlations, but with strong restriction on the contemporary correlation.
chen2021factor and han2020 studied the estimation procedures of the tensor factor model, assuming the ranks of the core tensor ${\cal F}_t$ is given, with some ad hoc rank determination suggestions. In this paper we formally propose two criteria for specifying the ranks of the core factor process, which we name “the information criterion” (IC) and “the eigenvalue ratio” (ER). They are all based on examining the eigenvalues of the sample cross-auto-moment of the observed tensor time series, utilizing the whiteness property of the noise process. The IC estimators aim at truncating eigenvalues, which is similar to the information criteria in vector factor models (e.g., bai2002 and hallin2007). The ER estimators are obtained by minimizing the ratio of two adjacent eigenvalues arranged in ascending order, extending the standard ER estimator in lam2012 and wang2019 with an added small penalty term in both the numerator and denominator of the ratio. The penalty term behaves like a lower bound correction to the true zero eigenvalues. We adopt similar ideas of the TOPUP and TIPUP procedures of chen2021factor, and their corresponding iterative versions, iTOPUP and iTIPUP, of han2020, to construct sample auto-cross-moments. Our theoretical and empirical investigations show that estimators based on the iterative algorithms are much better than that based on the non-iterative ones, as the iterative algorithms significantly improve the estimation accuracy of the eigenvalues. The finite sample properties of the IC and ER criteria are also good. The empirical evidences show that the best estimators in tensor factor model are the IC and ER estimators based on iTIPUP, under some mild conditions on the level of signal cancellation typically associated with the TIPUP based procedures.
This paper is organized as follows. Section (ref) briefly describes the tensor factor model and the corresponding estimation procedures proposed in chen2021factor and han2020. Section (ref) introduces the criteria for determining the ranks of the core tensor factor process, and their iterative versions. Section (ref) investigates theoretical properties of the proposed estimators. Section (ref) presents simulation studies of the finite sample properties of the proposed methods. Real data analysis is given in Section (ref). Discussions are provided in Section (ref). All technical details are relegated to the Supplementary Material.
In this section, we first propose two general order determination criteria based on the properties of the estimated eigenvalues of a semipositive matrix.
Let $\widehat W$ be a $p\times p$ symmetric and non-negative definite matrix, which is a sample version of a true $p\times p$ symmetric and non-negative definite matrix. We assume $W=\mathbb{E} \widehat W$. Also let $\hat \lambda_j$ be the eigenvalues of $\widehat W$ such that $\hat\lambda_1\ge \hat\lambda_2\ge\ldots \ge \hat\lambda_{p}$. Let $\lambda_1\ge \ldots\ge \lambda_r> \lambda_{r+1}=\ldots =\lambda_{p}=0$ be the eigenvalues of $W$. Note that the rank of $W$ is $r$.
Let $m^*<p $ be a predefined upper bound and functions $G(\widehat{W})$ and $H(\widehat{W})$ be some appropriate positive penalty functions. We propose the following two quantities
The first criterion in ((ref)) is similar to an information criterion as its first term mimics the residual sum of squares of using a rank $m$ matrix to approximate the matrix $\widehat W$ while the second term $mG(\widehat W)$ penalizes the model complexity $m$. We will call it {\it the information criterion (IC)}. The second criterion in ((ref)) uses the ratio of two adjacent eigenvalues of $\widehat W$, with a small penalty term $H(\widehat W)$ added to both the numerator and denominator. We will call it {\it the eigen-ratio criterion (ER)}.
Here we show a consistency result for the general estimator. To be more precise, let $\widehat W^{(n)}$ and $W^{(n)}$ be two sequences of semi-positive symmetric matrices, with $W^{(n)}=\mathbb{E} \widehat W^{(n)}$ and $n$ be an index associated with the sample size and dimension. Also let $\hat \lambda_j^{(n)}$ be the eigenvalues of $\widehat W^{(n)}$ such that $\hat\lambda_1^{(n)}\ge \hat\lambda_2^{(n)}\ge\ldots \ge \hat\lambda_{p}^{(n)}$. Let $\lambda_1^{(n)}\ge \ldots\ge \lambda_r^{(n)}> \lambda_{r+1}^{(n)}=\ldots =\lambda_{p}^{(n)}=0$ be the eigenvalues of $W$. We assume $\lambda_r^{(n)}\to\infty$ as $n\to\infty$. The following proposition provides the sufficient conditions for the consistency of the IC and ER estimators in (ref) and (ref), respectively. It provides a guideline of choosing proper penalty functions $G(\cdot)$ and $H(\cdot)$ in order determination for any generic $\widehat W^{(n)}$.
In the conditions of Proposition (ref), $\gamma_n$ represents the convergence rate of the sample eigenvalues corresponding to the non-zero eigenvalues of $W^{(n)}$, and $\beta_n$ represents the rate of the sample eigenvalues corresponding to the zero eigenvalues of $W^{(n)}$. For example, under the strong factor model of lam2012's setting, $\gamma_n=p^2T^{-1/2}$ and $\beta_n=p^2T^{-1}$, where $T$ is the sample size.
For the specific problems such as the tensor factor model problem we focus there, a detailed analysis of the rates $\gamma_n$ and $\beta_n$ is needed to construct the penalty functions $G(\cdot)$ and $H(\cdot)$ and to establish the consistency of the rank estimators. In fact one can establish the convergence rate with a more detailed analysis beyond the simple consistency results in Proposition (ref), as we will do for the tensor factor model.
Here we briefly introduce the tensor factor model setup in chen2021factor and han2020. A tensor factor model can be written as
where ${\cal X}_t\in\mathbb{R}^{d_1\times\cdots\times d_K}$ is the observed tensor at time $t$, the core tensor ${\cal F}_t$ is the unobserved latent tensor factor process of dimension $r_1\times\ldots \times r_K$, $A_k$ are the deterministic loading matrix of size $d_k\times r_k$ and $r_k\ll d_k$, and ${\cal E}_t$ is the idiosyncratic noise components of ${\cal X}_t$, which is assumed to be a white process. Here the $k$-mode product of ${\cal X}\in\mathbb{R}^{d_1\times d_2\times \cdots \times d_K}$ with a matrix $U\in\mathbb{R}^{d_k'\times d_k}$, denoted as ${\cal X}\times_k U$, is an order $K$-tensor of size $d_1\times \cdots \times d_{k-1} \times d_k'\times d_{k+1}\times \cdots \times d_K$ such that $$ ({\cal X}\times_k U)_{i_1,...,i_{k-1},j,i_{k+1},...,i_K}=\sum_{i_k=1}^{d_k} {\cal X}_{i_1,i_2,...,i_K} U_{j,i_k}.$$ The core tensor ${\cal F}_t$ is usually much smaller than ${\cal X}_t$ in dimension. We also assume that the rank of $A_k$ is $r_k$. Otherwise ${\cal X}_t$ in (ref) may be expressed equivalently with a lower-dimensional factor process. The parameters $r_1,...,r_K$ are assumed to be fixed but unknown. For more details of the tensor factor model (ref), see chen2021factor and han2020.
It is obvious that the loading matrices $A_k$ are not identifiable in Model (ref). Model (ref) is unchanged if we replace $(A_1,...,A_K, {\cal F}_t)$ by $(A_1H_1,...,A_KH_K, {\cal F}_t\times_{k=1}^K H_k^{-1})$ for any invertible $r_k\times r_k$ matrix $H_k$. However, the linear space spanned by the columns of $A_k$, called the factor loading space, is uniquely defined. Assume $A_k$ has a SVD representation $A_k=U_k\Lambda_k V_k^\top$. Then, the factor loading space of $A_k$ can be represented by the orthogonal projection $P_k$,
The two criteria introduced in Section (ref) can be used to estimate the number of factors in the tensor factor model ((ref)), using properly constructed matrices $W$ and $\widehat W$. Particularly, we will study the following four constructions.
Since they have been proposed and used for loading space estimation in chen2021factor and han2020, we adopt the same names to represent them.
{\bf (I) TOPUP:\ } Let
and
where $\otimes$ is the tensor product such that, for any ${\cal A}\in\mathbb{R}^{m_1\times m_2\times \cdots \times m_K}$ and ${\cal B}\in \mathbb{R}^{r_1\times r_2\times \cdots \times r_N}$, $$({\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} ,$$ and ${\rm mat}_k$ is the tensor unfolding (into a matrix) operation along mode-$k$ of a tensor. Here we emphasize that $\widehat W_k$ is constructed using ${\cal X}_{1:T}=({\cal X}_1,\ldots,{\cal X}_T)$. The constant $h_0$ is a (small) predetermined integer and the sum over $h$ in $\widehat W_k$ is to accumulate the information from different time lags $h$. The rank of its population version can be shown to be $r_k$ under certain conditions, hence we can use the IC and ER estimators presented in Section (ref) to determine $r_k$.
{\bf (II) TIPUP: \ } Define a $d_k\times (d_k h_0)$ matrix as
which replaces the tensor product in $\operatorname*{\text{mat}_1}({\rm{TOPUP}}_k({\cal X}_{1:T})$ by the inner product. Let
{\bf (III and IV) iTOPUP and iTIPUP: \ } han2020 proposed an iterative procedure to estimate $U_k, k=1,...,K$ in (ref), based on either TOPUP or TIPUP procedure. Briefly, at $i$-th iteration, suppose we have obtained an estimate of the ranks $\widehat{r}_k^{(i-1)}$ ($k=1,\ldots,K$) and their corresponding $\widehat U_{k,\widehat r_k^{(i-1)}}^{(i-1)}$ at $(i-1)$-th iteration, we calculate the orthogonal projections of ${\cal X}_t, 1\le t\le T$, to obtain
Note that ${\cal Z}_{k,t}^{(i)}$ uses projection of ${\cal X}_t$ on all modes, except mode-$k$. The initial ranks $\widehat{r}_k^{(0)}$ ($k=1,\ldots,K$) and their corresponding $\widehat U_{k,\widehat r_k^{(0)}}^{(0)}$ can be obtained through the non-iterative TOPUP and TIPUP procedure. Let ${{\cal Z}}_{k,1:T}^{(i)}=({{\cal Z}}_{k,1}^{(i)},\ldots,{{\cal Z}}_{k,T}^{(i)})$, and define
The iterative procedure is motivated by the observation that ${\cal Z}_{k,t}^{(j)}$ is a $r_1\ldots r_{k-1}d_k r_{k+1}\ldots r_K$ tensor, much smaller than ${\cal X}_t$, which is $d_1\ldots d_k$ tensor. Hence $A_k$ can be estimated more accurately if all $A_j(j\neq k)$ are given in advance or can be estimated accurately, since the convergence rate now depends on $(d_k/r_k)\prod_{i=1}^K r_i$ rather than $\prod_{i=1}^K d_i$.
The IC and ER estimators are constructed by replacing $\widehat W$ in ((ref)) and ((ref)) with $\widehat W_k$, $\widehat W_k^*$, $\widehat W_k^{(i)}$ and $\widehat W_k^{*(i)}$. This yields eight different criteria, summarized in Table (ref). Again, we use the same names of the procedures as that in chen2021factor and han2020 to represent the various constructions of $\widehat W$. For the iterative procedures, we start with an initial rank estimate $r_k^{(0)}$, $k=1,\ldots,K$ and estimate the ranks through iteration until convergence. See remark below for setting the initial ranks and the stopping criteria.
{\bf The choice of the penalty function $G(\cdot)$ and $H(\cdot)$: \ } Both criteria essentially try to distinguish the smallest (true) non-zero eigenvalue from the true zero eigenvalue using noisy estimators of the eigenvalues. Hence the penalty function is closely related to the amount of error in the eigenvalue estimation and the strength of the smallest (true) non-zero eigenvalue. We consider the following penalty functions $G(\cdot)=g_k(d,T)$:
where $d=\Pi_{k=1}^K d_k$ and $\nu$ is a tuning parameter. Ideally $\nu$ should be chosen to be the strength of the weakest factor (see Assumption (ref) in Section (ref)), though in practice we usually do not know its precise value. A more thorough discussion on this issue will be given later in Remark (ref). Note that only $g_{k,5}$ involves $k$.
For the eigen-ratio criterion, we consider the following penalty function $H(\cdot)=h_{k}(d,T)$:
where $c_0$ is a small constant, e.g. $c_0=0.1$. Note that the penalty functions scale with $h_0$, because the strength of divergent eigenvalues increases with $h_0$. Our theoretical analysis indicates that a better penalty function $H(\cdot)$ should also involve the strengths of the factors, similar to $G(\cdot)$. However, the function $G(\cdot)$ has a much wider allowable range, and in most of the situations a simple constant function $h_{k,1}$ is sufficient. A more detailed discussion will be given later in Remark (ref).
{\bf More considerations of the iterative procedure:} The non-iterative procedures estimate $r_k$ ($k=1,\ldots,K$) individually. The accuracy of $\hat{r}_k$ does not depend on the accuracy of the estimation of the ranks in other directions. On the other hand, the iterative procedures estimate all the ranks simultaneously, hence the accuracy of estimated rank in one direction depends on that in all other directions. The iterative algorithm improves the estimation accuracy of the eigenvalues and the principal subspace because the projected tensor ${\cal Z}_{k,t}^{(i)}$ in (ref) is of lower dimensional than ${\cal X}_t$. Figure (ref) numerically shows that the iterative algorithms improve the accuracy of estimated true zero eigenvalues over the non-iterative algorithms. In iterative algorithms, one would need to specify all the ranks $\widehat{r}_k^{(i)}$ ($k=1,\ldots, K$) in each iteration. Intuitively, an overestimated $\widehat{r}_k^{(i)}>r_k$ would still produce consistent estimators, since the non-iterative procedure, using $\widehat{r}_k^{(i)}=d_k$ for all other directions, is consistent. Our theoretical results shown later confirm that, if $\widehat r_k^{(i-1)}$ used is larger than the true $r_k$, the iterative algorithm warrant the consistency of the IC and ER estimators at $i$-th iteration. However, an underestimated $\widehat{r}_k^{(i)}<r_k$ would potentially result in loss of signal strength hence negatively impacting the estimation in other dimensions. A precise quantification of the impact requires a more detailed investigation. But the numerical studies show that for iteration $i>1$, the performance is relative robust by using the order obtained by the IC or ER criteria. This is partially due to the theoretical justification that, for fixed ranks $r_k$, iterative algorithm only needs one iteration to achieve the ideal convergence rate for the estimation of the eigenvalues (see Theorem (ref) later). In fact, if one has a priori information about a possible maximum (fixed) ranks of the core factor process, one could use such ranks in the iterative algorithms accordingly.
For iTOPUP procedure, we use $\widehat{r}_k^{(i)}=\widehat r_k^{(i)}({\rm IC})$ or $\widehat{r}_k^{(i)}=\widehat r_k^{(i)}({\rm ER})$ after the initial iteration $i\geq 1$. However, one needs to use a more conservative estimator of the rank for the initial step since $\widehat r_k^{(0)}({\rm IC})$ and $\widehat r_k^{(0)}({\rm ER})$ tend to be inaccurate. We suggest to use $\hat{r}_k^{(0)}=\min\{2\widehat r_k^{(0)}({\rm IC}), \widehat r_k^{(0)}({\rm IC})+3\}$ or $\hat{r}_k^{(0)}=\min\{2\widehat r_k^{(0)}({\rm ER}), \widehat r_k^{(0)}({\rm ER})+3\}$ by default, unless one has prior knowledge of the number of the factors. iTIPUP procedure is similar. Although it is safer to use larger initial ranks $r_k^{(0)}$, it is often not necessary to be extremely conservative, as the initial loss of signal strength of using a rank too small can be corrected later through iterations.
In the iterative algorithms, iteration is not stopped until the convergence of both the rank estimators and the loading space estimators. Theoretical properties in Section (ref) only state consistency results in each iteration step. This stopping rule is mainly suggested by simulation study. In practice, we may stop the algorithm when the estimated number of factors in current iteration is the same as that in previous iteration.
We introduce some notations first. Let $d=\prod_{k=1}^K d_k$ and $d_{-k}=d/d_k$. For a matrix $A = (a_{ij})\in \mathbb{R}^{m\times n}$, write the SVD as $A=U\Sigma V^\top$, where $\Sigma=\text{diag}(\sigma_1(A), \sigma_2(A), ..., \sigma_{\min\{m,n\}}(A))$, with the singular values $\sigma_1(A)\ge\sigma_2(A)\ge \cdots\ge \sigma_{\min\{m,n\}}(A)\ge 0$ in descending order. The matrix Frobenius norm can be denoted as $\|A\|_{\rm F} = (\sum_{ij} a_{ij}^2)^{1/2}=(\sum_{i=1}^{\min\{m,n\}}\sigma_i^2(A))^{1/2}$. Define the spectral norm $$ \|A\|_{\rm 2} = \max_{\|x\|_2=1,\|y\|_2= 1} \|x^\top A y\|_2=\sigma_1(A).$$ The tensor Hilbert Schmidt norm for a tensor ${\cal A}\in\mathbb{R}^{m_1\times m_2\times \cdots \times m_K}$ is defined as $$ \|{\cal A}\|_{{\rm HS}}=\sqrt{\sum_{i_1=1}^{m_1}\cdots\sum_{i_K=1}^{m_K}({\cal A})_{i_1,...,i_K}^2 }. $$ Define the tensor operator norm for an order-4 tensor ${\cal A}\in\mathbb{R}^{m_1\times m_2\times m_3\times m_4}$, $$ \| {\cal A}\|_{\rm{op}} =\max\left\{ \sum_{i_1,i_2,i_3,i_4} u_{i_1,i_2} \cdot u_{i_3,i_4}\cdot ({\cal A})_{i_1,i_2,i_3,i_4}:\|U_1\|_{\rm F}=\|U_2\|_{\rm F}=1 \right\},$$ where $U_1=(u_{i_1,i_2})\in\mathbb{R}^{m_1\times m_2}$ and $U_2=(u_{i_3,i_4})\in\mathbb{R}^{m_3\times m_4}$. Define order-4 tensors
where $\otimes$ is the tensor product and $U_k$ is from the SVD form of $A_k=U_k \Lambda_k V_k^\top$. We view $\Phi^{(\text{\footnotesize cano})}_{k,h}$ as the canonical version of the auto-covariance of the factor process. Similarly define
Write $\Phi_{k,1:h_0}=(\Phi_{k,h}, h=1,\ldots,h_0)$ and $\Phi_{k,1:h_0}^*=(\Phi_{k,h}^*, h=1,\ldots,h_0)$. Denote $\overline\mathbb{E}(\cdot)=\mathbb{E}(\cdot|\{ {\cal F}_1,...,{\cal F}_T\})$. Let $\tau_{k,m}$ be the $m$-th largest singular value of $\overline\mathbb{E} ({\text{TOPUP}}_k ({\cal X}_{1:T}))$, $$\tau_{k,m}=\sigma_{m}(\overline\mathbb{E} ({\text{TOPUP}}_k ({\cal X}_{1:T})))= \sigma_{m}\big({\text{mat}}_1(\Theta_{k,1:h_0})\big).$$ Similarly, let $$\tau_{k,m}^*=\sigma_{m}(\overline\mathbb{E} ({\text{TIPUP}}_k ({\cal X}_{1:T})))= \sigma_{m}\big({\text{mat}}_1(\Theta_{k,1:h_0}^*)\big).$$ For simplicity, we write $U_k=U_{k,r_k}$ and $\widehat U_k=\widehat U_{k,r_k}$.
To facilitate consistency properties of the proposed procedures, we impose the following assumptions.
Assumption (ref) is the same assumption used in chen2021factor and han2020. This assumption corresponds to the white noise assumption of lam2011,lam2012. It allows substantial contemporaneous correlation among the entries of ${\cal E}_t$. Note that the normality assumption, which ensures fast convergence rates in our analysis, is imposed for technical convenience. In fact we only need to impose the sub-Gaussian condition. Assumption (ref) allows a very general class of time series models, including causal ARMA processes with continuously distributed innovations; see also tong1990non, bradley2005, tsay2005analysis, fan2008nonlinear, rosenblatt2012markov, tsay2018nonlinear, among others. The restriction $\theta_1\le 1$ is introduced only for presentation convenience. Assumption (ref) requires that the tail probability of any orthonormal projection of ${\cal F}_t$ decay exponentially fast. In particular, when $\theta_2=2$, ${\cal F}_t$ is sub-Gaussian.
Assumption (ref) is similar to the signal strength condition of lam2012, and the pervasive condition on the factor loadings (e.g., stock2002 and bai2003). It plays a key role in identifying the common factors and idiosyncratic noises in (ref). Indices $\delta_0,\delta_1$ are measures of the strength of factors, or the rate of signal strength growth as the dimension $d_k$ grows. When $\delta_0=\delta_1=0$, the factors are called strong factors; otherwise, the factors are called weak factors. In particular, $\delta_0$ represents the strength of the strongest factors and $\delta_1$ the strength of the weakest factors.
Instead of Assumptions (ref) to (ref), chen2021factor and han2020 imposed conditions on $\|\Theta_{k,0}\|_{\text{op}},\|\Theta_{k,0}^*\|_{\rm 2},\tau_{k,r_k}$ and $\tau_{k,r_k}^*$ in order to allow $r_k$ to increase with $d_k$. The following proposition establishes a connection between these two types of assumptions when the rank $r_k$ is fixed.
In this section, we shall present theoretical properties of the IC and ER estimators using non-iterative TOPUP, non-iterative TIPUP, iTOPUP and iTIPUP. We first introduce some quantities related to the estimation errors of the estimated eigenvalues of the four different methods. For non-iterative TOPUP, define
For non-iterative TIPUP, define
Similarly, we use $\tilde\beta_{k}$, $\tilde\gamma_{k}$ and $\tilde\beta_{k}^*$, $\tilde\gamma_{k}^*$ for iTOPUP and iTIPUP, respectively, where
It is clear that $\beta_k$, $\gamma_k$, $\beta_k^*$ and $\gamma_k^*$ dominate $\tilde\beta_k$, $\tilde\gamma_k$, $\tilde\beta_k^*$ and $\tilde\gamma_k^*$, respectively. Actually, $\beta_k$ and $\gamma_k$ (resp. $\beta_k^*, \gamma_k^*$, or $\tilde\beta_k, \tilde\gamma_k$, or $\tilde\beta_k^*, \tilde\gamma_k^*$) correspond to the $\beta_n$ and $\gamma_n$ sequence in Proposition (ref).
We impose the following set of conditions to ensure that the estimation error of the divergent eigenvalue is much smaller than the true smallest non-zero eigenvalue, and the estimation error of the zero eigenvalues is relatively small. The $\gamma$'s above are part of the estimation errors of the divergent eigenvalues, and the $\beta$'s are the estimation errors of the (true) zero eigenvalues, for the four different estimation method. Note that $d^{2-2\delta_1}$ is the growth rate of the smallest non-zero eigenvalues corresponding to the weakest factors. See also Theorem (ref) below.
We will impose the following sufficient conditions on the penalty function $g_k(\cdot)$ and $h_k(\cdot)$.
The following is a set of different sample size conditions for different settings. They ensure sufficiently large sample size $T$ so that the non-iterative (true rank) factor loading space estimator based on TOPUP or TIPUP is consistent (for (a) and (b)), or has a relatively small error (for (c) and (d)).
Theorem (ref) presents the asymptotic properties of the IC and ER estimators in (ref) and (ref) based on non-iterative TOPUP, non-iterative TIPUP, iTOPUP and iTIPUP.
In addition to the consistency of the rank estimators, we also have the following more detailed properties of the estimated eigenvalues. Let $\hat\lambda_{k,j}$ be the eigenvalues of $\widehat W_k$ (defined in Section (ref)) such that $\hat\lambda_{k,1}\ge \hat\lambda_{k,2}\ge...\ge\hat\lambda_{k,d_k}$, $1\le k\le K$. Also let $\lambda_{k,j}$ be the eigenvalues of population version $\mathbb{E}\widehat W_k$ such that $\lambda_{k,1}\ge ... \ge \lambda_{k,r_k}> \lambda_{k,r_k+1}=...=\lambda_{k,d_k}=0$. Similarly, define $\hat\lambda_{k,j}^*$, $\lambda_{k,j}^*$, $\hat\lambda_{k,j}^{(i)}$, $\lambda_{k,j}^{(i)}$, $\hat\lambda_{k,j}^{*(i)}$ and $\lambda_{k,j}^{*(i)}$ as the eigenvalues of $\widehat W_k^*$, $\mathbb{E}\widehat W_k^*$, $\widehat W_k^{(i)}$, $\mathbb{E}\widehat W_k^{(i)}$, $\widehat W_k^{*(i)}$, $\mathbb{E}\widehat W_k^{*(i)}$, $i\ge 1$, respectively.
Combining the discussion in Remarks (ref) and (ref), we can conclude that in general the iterative procedues are better than the non-iterative ones and the TIPUP based procedures are better than the TOPUP ones, assuming no signal cancellation when TIPUP is used (see Remark (ref)).
In this section, we compare the empirical performance of the proposed methods and their variants under various simulation setups. We consider the identification of the number of factors based on the non-iterative TIPUP and TOPUP methods (denoted as initial estimators), the one step iterative methods (denoted as one-step estimators) and the iterative procedures after convergence (denoted as final estimator). In the iterative algorithm, at $i$-th iteration ($i>1$), we use the order obtained by the IC or ER criteria. Iteration is stopped when both the rank estimators and the loading space estimators converge. We also check the performance of different choices of penalty function $g_k(d,T)$ and $h_k(d,T)$. Specifically, we consider penalty functions (ref) and (ref), and denote them as IC1-IC5, ER1-ER5, respectively. The empirical performance of IC1-IC5 (resp. ER1-ER5) are very similar, thus we only present IC2 and ER1 in this section. The detailed comparison are shown in Appendix B.
The simulation study consists of three parts. The first part is designed to investigate the overall performance of our methods and their comparisons under models with different factor strength. As the strength of the weakest factors is unknown, we by default set $\nu=0$ for all the penalty function $g_k(d,T)$ in (ref). In the second part, we investigate the case in which some factors have a dominantly strong explanatory power. The third part is the case in which we use $\nu=\delta_1$ in (ref) for all IC estimators, when some factors are weak. For each case, we compute the proportion of correct identification of the rank of the factor processes or the root mean squared errors (RMSEs) of the rank estimates from 1000 simulated data sets. In Section (ref), we study the selection of the optimal constant $c$ in IC criteria, using the method proposed in Remark (ref).
The simulation uses the following matrix factor model: \[ X_t=A_1 F_t A_2^\top+ E_t. \] Here, $E_t$ is white and is generated according to $E_t=\Psi_1^{1/2} Z_t\Psi_2^{1/2}$, where $\Psi_1, ~\Psi_2$ are the column and row covariance matrices with the diagonal elements being $1$ and all off diagonal elements being $0.2$. All of the elements in the $d_1\times d_2$ matrix $Z_t$ are i.i.d $N(0,1)$. This type of model of $E_t$ has been proposed and studied in the literature, see, for example Hoff2011,Hafner&2020,linton2020estimation. The entries $f_{ijt}$ in the factor matrix $F_t$ were drawn from independent univariate AR(1) model $f_{ijt}=\phi_{ij}f_{ij(t-1)}+\epsilon_{ijt}$ with standard $N(0,1)$ innovation.
In the first part, the following three models are studied:
All the factors in Models M1 are strong factors. Model M2 is the case in which four ($2\times 2$) factors are strong ($\delta_0=0$), twelve factors are weak factor with strength $0.2$ and the rest nine factors are weak factor with strength $\delta_1=0.4$. Model M3 is the case in which all the factors are very weak factors with strength $\delta_0=\delta_1=0.6$. Models M2 and M3 are designed to examine the effects of weak factors on the estimators. We choose a set of data dimensions to be $(d_1,d_2)=(20,20),(40,40),(80,80)$ and the sample size to be $T=100,300,500,1000$. Again, in this first part of simulation, we fix $h_0=1$ and set $\nu=0$ in the penalty function in (ref), under the assumption that all factors are strong, even though some of factors simulated are weak (e.g., true $\delta_1=0.4$ in Model M2).
For Model M1, the results in Table (ref) show clearly that, using TOPUP and IC, the initial estimator behaves very poorly even for large sample sizes. On the other hand, the one step estimator uniformly and significantly outperforms the non-iterative initial estimator. In addition, the final estimator performs the best over all choices of $d_1$, $d_2$ and $T$. We also observe that the performance improves as the dimension increases, except that $d_1=40$ is not as good as $d_1=20$ when $T=100$. This improvement is due to the fact that when with strong factors ($\delta_0=\delta_1=0$), larger dimension $(d_1,d_2)$ provides more data points and information on the rank $r_k$. With the same settings, IC criterion based on TIPUP determines the ranks perfectly, indicating that it is uniformly better than IC criterion based on TOPUP. This is partially due to the fact that non-iterative and iterative TIPUP procedures estimate the loading matrices and the eigenvalues more accurately than the corresponding TOPUP procedures. More interestingly, with the same setting, the ER estimators based on both TOPUP and TIPUP procedures perform perfectly.
For Model M2, Tables (ref) reports the proportion of correct rank identification using IC and ER estimators based on TOPUP and TIPUP procedures. It is seen that, for small sample sizes ($T=100$ and $300$), the performance of IC estimators deteriorate when we use iterative procedures, which may indicate that the sample size $T$ and dimension $d$ do not meet the required Assumption (ref)(a) in Theorem (ref). In addition, for $T=100$, the IC estimators using TIPUP procedures do not work at all, though they are better when $T=300$. This is due to the existence of weak factors and the fact that we use the default $\nu=0$ in (ref). When the sample size is small, the estimators tend to identify the strong factors while miss the weak factors as their corresponding eigenvalues are relatively small and comparable to the penalty function $g_k(\cdot)$. The IC estimators using TOPUP procedures performed much better for small sample sizes. The performance also becomes worse as $d_1$ and $d_2$ increases, also due to the existence of weak factors. The ER estimators based on TIPUP procedures show almost perfect accuracy, even there are weak factors in the model. We note that, even with weak factors, the performance remains almost the same with larger $(d_1,d_2)$. Again, accuracy improves by using the iterative procedure. The ER estimator based on TIPUP procedures are much better than all the other estimators for Model M2.
Table (ref) shows the more detailed identification results using IC2 and ER1 estimators for Model M2 over 1000 replications, based on TOPUP and TIPUP procedures. The sample size is $T=300$ and the data dimension is $(d_1,d_2)=(80,80)$. The true rank pair is (5,5) for Model M3. From the table, it is seen that the IC procedures tend to under-estimate the number of factors, with the correspoding iterative procedures perform the worst. The ER estimators using the TOPUP procedures are likely to pick up only the strong factors, as the gap between strong factors ($\delta_1=0$) and weak factors ($\delta_1=0.4$) may be larger than that between weak factors and true zero eigenvalue estimations. We note that the outstanding performance of the ER estimators using the TIPUP procedures is quite different from the performance of a similar ER estimator in vector factor models under similar mixed strong and week factor cases (see e.g. lam2012). The main reason is that the other tensor modes provide additional information and in certain sense serve as additional samples. Then, for each $k\le K$, the signals of all divergent eigenvalues depend on $d$ instead of $d_k$, leading to larger gap between weak factors and true zero eigenvalue estimations.
For Model M3 with all very weak factors ($\delta_0=\delta_1=0.6$), Table (ref) reports RMSEs of the ER1 estimators. The results of using IC estimators (not shown here) are significantly worse than that of the ER estimators due to the difficulty of IC estimators in dealing with weak factors when $\nu=0$ is used. It is seen from the table that ER estimators based on TIPUP procedure outperform that based on TOPUP procedure. We also see that the iterative algorithms improves the performance very significantly under this very weak factor case. The performance varies with the change of $(d_1,d_2)$ in a non-standard way, as the performance with $(40,40)$ seems to be better than that with $(20,20)$ and $(80,80)$. We note that in weak factor cases, a large $d_k$ will potentially reduce the accuracy of the estimator of $r_k$, since the signal level on the $k$-th dimension becomes weaker. On the other hand, a larger $d_i$ ($i\neq k$) in the other dimension potentially improves the estimation of $r_k$, since we have more “repeated” observations to be used for estimating $r_k$.
Table (ref) shows the relative frequency of different estimated ranks of the ER1 estimator based on both TIPUP and TOPUP procedures, for the case of $T=300$ and $(d_1,d_2)=(80,80)$ under Model M3 with the true rank $(5,5)$. It is seen that the ER estimators tend to overestimate the number of factors, when all the factors are weak. All ER1-TOPUP estimators essentially identify $(6,6)$ as the rank. The non-iterative ER1 estimator based on TIPUP procedure can overestimate the ranks by a large margin. However, the iterations can gradually correct the over-estimation.
In summary, the first part of simulation shows that the ER estimators and IC estimator based on TIPUP procedures perform very well when all the factors are strong. The ER estimators significantly outperform the IC estimators when some or all factors are weak. Different from the results shown in lam2012 for vector factor models with both strong and weak factors, in tensor factor models, the ER estimators (based on TIPUP procedure) are able to determine the correct number of factors in many cases. The results also show that the iterative procedure significantly improves the performance except the IC estimators based on TOPUP in Model M2, which may due to that the sample size $T$ and dimension $d$ do not meet the required Assumption (ref)(a) in Theorem (ref). It also shows that the estimators based on TIPUP perform better than that based on TOPUP in general. Hence, when some factors are weak, the iterative ER estimators based on TIPUP are the choice.
The second part of our simulation examines the effects of dominate strong factors on the IC and ER estimators. The data are generated from the following model,
We fix $d_1=d_2=40$ and $T = 200$. Again, we use IC2 and ER1 for demonstration, and assume $\nu=0$ in the penalty function in (ref). Although all of the four factors are strong factors, the strongly imbalanced signal strength in ${\cal F}_t$ makes one of factors dominating the others in explanatory power. Table (ref) reports the relative frequencies of estimated rank pairs over 1000 replications. It is seem that the ER estimators are very likely to pick up only the dominate factor, although iterations significantly improves the accuracy. The IC estimators performs much better in this case. And over all, estimators based on TIPUP perform better than the corresponding estimators using TOPUP. Overall, IC-TIPUP performs the best in this case.
The third part of the simulation considers the impact of penalty function selection. We consider Model M2 and M3 again, with $d_1=d_2=40$, but use the true weakest factor strength $\delta_1$ as $\nu$ in the penalty function $g_k(T,d)$ in (ref) instead assuming strong factor and use $\nu=0$ as in Section 5.1. To compare with Tables (ref) for M2 and (ref) for M3, we report the proportion of correct rank identification using IC2 estimators in Table (ref) for M2 and RMSEs of the IC2 estimators in Table (ref) for M3.
Comparing Tables (ref) and (ref), it is seen that the performance of IC2 estimators using TIPUP improve by using the correct $\nu$ in the penalty function. We notice that the initial and one step IC estimators using TOPUP with the correct $\nu$ in the penalty term actually under-perform the ones with the incorrect $\nu$. The reason for this unusual behavior of TOPUP is unclear. It might be that the magnitude of the estimation of true zero eigenvalue is much larger in this mixed weak and strong factor case, as shown in Figure (ref). Hence, reducing the penalty term by using $\nu=\delta_1=0.4$ resulting in severe over-estimation. However, the TIPUP estimators estimate the true zero eigenvalues more accurately, hence are able to take advantage of the more accurate penalty term.
For Model M3 with all weak factors ($\delta_0=\delta_1=0.6$), all IC estimators using the wrong penalty function $g_k(T,d)$ with $\nu=0$ in (ref) identified $(1,1)$ rank pair in all 1000 simulations for all sample sizes. Comparing it with the result shown in Table (ref), the importance of using the right $\nu$ in the penalty function is obvious in this all-weak factor case. Table (ref) also shows that IC estimator using TIPUP outperforms that using TOPUP, when the right penalty function is used. Comparing Tables (ref) and (ref), the iterative IC estimators out-perform the ER estimators when $T=300,500,1000$. It shows the great potential of IC estimators using a proper $\nu$ in the penalty function under the weak factor cases. As mentioned earlier, more investigation is needed to determine the proper value for $\nu$.
To study the empirical property of the optimal constant $c$ discussed in Remark (ref), we simulate data from Model M1 with $r_1=r_2=5$, and set $d_1=d_2=80$, $T=300$ and set the upper bound of the rank as $m^*=10$. In Figure (ref), we show respectively the behavior of $\widehat r_{k,c,j}$ as a function of $(d_{1,j},d_{2,j})$, and of $\widehat r_{k,c,J}$ and of $S_{k,c}$ as functions of $c$, when setting $T_1=T_2=\cdots=T_J=T$ and $d_{1,j}=d_{2,j}$. The rank $r_k$ of the factor process can be determined by considering the mapping $c \to S_{k,c}$ and choosing $\widehat r_{k,c,J}=\widehat r_{k,\widehat c,J}$, where $\widehat c$ belongs to an interval of $c$ implying $S_{k,c} \approx 0$ and therefore the value of $\widehat r_{k,c,J}$ is a constant function of $c$. Similar to hallin2007 and alessi2010improved, we see that the second stability interval always delivers an estimated number $\widehat r_{k,c,J}$ which is closer to the true $r_k = 5$ than the number suggested by the other intervals. That is, the smallest values of $c$ for which $\widehat r_{k,c,j}$ is also close to a constant function of $j$, $j\le J$. Note that the first stability interval always corresponds to the predefined upper bound $m^*$ and it is thus a non-admissible solution.
In this section, we illustrate the proposed procedures using the Fama–French 10 by 10 monthly return series as an example. According to ten levels of market capital (size) and ten levels of book to equity ratio (BE), stocks are grouped into 100 portfolios. The sampling period used in this excises is from January 1964 to December 2015 for a total of 624 months. There are overall 62,400 observed monthly returns used in this analysis. The data is from \\ \centerline{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html.} \\ Similar to wang2019, we subtract from each of the series their corresponding monthly excess market return.
The data was used by wang2019 to demonstrate the estimation of a matrix factor model using a non-iterative TOPUP procedure. They used an estimator similar to our non-iterative TOPUP based ER estimator without a penalty term to estimate the number of factors and found the estimator suggested $(1,1)$ as the rank -- a single factor in the model. In the end they demonstrate the model using $(2,2)$ as the ranks for demonstration of the estimation procedures.
Figure (ref) shows the 1st to 3rd largest eigenvalues of using the non-iterative TIPUP and TOPUP procedures, $h_0^{-1}\hat\tau_{k,m_k}^{*2}$ and $h_0^{-1}\hat\tau_{k,m_k}^{2}$ ($m_k\le 3$), under different lag values $h_0$. The pattern of eigenvalues using the iterative TIPUP and TOPUP are similar, thus is omitted. It can be seen from panel (a) of Figure (ref) that, using TIPUP procedure, the 1st and 2nd largest eigenvalues reach their maximum value at $h_0=1$, and tends to decrease as $h_0$ increases. However, the 3nd largest eigenvalue reaches its maximum at $h_0=2$. In contrast, from panel (b) of Figure (ref), using TOPUP procedure, the 1st to 3nd largest eigenvalues reach the maximum at $h_0=1$. The difference of the patterns of estimated singular values indicates possible severe signal cancellation when using $h_0=1$, according to the suggestions in han2020. Hence, we choose $h_0=2$, and consider IC and ER estimators based on the iterative TIPUP procedure.
Figure (ref) shows the estimated rank $\hat r_k$ of the core factor process with different number of iterations, using IC2(TIPUP) with $\nu=0$ (left figure) and ER1(TIPUP) (right figure). It is seen that the iterative algorithms converge very quickly.
Table (ref) shows the estimated rank pairs using different IC(TIPUP) estimators and different $\nu$ parameter for the penalty function. It is seem that IC1 and IC3 tend to select larger models. These rank estimates do not change when we use $h_0=3$ and $4$.
On the other hand, ER1-ER5 in (ref) produce exactly the same rank estimate $(1,3)$ using $h_0=2$. But these rank estimates change to $(1,1)$ when we use $h_0=3$ and $4$. Figure (ref) shows the estimated eigenvalues $\tau_{k,m_k}^{*2}$, $1\le m_k\le d_k$, using the non-iterative initial TIPUP procedure, for $k=1$ (size factor) and $k=2$ (BE factor). It is seen that for the size factor, the largest eigenvalue is more than 20 folds larger than the second largest eigenvalue. As simulation results in Section 5.2 show, in such an unbalanced case, the ER estimator may find it difficult to find the gap between the true non-zero eigenvalues and the true zero eigenvalues, based on the ratio of the eigenvalues. On the other hand, the IC estimator may fare better in such cases since it is based on the level of estimation error of the true zero eigenvalue.
Overall, it seems that $(2,2)$ or $(3,3)$ are possibly good choices. More detailed analysis, include goodness-of-fit measures, prediction performance and result interpretation, is needed.
In this paper, we develop two rank identification estimators, in an attempt to fill a gap on modelling tensor factor model in the literature. Non-iterative and iterative IC and ER estimators, based on similar ideas of the TOPUP and TIPUP procedures of chen2021factor and han2020 are considered. Theoretical analysis shows that the iterative estimators are much better than the non-iterative estimators. We show that in general the estimators based on TIPUP procedures are better than that based on the TOPUP procedure, due to its fast convergence rate of the estimated eigenvalues under proper conditions. However, in situations when TIPUP procedures also lead to significant signal cancellation, extra care needs to be taken, including increasing the maximum lag $h_0$ in the procedure.
Simulation studies are conducted to compare the finite sample performance of the estimators using the non-iterative and iterative estimation procedures. The results show that the ER estimators based on both TIPUP and TOPUP procedures, and the IC estimators based on TIPUP procedures generally perform very well when all the factors are strong. The ER estimators are better than the IC estimators when some factors are weak, unless one chooses the precise tuning parameter $\nu$ in the IC penalty function, which is a difficult task. When some dominant factors have unrealistically high explanatory power, the ER estimators may not perform well. But the IC estimators still work very well, since the factors are strong. In summary, IC estimator based on iTIPUP shall be used to estimate the number of strong factors, while ER estimators based on iTIPUP are likely to capture weak factors.