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.
90,233 characters · 18 sections · 102 citation commands
Sparse Asymptotic PCA: Identifying Sparse Latent Factors Across Time Horizon in High-Dimensional Time Series
{
}
In the realm of big-data analysis, the study of large-dimensional panel data has gained significant prominence in recent years. Large-dimensional panel data refers to datasets where observations are recorded for multiple individuals or entities over time, resulting in a wealth of information over the spatial and time horizons. Analyzing such datasets is essential for understanding complex economic, financial, and social phenomena. However, as the dimensionality of the data increases, traditional statistical techniques encounter numerous challenges, including multicollinearity, computational complexity, and difficulties in extracting meaningful insights. To address these challenges, latent factor modeling, also referred to as statistical factor modeling in campbell1997econometrics, has emerged as a powerful and versatile approach in the field of panel data analysis. Latent factor models offer a structured framework for dimension reduction and capturing common sources of variation in high-dimensional data. For example, asset returns in finance are often modeled as functions of a small number of factors; see Ross1976 and connor1986performance,connor1988risk. Macroeconomic variables of multiple countries are found to have common components; see stock1989new, gregory1999common, and forni2000reference. In demand systems, the Engle curves can be expressed in terms of a finite number of factors; see lewbel1991rank. As the dimensions of the panel data systems increase, various factor models have been developed to reduce the dimensionalities under different scenarios. See, for example, the approximate factor models in Chamberlain1983, Stock2002a,Stock2002b, bai2002determining, bai2003inferential, fan2013large, lettaupelger2018, pelger2019, and gaotsay2022,gao2023supervised, and the dynamic factor models in forni2000generalized, among others.
One of the most profound challenges in factor modeling is to handle high-dimensional time series data efficiently and effectively. Sparse factor modeling represents a groundbreaking approach to this challenge, offering a powerful technique to extract meaningful insights from complex datasets characterized by a multitude of variables. Sparse factor modeling seeks to identify and capture the underlying structure of data while simultaneously promoting simplicity by selecting only a subset of relevant variables (or features) from the original set. The central idea behind sparse factor modeling is to uncover latent factors that drive the observed data's variation while enforcing a degree of sparsity, meaning that only a few variables are deemed essential to explain the data's structure. This approach is particularly valuable in scenarios where the number of variables far exceeds the number of observations, as it not only reduces computational demands but also enhances the interpretability of the results. Examples of articles concerning this approach include the sparse PCA in jolliffe2003modified, zou2006, shen2008sparse, johnstone2009consistency, witten2009penalized, and ma2013sparse, and the sparse factor analysis with sparse loadings in kristensen2017diffusion and uematsu2022estimation, among others, where all of the aforementioned works are in line with factor analysis with sparse loadings, implying that each factor process is a linear combination of a small subset of the original panel series only.
It is widely known that the key idea in a large-dimensional factor model is that the dimensions of both the timeline and the cross-section of the data are large, and most of the co-movements can be explained by a few factors. These factors and loadings are usually estimated by the conventional PCA method due to its ability to parsimoniously capture much of the information in a large number of variables. However, the resulting latent PCs or factors are typically linear combinations of all cross-sectional units/variables, which are usually hard to interpret. The sparse PCA restricts the cardinality of the weight vectors for the PCs, or equivalently, the loadings for the factor processes, so the PCs are sparse linear combinations of the underlying variables. Consequently, the PCs or factors are only linear combinations of a small subset of the cross-sectional units, which facilitates the interpretations of the resulting PCs or factors. Nevertheless, there is little literature concerning the justificaton of the sparsity assumption in the loadings. For example, the assumption may not be appropriate for financial returns, where the exposure of the returns to a market factor is universal and non-sparse, as pointed out in pelger2022interpretable. {To illustrate, we examine the risk factors of the daily return data of the S&P 500 stocks studied in pelger2019, with \( N = 332 \) stocks over \( T = 3273 \) time points. We first apply PCA to the panel, and the estimated factors \( f_t \) are shown in Figure S.1 of the Supplement. We then regress each return series \( x_{i,t} \) on \( f_t \), and the resulting \( p \)-values for testing the significance of the regression coefficients are all close to zero, implying that the dependence of the returns on the factors may be non-sparse. Next, we define a sparse factor \( f_t^s \), where \( f_t^s = f_t \) if \( f_t \) is among the 500 largest values in absolute magnitude, and \( f_t^s = 0 \) otherwise. We find that ${\sum_{t=1}^{3273}(f_t^s)^2}/{\sum_{t=1}^{3273} f_t^2} = 79.4\%$, indicating that the top 500 factor points account for 79.4% of the total factor variance over the entire time horizon. Finally, we calculate the \( R^2 \) values for regressing \( x_{i,t} \) on \( f_t \) and \( f_t^s \), respectively. The results are presented in Figure (ref). From the figure, we observe that the \( R^2 \) values based on only 500 factor points represent a substantial proportion of those based on the full factor series. The average ratio of \( R^2 \) values from the sparse to the full factor regressions is 77.5%, suggesting that the explanatory power of the sparse factor over time is substantial and deserves further attention. }
In view of the above discussion, this paper marks a further development in the sparse factor modeling of large-dimensional panel data from a different perspective. Motivated by the empirical success of a general approximate factor model, see pelger2019 and the references therein, we also adopt such a similar framework. Unlike the sparse factor modeling approach with sparse loadings in existing literature, we assume the latent factor processes are sparse over time, while the loadings can be nonzero in general settings. {Specifically, over a given time period, cross-sectional economic or financial series may exhibit varying degrees of co-movement. During certain periods, weak co-movement allows idiosyncratic components to dominate, resulting in relatively lower systematic risk—see, for example, the empirical evidence in campbell2001have and goyal2003idiosyncratic. Since the loading vector for each series reflects a linear combination of the underlying factor series, sparsity—where only the most significant factors are retained—implies that loadings are primarily driven by a subset of factors. This formulation is intuitive in financial applications, where asset returns often alternate between periods of high volatility and relative calm. Co-movements tend to intensify during major events such as interest rate hikes by the Federal Reserve or financial crises related to debt, currency, or mortgage markets. Our proposed method offers a useful framework for uncovering connections between major economic or policy events and the evolving co-movement structure of economic or financial systems.}
In this paper, we present a sparse APCA technique by formulating a sparse factor modeling framework with sparse factors over time horizons. {Under the $\ell_0$-constraint imposed on the factors, we propose a truncated power method to estimate the sparse factors in the one-factor case, corresponding to the hard-thresholding technique widely used in the statistics and machine learning literature. As noted in johnstone2009consistency, a key advantage of hard-thresholding over soft-thresholding is that it preserves the magnitude of the retained signals, whereas soft-thresholding may introduce bias into the estimators.} While hard-thresholding has been employed in studies such as johnstone2009consistency and ma2013sparse, their focus is primarily on independent and identically distributed (i.i.d.) data with independent/white noises, where the distinction between factors and loadings therein is minimal. {From a theoretical perspective, the presence of serial correlation in both the data and the factors poses additional challenges in establishing the convergence rate of the sample covariance matrix and conducting inference on the estimated loadings. In particular, large deviation inequalities such as Bernstein’s inequality cannot be directly applied to time series data. Furthermore, the thresholding values proposed in johnstone2009consistency and ma2013sparse, which rely on the assumption of normality, are not applicable in our setting, as time series data are generally not normally distributed.} For multi-factor cases, we develop a sequential deflation estimation procedure that integrates the truncated power method with projection techniques. Additionally, we introduce a data-driven approach to determine the sparsity of risk factors over time using a novel cross-sectional cross-validation method. Theoretical properties of the proposed estimators are established under mild conditions for dependent data as both the dimension and sample size tend to infinity. Simulated examples are used to illustrate the proposed method. We empirically examine specific events associated with the identified sparse factors that systematically influence the stock market. Our findings reveal nine key time factors that may systematically affect the financial market.
{The proposed framework is fundamentally different from existing sparse PCA and sparse factor models in both formulation and interpretation, and it also employs a distinct estimation technique. In the sparse PCA literature (e.g., zou2006, shen2008sparse, witten2009penalized), the true sparse eigenvectors are treated as constant vectors, and regularized algorithms are used for their estimation. In sparse factor models (e.g., kristensen2017diffusion, uematsu2022estimation), the factor loadings are assumed to be sparse and non-random, with estimation based on $\ell_1$-regularization. The sparse proximate factor model of pelger2022interpretable assumes sparsity only in factor weights, while the resulting factor and loading processes remain dense. Time-varying factor models, such as su2017time, allow for temporal dynamics of the loadings but do not impose sparsity. In contrast, our sparse APCA framework assumes that the factor processes themselves exhibit sparsity—i.e., each factor entry is either an active random variable or zero at a given time. A zero factor value suggests weaker co-movement at that time point. Our method imposes $\ell_0$-constraints directly on the factors and is thus theoretically distinct from prior approaches, which focus on cross-sectional sparsity. This temporal sparsity structure enables us to identify specific time points where systematic risk factors significantly influence the panel.}
The contributions of this paper are multi-fold. First, to the best of our knowledge, this is the first study to consider sparse factors over the time horizon, providing economists with a novel approach to bridging the gap between time and risk in economic and financial systems. Second, unlike the $\ell_1$-penalized methods in most sparse PCA or factor modeling literature, we propose a truncated power method based on an $\ell_0$-constraint for the one-factor case, along with a sequential optimization procedure for estimating multiple factors. Our approach differs from the hard-thresholding methods for i.i.d. data in johnstone2009consistency and ma2013sparse by explicitly accounting for time series data with both cross-sectional and dynamic dependence. Consequently, our method extends the asymptotic PCA technique developed by connor1986performance,connor1988risk to high-dimensional factor analysis. Third, instead of relying on traditional cross-validation to determine sparsity structures, we introduce a novel cross-validation approach that partitions the spatial dimensions. Fourth, most thresholding techniques require specifying a threshold value, as seen in johnstone2009consistency and ma2013sparse. In contrast, our method introduces a selection criterion that determines the number of nonzero factor points without requiring pre-specified threshold values. Theoretically, we establish the consistency of the proposed estimators as the dimension and sample size approach infinity. Importantly, we provide statistical guarantees for the proposed cross-sectional cross-validation method, demonstrating its ability to consistently estimate the sparsity structure over the time horizon of factor models.
The rest of the paper is organized as follows. Section (ref) introduces the sparse factor modeling framework and its estimation procedure. Section (ref) presents asymptotic properties of the estimators obtained in Section (ref). Section (ref) studies the finite-sample performance of the proposed approach via simulation, and (ref) illustrates the proposed procedure with an empirical application. Section (ref) concludes. Some illustrative examples of sparse factors over time horizons, additional tables and figures used in the numerical analyses, all theoretical proofs of the main theorems, and a description of the text data extracted from {\it CNN} are relegated to an online Supplement.
{\bf Notation:} We use the following notation. For a $p\times 1$ vector ${\mathbf u}=(u_1,..., u_p)'$, $\|{\mathbf u}\|_1=\sum_{i=1}^p|u_i|$ is the $\ell_1$-norm and $\|{\mathbf u}\|_\infty=\max_{1\leq i\leq p}|u_i|$ is the $\ell_\infty$-norm. ${\mathbf I}_p$ denotes the $p\times p$ identity matrix. For a matrix ${\mathbf H}$, its Frobenius norm is $\|{\mathbf H}\|=[\mathrm{trace}({\mathbf H}'{\mathbf H})]^{1/2}$ and its operator norm is $\|{\mathbf H} \|_2=\sqrt{\lambda_{\max} ({\mathbf H}' {\mathbf H} ) }$, where $\lambda_{\max} (\cdot) $ denotes the largest eigenvalue of a matrix, and $\|{\mathbf H}\|_{\min}$ is the square root of the minimum non-zero eigenvalue of ${\mathbf H}{\mathbf H}'$. $|\mathbf{H}|$ denotes the element-wise absolute value of $\mathbf{H}$. The superscript ${'}$ denotes the transpose of a vector or matrix. We also use the notation $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$ or $a$ and $b$ have the same order of stochastic bound when they are random variables.
Let $x_{i,t}$ be the $i$-th unit of the cross-sectional panel of time series ${\mathbf x}_t=(x_{1,t},...,x_{N,t})'\in R^N$ at time $t$, for example, $x_{i,t}$ can be the stock return of the $i$-th asset at time $t$, we consider the following approximate factor model:
where $x_{i,t}$ is the only observed datum for the $i$-th cross-section at time $t$ ($i=1,...,N$; $t=1,...,T$), ${\mathbf f}_t$ is an $r$-dimensional vector of common or systematic risk factors, $\boldsymbol{\lambda}_i$ is an $r$-dimensional vector of factor loadings, and $e_{i,t}$ is the idiosyncratic component of $x_{i,t}$ that cannot be explained by the common risk factors. We may combine all the cross-sectional units together and write the above equation as
where $\boldsymbol{\Lambda}=(\boldsymbol{\lambda}_1,...,\boldsymbol{\lambda}_N)'$ and ${\mathbf e}_t=(e_{1,t},...,e_{N,t})'$. Assume a panel data set of $T$ time-series observations and $N$ cross-sectional observations, denoted as ${\mathbf X} \in R^{T\times N}$, has a factor structure with $r$ common factors. Let ${\mathbf F}=({\mathbf f}_1,...,{\mathbf f}_T)'$ and ${\mathbf e}=({\mathbf e}_1,...,{\mathbf e}_T)'$, then Models ((ref)) and ((ref)) can be written in the following matrix form
For ease of notation, we also denote ${\mathbf X}=(\underline{{\mathbf x}}_1,...,\underline{{\mathbf x}}_N)$, ${\mathbf F}=(\underline{{\mathbf f}}_1,...,\underline{{\mathbf f}}_r)$, and $\boldsymbol{\Lambda}=(\underline{\boldsymbol{\lambda}}_1,...,\underline{\boldsymbol{\lambda}}_r)$ when referring to their columns. In the econometrics/statistical and finance literature, Model ((ref)) is usually estimated using the Principal Components or Asymptotic Principal Components estimation method (see connor1986performance,connor1988risk, bai2002determining, bai2003inferential, and fan2013large, among others). As a result, all the elements in the estimators $\widehat{\mathbf F}$ and $\widehat\boldsymbol{\Lambda}$ are usually nonzero, making it difficult to interpret the dynamic and cross-sectional relationships of the data. Moreover, the role of the idiosyncratic terms is rarely considered or even ignored in characterizing the individual dynamics of the series.
{ In Model ((ref)), only the panel ${\mathbf X}$ is observed, while the factor ${\mathbf F}$ and loading $\boldsymbol{\Lambda}$ are latent and unidentifiable up to a rotation: for any invertible $r \times r$ matrix ${\mathbf H}$, the pair $({\mathbf F}{\mathbf H}', \boldsymbol{\Lambda}{\mathbf H}^{-1})$ yields the same model. To ensure identification, two commonly used sets of conditions are adopted in the PCA framework:
See bai2002determining, fan2013large, and jiang2023revisiting, among others. These identification restrictions help resolve the rotational indeterminacy when the leading eigenvalues of the covariance matrix are distinct; see bai2013principal. Under either set of conditions, PCA provides consistent estimates of the common component $\widehat{\mathbf F} \widehat\boldsymbol{\Lambda}'$.
In this paper, we adopt the first set of conditions in ((ref)) and apply the Asymptotic Principal Component Analysis (APCA) approach to uniquely identify ${\mathbf F}$ and $\boldsymbol{\Lambda}$. Suppose the true model is ${\mathbf X} = {\mathbf F}^*\boldsymbol{\Lambda}^{*'} + {\mathbf e}$, jiang2023revisiting shows that there exists a rotation matrix ${\mathbf H}$ such that the estimated APCA factors $\widehat{\mathbf F}$ satisfy $\widehat{\mathbf F}'\widehat{\mathbf F}/T = {\mathbf I}_r$ and are consistent for ${\mathbf F}^*{\mathbf H}'$. Redefining ${\mathbf F} = {\mathbf F}^*{\mathbf H}'$ and $\boldsymbol{\Lambda} = \boldsymbol{\Lambda}^*{\mathbf H}^{-1}$ transforms the model to ${\mathbf X} = {\mathbf F}\boldsymbol{\Lambda}' + {\mathbf e}$, which allows us to impose sparsity assumptions directly on ${\mathbf F}$ without encountering identification issues.
Our first goal is to extend the APCA framework of connor1986performance, connor1988risk to high-dimensional settings using modern machine learning techniques. Moreover, since factors and loadings in approximate factor models are typically difficult to interpret, our second goal, building on the framework of bai2013principal, is to impose additional sparse structure on the factors. This enables interpretation by linking the estimated factors to specific time events in the data.
}
In this paper, unlike the existing approach in sparse factor modeling where the loadings are sparse, we assume the systematic risk factor process ${\mathbf f}_t$ is sparse along the time horizon as it characterizes the co-movement of the panel dynamically. When all components of ${\mathbf f}_t$ are zero at some certain time point, it implies that the panel is driven dominantly by individual factors or the corresponding idiosyncratic terms rather than a systematic one. On the other hand, for the economic and financial panel series, the nonzero factor ${\mathbf f}_t$ at certain data points can be treated as a systematic response of the panel to some important and influential economic or financial events or outcomes.
We first introduce some further notation used in the estimation procedure. For an $N\times 1$ vector ${\mathbf v}=(v_1,...,v_N)'$, $\|{\mathbf v}\|_2=\sqrt{\sum_{i=1}^2v_i^2}$ is the Euclidean norm and $\|{\mathbf v}\|_0=\text{card}\{\text{support}({\mathbf v})\}$ is the cardinality or number of non-zero elements in ${\mathbf v}$. Let $\mathbb{V}$ be the set of $N\times r$ semi-orthogonal matrices and $\mathbb{V}^{\perp}$ be the set consisting of the orthogonal complements of the ones in $\mathbb{V}$. {For a matrix ${\mathbf V} = ({\mathbf v}_1, \ldots, {\mathbf v}_r) \in \mathbb{V}$ and a vector ${\mathbf s} = (s_1, \ldots, s_r)'$, $\|{\mathbf V}\|_0 \leq {\mathbf s}$ means that each column ${\mathbf v}_i$ satisfies $\|{\mathbf v}_i\|_0 \leq s_i$ for $1 \leq i \leq r$.}
In the following subsection, we first formulate the estimation procedure in the one-factor case, i.e. the number of factors $r=1$ in Model ((ref)), and the multi-factor case can be carried out based on the one-factor estimation procedure.
We consider the case when $r=1$ and the first factor process over the timeline is $\underline{{\mathbf f}}_1=(f_{1,1},...,f_{1,T})'$. Let ${\mathbf S}={{\mathbf X}{\mathbf X}'}/{(NT)}$, it follows that $\underline{{\mathbf f}}_1/\sqrt{T}$ can be approximated by the first normalized eigenvector of the $T\times T$ positive semi-definite matrix ${\mathbf S}$ according to Model ((ref)). We assume $\|{\mathbf v}\|_0\leq s_1$, it suffices to solve the following optimization problem:
and consequently, the estimator for $\underline{{\mathbf f}}_1$ is denoted as $\underline{\widehat{\mathbf f}}_1=\sqrt{T}\widehat{\mathbf v}$ such that $\underline{\widehat{\mathbf f}}_1'\underline{\widehat{\mathbf f}}_1/T=1$, which satisfies the first condition in ((ref)). Note that the solutions in ((ref)) are the same as those obtained by PCA in bai2002determining and gaotsay2022, among others, if the constraint $\|{\mathbf v}\|_0\leq s_1$ is removed. Furthermore, the problem in ((ref)) is also equivalent to finding the largest eigenvalue associated with an $s_1$-sparse eigenvector of the matrix ${\mathbf S}$:
which is a non-convex problem in general. In fact, it is NP-hard because it can be reduced to the subset selection problem for ordinary least-squares problem; see moghaddam2006generalized.
The optimization in ((ref)) is still a sparse PCA problem symbolically, and numerous methods have been developed to obtain sparse eigenvectors during the past decades. See the regularization method with elastic net in zou2006 and the penalized matrix decomposition (PMD) algorithm using $\ell_1$-penalty in witten2009penalized, among others.
In this paper, we extend the APCA method of connor1986performance,connor1988risk to high dimensions by adopting an $\ell_0$-constraint on the sparse factors. We propose to use the truncated power method introduced in yuan2013 to solve the first normalized sparse eigenvector in problem ((ref)). The method is similar to the classical power method but includes an additional truncation operation to ensure sparsity. The rationale for this is as follows. First, the conventional APCA is conducted based on the following matrix perturbation formulation:
where ${\mathbf S}$ is the matrix obtained by the noisy observation ${\mathbf x}_t$, $\boldsymbol{\Sigma}$ is a symmetric matrix whose eigenvectors are the true ones, and ${\mathbf E}$ is a random perturbation. The decomposition of ((ref)) is formulated in (S.1) of the online Supplement. If the true largest eigenvector ${\mathbf v}_1$ of $\boldsymbol{\Sigma}$ is sparse, then it is natural to recover ${\mathbf v}_1$ from the noisy matrix ${\mathbf S}$. This recovery is guaranteed when the error ${\mathbf E}$ is of a smaller order using the well-known $\sin \theta$ theorem in davis1970rotation under the approximate factor model. Second, for any given vector ${\mathbf u}_0\in R^{N}$, the power method estimates the first eigenvector by
and a normalized ${\mathbf u}_{k}$ converges to the first eigenvector of ${\mathbf S}$. To see this, note that there exist $c_1$,..., $c_N$ such that \[{\mathbf u}_0=c_1\widehat{\mathbf v}_1+c_2\widehat{\mathbf v}_2+...+c_N\widehat{\mathbf v}_N,\] where $\widehat{\mathbf v}_1,....,\widehat{\mathbf v}_N$ are the $N$ eigenvectors associated with the eigenvalues $\{\widehat\lambda_i,i=1,...,N\}$ of ${\mathbf S}$. Then, \[{\mathbf u}_k={\mathbf S}^k{\mathbf u}_0=\sum_{i=1}^Nc_i{\mathbf S}^k\widehat{\mathbf v}_i=\widehat\lambda_1^k\left[c_1\widehat{\mathbf v}_1+\sum_{i=2}^N c_i(\frac{\widehat\lambda_i}{\widehat\lambda_1})^k\widehat{\mathbf v}_i\right].\] Assuming that $|\widehat\lambda_1|>|\widehat\lambda_i|$ for $i\geq 2$, it follows that the direction of ${\mathbf u}_k$ converges to that of $\widehat{\mathbf v}_1$.
Motivated by the above discussion, we modify the power method by adding a truncation in each iteration, and the procedure is given in Algorithm (ref) below.
The procedure, presented in Algorithm (ref), generates a sequence of intermediate $s_1$-sparse eigenvectors $\{{\mathbf u}_1, {\mathbf u}_2,...\}$ from an initial sparse approximation ${\mathbf u}_0$. At each time stamp $t$, the intermediate vector ${\mathbf u}_{t-1}$ is multiplied by ${\mathbf S}$, and then the entries are truncated to zeros except for the largest $s_1$ entries. The resulting vector is then normalized to unit length. Finally, the estimated first factor process is $\underline{\widehat{\mathbf f}}_1=\sqrt{T}\widehat{\mathbf v}_1$, where $\widehat{\mathbf v}_1$ is the output one in Algorithm (ref). Instead of pre-specifying threshold values as in johnstone2009consistency and ma2013sparse, we rely on the sparsity parameter $s_1$, which is unknown in practice. To address this, we will propose a cross-validation procedure to determine the optimal value of $s_1$ later.
Algorithm (ref) can consistently estimate sparse eigenvectors under a factor structure when the corresponding eigenvalues diverge. However, this property may not extend to general covariance matrices. For instance, it is possible to construct a diagonal covariance matrix with sparse eigenvectors, where the algorithm may fail to correctly identify the leading sparse eigenvectors for some initial vectors. For instance, consider ${\mathbf S}=\mbox{diag}(2,1)$. In this case, the leading eigenvector is ${\mathbf v}_1=(1,0)'$, but the algorithm will incorrectly output ${\mathbf v}_1=(0,1)'$ if the initial vector is ${\mathbf u}_0=(1,3)'$. This occurs because the top eigenvalue is not sufficiently large. However, if ${\mathbf S}=\mbox{diag}(4,1)$, where the top eigenvalue is significantly larger—as is common in a factor model, the algorithm correctly identifies the leading eigenvector.
In practice, the convergence criterion in Algorithm (ref) needs to be given first. A common criterion for the convergence of Algorithm (ref) is that the two eigenvector iterates ${\mathbf u}_{t}$ and ${\mathbf u}_{t-1}$ for $t\geq 1$ satisfy
for a small ${{\varepsilon}}>0$. Simulation results suggest that the convergence is not sensitive to a sufficiently small ${{\varepsilon}}$, and the numerical results in Section (ref) indicate that the algorithm works well when we take ${{\varepsilon}}=10^{-3}$ in ((ref)).
{In this section, we consider the general case where multiple factor processes are present in Model ((ref)), i.e., when $r > 1$. We extend the truncated power method to accommodate the multi-factor setting.} In the presence of more than one common factor process, we need to apply Algorithm (ref) multiple times and extract all the common factors in a sequential way. Specifically, when $\widehat{\mathbf v}_1$ is given via Algorithm (ref), we subtract the projection on first factor component of the panel and the residual would be $\widetilde{\mathbf X}={\mathbf X}-\widehat{\mathbf v}_1\widehat{\mathbf v}_1'{\mathbf X}=({\mathbf I}_T-\widehat{\mathbf v}_1\widehat{\mathbf v}_1'){\mathbf X}$. The resulting scaled covariance $\widetilde{\mathbf S}_1=({\mathbf I}_T-\widehat{\mathbf v}_1\widehat{\mathbf v}_1'){\mathbf S}({\mathbf I}_T-\widehat{\mathbf v}_1\widehat{\mathbf v}_1')$. Then the second eigenvector $\widehat{\mathbf v}_2$ can be obtained by solving the following optimization problem ((ref)) below:
and the second eigenvector is $\widehat{\mathbf v}_2=\widetilde{\mathbf v}_2/\|\widetilde{\mathbf v}_2\|_2$. Consequently, the second estimated factor process is $\underline{\widehat{\mathbf f}}_2=\sqrt{T}\widehat{\mathbf v}_2$. We use the normalization condition ${\mathbf v}'({\mathbf I}_T-\widehat{\mathbf v}_1\widehat{\mathbf v}_1'){\mathbf v}=1$ to maximize the additional variance of the original matrix ${\mathbf S}$, which is the same as the deflation framework in mackey2008. In other words, we need to eliminate the effect of the previous eigenvector directions by a projection method. Consequently, we formulate the procedure in Algorithm (ref) below.
From the optimization problem in ((ref)) and the procedure in Algorithm (ref), we cannot enforce the orthogonality and sparsity at the same time, which is similar to the case in the sparse PCA framework. When there is no $\ell_0$-constraint in ((ref)) and ((ref)), we can easily obtain that $\widehat{\mathbf v}_1$ and $\widehat{\mathbf v}_2$ are the two normalized and orthogonal eigenvectors associated with the top two eigenvalues of ${\mathbf S}$, and the approach reduces to the traditional APCA method. {One of the differences between the one-factor case in ((ref)) and the multi-factor case in ((ref)) is that the estimated factors and loadings in the one-factor case satisfy the conditions in ((ref)), whereas those obtained in the multi-factor case may not. It is important to note that this discrepancy is only a finite sample outcome. Asymptotically, we can demonstrate that the multi-factor case still satisfies the identification conditions in ((ref)). Similar outcomes are also observed in the sparse PCA literature, including zou2006 and witten2009penalized, as well as in the sparse factor modeling frameworks discussed in kristensen2017diffusion.}
For each sparse vector $\widehat{\mathbf v}_i$ obtained in Algorithm (ref), the estimated factor process is obtained as $\underline{\widehat{\mathbf f}}_i=\sqrt{T}\widehat{\mathbf v}_i$. It is not hard to see that the number of iterations is $r$, and the $i$th iteration outputs a sparse estimator (multiplied by $\sqrt{T}$) of the $i$th column of ${\mathbf F}$. For the optimization problem in Step 4 of Algorithm (ref), by the argument in Lemma 1 of the Supplement, the matrix ${\mathbf B}_i$ in each step is a symmetric one. Therefore, we modify the truncated power method in Algorithm (ref) and propose the following algorithm to obtain a sparse eigenvector for Step 4 of Algorithm (ref).
Note that the symmetric matrix ${\mathbf B}_i$ in Algorithm (ref) is not strictly positive definite, the half-inverse ${\mathbf B}_i^{-1/2}$ should be taken as a generalized one in line with the Moore-Penrose generalized inverse matrix.
The estimation method in Section (ref) depends on a known number of factors $r$. In practice, $r$ is unknown and we need to develop a data-driven method to estimate it. There are several useful methods developed in the past decades to estimate the number of factors including the information criterion in bai2002determining, the random matrix theory method in onatski2010determining, the ratio-based method in lam2012factor and Ahn2013, and the white noise testing approach in gao2022modeling, among others.
We introduce two commonly used methods. The first is an information criterion (IC) approach adapted from bai2002determining, with an added $\log(T)$ penalty to account for the scale of $s$. The number of factors $r$ is estimated by
where $K$ is a pre-specified upper bound, and $\widehat{\mathbf F}_k$, $\widehat\boldsymbol{\Lambda}_k$ are the estimated factors and loadings for $k$ factors. The other method is based on the eigenvalue ratios introduced by lam2012factor and Ahn2013. Specifically, let $\widehat\lambda_1\geq\widehat\lambda_2\geq...\geq\widehat\lambda_T$ be the $T$ sample eigenvalues of ${\mathbf X}{\mathbf X}'$, we adopt the ratio-based method to estimate $r$ by
where $K$ is a prescribed upper bound as in ((ref)) to control the stability of the ratios. In practice, we may choose $K=\min(T,N)/3$ under the assumption that the number of factors $r$ is usually not large.
For a given number of factors $r$, which can be estimated by the methods in Section (ref) above, if the cardinality of each column in ${\mathbf F}$ is known, we can obtain the estimated factors $\widehat {\mathbf F}$ by solving the optimization problem using the proposed algorithms in Section (ref). The resulting estimated loading matrix is $\widehat\boldsymbol{\Lambda}'=(\widehat{\mathbf F}'\widehat{\mathbf F})^{-1}\widehat{\mathbf F}'{\mathbf X}$, which can be obtained by the Ordinary Least-Squares (OLS) method. When the true factors are sparse and orthogonal, the theoretical results in Section (ref) below suggest that $\widehat\boldsymbol{\Lambda}\approx {\mathbf X}'\widehat{\mathbf F}$, implying that each column loading vector is a sparse linear combination of the panel data over the timeline.
{For theoretical purposes, we may assume \( s_i / T = \delta_i \in (0,1) \), analogous to the dimensional asymptotics commonly considered in random matrix theory.} In practice, the cardinality of each column in ${\mathbf F}$ is unknown, and we may choose them by the cross-validation method which is commonly used in machine learning literature. In fact, even though the sparsity parameters $s_1, ..., s_r$ of $r$ factor processes may not be the same, it is still a convenient way to assume {that $s_1=...=s_r=s_0\asymp T$} which can simplify the cross-validation procedure significantly. According to the theoretical results in Section (ref) and the proofs in the Supplement, when the $r$ sparsity parameters are distinct, the convergence of the estimator for $s^*$ can be achieved, where $s^*=\max(s_1,...,s_r)$. Once $s^*$ is estimated in this scenario, it can be fixed, allowing the estimation of the second-largest sparsity parameter using the proposed method. For simplicity, in this section, we assume that the sparsity parameters for each column are identical. The cross-validation procedure is outlined as follows.
For each fixed $s\in \mathbb{S}$ where $\mathbb{S}$ is a chosen candidate set for the sparsity parameter $s$, we first randomly divide the panel into two segments consisting of one training sample with $N_1$ components and the other a testing one with $N_2$ components, where $N_1\asymp N_2\asymp N$ and $N_1+N_2=N$. Equivalently, we partition the data matrix ${\mathbf X}\in R^{T\times N}$ into $\widetilde{\mathbf X}_1$ and $\widetilde{\mathbf X}_2$, where $\widetilde{\mathbf X}_1\in R^{T\times N_1}$ and $\widetilde{\mathbf X}_2\in R^{T\times N_2}$. Let $\widetilde{\mathbf F}_1^{s}$ be the estimated factors based on the training sample $\widetilde{\mathbf X}_1$ and $\widetilde\boldsymbol{\Lambda}_1$ be the estimated loading matrix. Define the testing error as
When we range the parameter $s$ over the candidate set $\mathbb{S}$, the optimal number of nonzero factors, denoted by $\widehat s$, is the one that produces the smallest errors in ((ref)). On the other hand, we do not expect $s$ to be a large one to avoid the overfitting issue. Therefore, we define $g(N_1,T)$ as a penalty function, and let \[PC(s)=R(s,\widetilde{\mathbf F}_1^{s})+rC_T\frac{s}{T}g(N_1,T),\] for some $C_T>0$ which will be specified later. Then, $\widehat s$ is defined as
implying that we choose the optimal $s$ which results in minimal testing errors in the cross-validations. There are many choices for the penalty function $g(N_1,T)$ as shown in Theorem (ref) of Section (ref) below. We will discuss the required properties of $g(N_1,T)$ such that we can establish the consistency of $\widehat s$ in theory.
{In practice, we may restrict the candidate set $\mathbb{S}$ to the range $c_1T \leq s \leq c_2T$, where $c_1, c_2 \in (0,1)$ are prescribed constants. We may choose $C_T=\log(T)$ in order to balance the magnitude of the scale of $s$.} Note that the cross-validation method is different from the classical one in the machine learning literature where the samples are partitioned along the time horizon, whereas we partition the samples over the cross-section in this paper. In fact, we can even perform more cross-validation experiments in order to obtain an optimal parameter $s$. For example, we can choose a sufficiently large integer $J>0$ and obtain $J$ testing errors as that in ((ref)) for each fixed $s$ by partitioning the samples for $J$ times. We may denote the testing error by $R_j(s,\widetilde{\mathbf F}_{1,j}^{s})$ for the $j$-th random partition and the average testing error for $s$ is defined as
where $\widetilde{\mathbf F}_{1,j}^{s}$ is the estimated factors based on the training sample in the $j$-th cross-validation. Then, the optimal number of nonzero factors $\widehat s$ is the one that produces the smallest error of the information criterion in ((ref)) by replacing the $R(s,\widetilde{\mathbf F}_1^{s})$ therein by $R^J(s)$ in Equation ((ref)) for $s\in \mathbb{S}$.
In this section, we present the consistency and asymptotic bounds of the estimators from Section (ref) as $N,T\rightarrow\infty$. We first derive the asymptotic results assuming known $r$ and $s_i$, and then establish the consistency of their estimators. Mathematical proofs are provided in the Supplement.
To derive the asymptotic properties of the proposed method and estimators from Section (ref), we need the following high-level assumptions. These assumptions can be readily verified by some standard and more primitive assumptions, which will be discussed in detail later. Most of the assumptions below are commonly used in the PCA or approximate-factor modeling literature, and some of which are stronger than those in bai2002determining in order to show that the estimated sparse factor $\widehat{\mathbf F}$ is close to ${\mathbf F}$. We will use $C$ or $c$ to denote a generic positive constant the value of which may change at different places.
Assumption (ref) is standard to characterize the dynamic dependence of the factor processes. See, for example, the Lemma 1 in the appendix of gao2019banded. Assumptions (ref)-(ref) establish an identification condition and a sparsity structure for the true factor process ${\mathbf F}$. {While the factors are sparse, Assumption (ref) is not restrictive, as it is reasonable to assume that \( s_i \) diverges. Moreover, the factors can be normalized such that Assumption (ref) holds under the condition \( s_i / T = \delta_i \in (0,1) \), as specified in Assumption (ref), which resembles the setting commonly adopted in random matrix theory such as those in bai2010spectral. Alternatively, one could adopt the framework of uematsu2022estimation, which models factors as weak with varying strengths. However, specifying the strength levels in practice can be challenging. Instead, we adopt the random matrix theory perspective and assume that the sparsity parameter grows at a fractional rate relative to the time dimension \( T \). This is purely a theoretical device and does not impact our primary goal of identifying nonzero factors over time.} Assumption (ref) specifies that the top eigenvalues of the covariance of the panel data are distinct in order to avoid multiplicity of eigenvalues. Assumptions (ref) and (ref) provide identifiability for the Model ((ref)) so that the factors and the loading matrix can be uniquely determined asymptotically. See bai2013principal for a detailed argument. {Assumption (ref) is primarily imposed to simplify theoretical derivations. It can be relaxed to weaker conditions, such as Conditions (c.i)–(c.iii) in bai2013principal, or replaced with strong mixing assumptions in both spatial and temporal dimensions, allowing the use of the Bernstein-type inequality in merlevede2011bernstein for the analysis. Similar independence assumptions are also found in onatski2012asymptotics and Assumption 3 of huang2022scaled. Assumption (ref) facilitates the derivation of a Bernstein-type concentration bound for establishing convergence rates of the estimated factors, and can likewise be weakened following merlevede2011bernstein.}
The following theorem establishes the consistency of the first estimated factor process.
{To maintain clear interpretation, we continue to use $s_1$ instead of $s_1\asymp T$ in the theorem. From Theorem (ref), we see that the angle between the estimated direction $\widehat{\mathbf v}_1$ and the true one ${\mathbf v}_1$ is asymptotically equal to zero or $\pi$ if $s_1\log(T)/(NT)\asymp \log(T)/N\rightarrow 0$. In fact, we can rewrite $\widehat{\mathbf v}_1'{\mathbf v}_1$ as $\cos(\theta)$ where $\theta$ is the angle between $\widehat{\mathbf v}_1$ and ${\mathbf v}_1$. Some remarks for Theorem (ref) are as follows.}
Next, we present the theorem concerning the consistency of all the estimated factors along the time horizon. For the distance between two matrices ${\mathbf H}_1$ and ${\mathbf H}_2$, there are several measures that are used in the literature. For example, we adopt the discrepancy measure used by pan2008modelling: for two $T\times r$ semi-orthogonal matrices ${\bf H}_1$ and ${\bf H}_2$ satisfying the condition ${\bf H}_1'{\bf H}_1={\bf H}_2'{\bf H}_2={\mathbf I}_{r}$, the difference between the two linear spaces $\mathcal{M}({\bf H}_1)$ and $\mathcal{M}({\bf H}_2)$ is measured by
Note that $D(\mathcal{M}({\bf H}_1),\mathcal{M}({\bf H}_2)) \in [0,1].$ It is equal to $0$ if and only if $\mathcal{M}({\bf H}_1)=\mathcal{M}({\bf H}_2)$, and to $1$ if and only if $\mathcal{M}({\bf H}_1)\perp \mathcal{M}({\bf H}_2)$. By Lemma A1(i) in pan2008modelling, $D(\cdot,\cdot)$ is a well-defined distance measure on some quotient space of matrices. Alternatively, we may also adopt the measure
which is the Frobenius norm of the difference between the projection matrices of two spaces and is also a well-defined distance between linear subspaces. In addition, if we denote the singular values of ${\mathbf H}_1'{\mathbf H}_2$ by $\{\sigma_i\}_{i=1}^r$, in descending order, then the principal angles between $\mathcal{M}({\mathbf H}_1)$ and $\mathcal{M}({\mathbf H}_2)$, $\boldsymbol{\Theta}({\mathbf H}_1,{\mathbf H}_2)=\mbox{diag}(\theta_1,...,\theta_r)$, are defined as $\mbox{diag}\{\cos^{-1}(\sigma_1),...,\cos^{-1}(\sigma_r)\}$; see, for example, Theorem I.5.5 of stewart1990matrix. The squared Frobenius norm of the so-called $\sin\boldsymbol{\Theta}$ distance, defined as
can also be used to measure the distance between two linear spaces. In fact, if $r$ is finite, the distances in ((ref))--((ref)) are equivalent, because
Therefore, we shall only use the distance in ((ref)) to present our theoretical results in the main article.
Next, we establish the theoretical results for the estimated loading matrix $\widehat\boldsymbol{\Lambda}$ in the following theorem.
The following corollary establishes the consistency of the estimated number of factors obtained via ((ref)) or ((ref)).
Finally, we provide the consistency of the cross-validation method in estimating the sparsity parameters $s_1,...,s_r$. We consider the case where the sparsity parameters of each column of ${\mathbf F}$ are the same, i.e., $s_1=...=s_r=s_0$, as discussed in Section (ref). Otherwise, the following theorem states the consistency for estimating the largest sparsity parameter among the $r$ columns. {Since $s_0 \asymp T$, it is generally challenging to derive consistency of $\widehat{s}$ when both $N$ and $T$ grow. Therefore, we restrict attention to the asymptotic regime where $N \to \infty$ and $T$ is large but fixed.}
In this section, we illustrate the finite-sample properties of the proposed methodology under different choices of $N$ and $T$. To ensure our simulation results are reproducible, we set the seed to 1234 in R programming.
First, we consider the one-factor case in Model ((ref)), i.e., $r=1$ is used in this section. The factor process is generated by \[f_t=\phi f_{t-1}+\eta_t,\quad \eta_t\overset{i.i.d.}{\sim} N(0,1),\quad t=1,...,T,\] where we set $\phi=0.5$. We consider the dimensions $N=50,100,150,300$ and $500$, with sample sizes $T=200,500,800,1000$, and $1200$ for each $N$. Each loading value is generated independently from the uniform distribution $U(-2,2)$, and the $N$-dimensional loading vector is then re-normalized to have an $l_2$-norm strength of $\sqrt{N}$. For each configuration of $(N,T)$, {we choose the sparsity $s=T/10$}, and we randomly generate $s$ integers from $\{1,...,T\}$ and keep those $s$ elements in $\underline{{\mathbf f}}_1=(f_1,...,f_T)'$ to be nonzero, re-normalizing the sparse vector $\underline{{\mathbf f}}_1=(f_1,...,f_T)'$ to have unit variance. We consider two scenarios for the idiosyncratic terms ${\mathbf e}_t$'s in each experiment: \[(1)\,\, {\mathbf e}_t\overset{i.i.d.}{\sim} N({\bf 0},{\mathbf I}_N)\quad\text{and}\quad (2)\,\, {\mathbf e}_t=\boldsymbol{\Phi}{\mathbf e}_{t-1}+\mbox{\boldmath$\varepsilon$}_t,\mbox{\boldmath$\varepsilon$}_t\overset{i.i.d.}{\sim}N({\bf 0},{\mathbf I}_N), t=1,...,T,\] where $\boldsymbol{\Phi}$ is a diagonal matrix with the diagonal elements generated from $U(0.5,0.9)\cup U(-0.9,-0.5)$. A total of 500 replications are used throughout the experiments.
We first study the accuracy of the estimated factors by defining the estimation error in each replication as
{Table (ref) reports the average estimation errors over 500 replications under two scenarios for the idiosyncratic components. {\it Proposal} refers to our proposed method, while {\it Lasso} denotes the Lasso-regularized approach applied to the factor process, following a similar procedure to Algorithm 1 in kristensen2017diffusion, with the tuning parameter set to $\psi_T = 0.6f_0$, where $f_0$ is the minimum nonzero component (in absolute value) of $\underline{\bf f}_1$. As shown in Table (ref), for each scenario, the estimation error decreases as the dimension $N$ increases for a fixed sample size $T$. In contrast, for fixed $N$, the error remains relatively stable or slightly increases as $T$ grows, which is consistent with the theoretical rate $\sqrt{\log(T)/N}$—a term that dominates $1/T$ for large $T$—as established in Theorem (ref). These results are consistent with our asymptotic theory. Moreover, the results clearly indicate that the Lasso procedure tends to incur additional estimation errors across all configurations of $(N, T)$. This highlights the advantage of the proposed hard-thresholding technique. }
Next, we study the accuracy of the empirical recovery (ER) of the sparse factors. Let $S$ be the set of true indexes of the $s$ nonzero elements in $\underline{{\mathbf f}}_1$, and $\widehat S$ the set of indices of the $\widehat s$ nonzero elements in $\underline{\widehat {\mathbf f}}_1$. For each replication, define the empirical recovery rate of the non-sparse indices as
where $\#\{S\cap\widehat S\}$ is the cardinality of the intersection between the estimated indices and the true ones. Table (ref) presents the empirical recovery rate of the sparsity in the factor process using the measure defined in ((ref)). {Similar to the comparison in Table (ref), we also compare the proposed method with the Lasso approach. From Table (ref), we observe that empirical accuracy improves as the dimension $N$ increases for each fixed sample size $T$, which is consistent with the asymptotic results presented in Section (ref). In addition, the Lasso procedure performs comparably to the proposed method in recovering the sparsity structure of the factor process, with the proposed method only slightly outperforming Lasso in certain cases. Taken together with the results in Table (ref), this suggests that the hard-thresholding technique provides an improvement in accuracy while remaining easy to implement. }
Finally, we evaluate the distribution of the estimated loadings as described in Theorem (ref). For simplicity, we consider the case when the noises are $i.i.d.$ and plot the empirical histogram of the estimated loading associated with the first series in Figure S.2 of the Supplement. The variance of the normal curve is based on the average of the estimated variance across 2000 replications. From Figure S.2, we see that the estimators behave closely to normal, which is in agreement with our asymptotic theory.
In this section, we further verify the efficacy of the proposed algorithms when there are multiple factors. The number of factors is set to be $r=3$, and the factors are generated by \[{\mathbf f}_t=\boldsymbol{\Psi}{\mathbf f}_{t-1}+\boldsymbol{\eta}_t,\boldsymbol{\eta}_t\overset{i.i.d.}{\sim}N({\bf 0},{\mathbf I}_r),t=1,...,T,\] where $\Psi=\mbox{diag}(0.5,-0.6,0.7)$ is a diagonal matrix. We consider the dimensions $N=50,100,150,200$, and $300$, with the sample size $T=100,200,300,500$, and $800$ for each $N$ in this experiment. For each configuration of $(N,T)$, we set the sparsity parameters to $s_1=s_2=s_3=T/10$. We first randomly generate $s_1$ indexes from $\{1,...,T\}$ to form an index set $S_1$, such that we only keep the corresponding $s_1$ time points of $\underline{{\mathbf f}}_1$ and the remaining ones are set to zero. The index set $S_2$ is formed by generating $s_2$ indexes $\{1,...,T\} \setminus S_1$, keeping only those $s_2$ locations of $\underline{{\mathbf f}}_2$. the sparse $\underline{{\mathbf f}}_3$ is obtained by repeating the above procedure. Then each $\underline{{\mathbf f}}_j$ is normalized to have unit variance. For the generation of the loading matrix $\boldsymbol{\Lambda}\in R^{N\times r}$, we first generate an $N\times r$ matrix ${\mathbf M}$ with elements independently generated from $U(-2,2)$, we then perform a singular-value decomposition on ${\mathbf M}$, and the left singular matrix is ${\mathbf U}$. The loading $\boldsymbol{\Lambda}$ is taken as ${\mathbf U}$ multiplied by $\sqrt{N}\mbox{diag}(3,2,1)$ on its right. We also consider the two scenarios for the idiosyncratic terms as in Section (ref). A total of 500 replications are used throughout the experiments.
Now, we first study the estimation accuracy of the factor processes. We define the measure of the errors of the estimated factors as
The average estimation errors of the factors over 500 replications are reported in Table S.I of the Supplement. From Table S.I, we observe a similar pattern to that in Section (ref). For each fixed sample size $T$, the error decreases as the dimension $N$ increases, which is in agreement with our asymptotic theory in Theorem (ref), regardless of whether the idiosyncratic terms are i.i.d. or dynamically dependent.
Furthermore, we study the empirical recovery (ER) of the sparsity in the factor processes. Similar to the measure in ((ref)), we define
{where $s=T/10$} and $\widehat S_j$ is the the estimated nonzero locations in $\underline{\widehat {\mathbf f}}_j$ for $j=1,2$, and $3$. Table S.II of the Supplement reports the empirical accuracy of the sparsity locations in the factor processes when $r=3$. From Table S.II, we see that the pattern is also similar to the case when $r=1$ in Section (ref). The empirical results are in line with our asymptotic ones in the sense that the estimation accuracy will increase as the dimension increases for each fixed $T$.
Although the ratio-based method of ((ref)) for determining the number of factors has been shown to be valid in many previous studies, we conduct an auxiliary experiment to verify its efficacy under our setting. Table S.III of the Supplement reports the empirical probabilities of $P(\widehat r=r)$ in 500 replications under the aforementioned setting. From Table S.III, we see that the ratio-based method performs well, with most of the empirical probabilities being close to one. This is understandable since all the factors are strong ones. Similar results can be found in gaotsay2022.
In this section, we study the estimation accuracy of the information criterion discussed in Section (ref) and Remark (ref) in estimating the sparsity of the factors. For simplicity, we only consider the case when $r=1$, but similar results can be obtained for $r>1$. The generation of the factors and the data is similar to those in Section (ref), but we only keep the largest $s=T/10$ elements of the factor process in absolute value for each sample size $T$. We consider the dimension $N=50,150,200$, and $300$, and the sample size $T=100,200,300,500$, and $800$ for each $N$ in this section. For each configuration of $(N,T)$, we set the number of partitions $J=1$ for simplicity and set $N_1=N_2=N/2$ in the partition. Since $s=T/10$ is diverging with the sample size $T$, we will use the information criterion defined in ((ref)) to estimate the sparsity. 500 replications are used throughout the experiment. Table S.IV of the Supplement reports the Empirical probabilities (EP) of determining the sparsity parameter using the information criterion in ((ref)) when the number of factors $r=1$. The empirical probabilities are calculated based on the 500 experiments. From Table S.IV we see that that the proposed information criterion estimates the sparsity parameter effectively. While dynamic dependence in idiosyncratic terms reduces accuracy, its performance generally improves with increasing $N$ for each fixed $T$, which is consistent with Theorem (ref).
In this section, we estimate the sparse latent risk factors across the time horizons in daily returns of individual stocks. The nonzero factors over certain time period can provide us one way to bridge time and certain type of risks in the financial market. The daily returns are downloaded from the CRSP daily security database and adjusted for dividend and stock splits. The data set used is the same as that in pelger2019 and consists of the daily stock returns for the balanced panel of S&P 500 stocks from January 1st 2004 to December 31st 2016. The daily interest rates from Prof. Kenneth French's website (\url{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}) are used to adjust the daily returns of individual stocks. The full data set is also available at \url{https://mpelger.people.stanford.edu/data-and-code}, where only the stocks with returns available for the full-time horizon are included, leaving us with a panel of $N=332$ and $T=3273$.
Similar to pelger2019, we group these 332 stocks into 14 categories by industry sector, as shown in Table S.V of the Supplement. The data encompasses a wide range of industry sectors within the stock market, including Oil, Finance, Electricity, Technology, Food, Manufacturing, Pharma & Chemicals, Primary Manufacturing, Machinery, Health, Transportation, Trade, Services, and Mining.
We first determine the number of factors for the centered data ${\mathbf X}$ using the eigenvalue-ratio method in ((ref)). Figure S.3 of the Supplement plots the ratios of eigenvalues of ${\mathbf X}{\mathbf X}'$ and shows that the largest gap between the eigenvalues occurs between $\widehat\lambda_1$ and $\widehat\lambda_2$, implying that the number of factors is $\widehat r=1$. Next, we apply Algorithm (ref) and the cross-validation method in Section (ref) with the number of random partitions $J=10$, where the seed number is set to be 1234 in R. {For the sparsity of the factors over the time horizon, we conduct a grid search over the time interval $s\in[\lceil \frac{T}{18} \rceil-50,\lceil \frac{T}{12} \rceil+50]=[132,232]$, where $\lceil x \rceil$ denotes the smallest integer that exceeds $x$. Figure (ref) plots the information criterion $IC(s)$ defined in ((ref)) and the estimated $\widehat s=211$. This suggests that there are 211 days during which stock returns were more strongly influenced by significant systematic risks. We also searched over the neighborhoods of \( T/10 \) and \( T/6 \), and found that the minimum of \( IC(s) \) in these regions is greater than \( IC(211) \). Moreover, the entire path of \( IC(s) \) increases after \( s = 211 \), indicating that \( \widehat s = 211 \) is a suitable estimate. It is important to note that this does not imply the absence of systematic risk on the remaining days; rather, the method is designed to highlight periods with stronger and more dominant risk signals. Compared to the full time span of $T=3273$, this represents a substantial reduction in temporal dimensionality, thereby simplifying subsequent analysis.}
We then plot the estimated sparse factor returns in Figure (ref). Similar to the example in the {\it Introduction}, we find that ${\sum_{t=1}^{3273}(f_t^s)^2}/{\sum_{t=1}^{3273} f_t^2} = 62.5\%$, indicating that the 211 nonzero factor values—representing only 6.4% of the total 3273 time points—account for 62.5% of the variance in the non-sparse factor process ${\mathbf f}_t$. Moreover, the plot reveals approximately three distinct clusters of time periods associated with significant systematic risks affecting the stock market. The first cluster is centered around 2008, driven by the 2007--2008 financial crisis. The second and third clusters are closely related, both associated with the European sovereign debt crisis. The most prominent cluster is centered around 2008, driven by the 2007--2008 financial crisis. This period saw severe disruptions in the financial markets, with significant declines in stock returns across the board. For the two closely related clusters associated with the European sovereign debt crisis, the first cluster corresponds to the initial phase of the crisis in 2010, while the second cluster relates to the intensification of the crisis in 2011 and 2012.
Next, we estimate the loadings of the 332 stock returns given the sparse common factors. We present the estimated loadings, indicating the level of dependence of 14 different industry sectors on the common risk factor in Figure (ref). The loadings quantify the extent to which each sector is influenced by these common factors. Here are the detailed observations based on the average dependence values: {\bf Oil Sector}: Exhibits moderate dependence on common risk factors, with an average loading of 0.0125. {\bf Finance Sector}: Shows significant dependence, with an average loading of 0.0165. {\bf Electricity Sector}: Demonstrates relatively low dependence, with an average loading of 0.0062. {\bf Technology Sector}: Reflects moderate dependence, with an average loading of 0.0097. {\bf Food Sector}: Indicates low dependence on common factors, with an average loading of 0.0046. {\bf Manufacturing Sector}: Displays moderate dependence, with an average loading of 0.0116. {\bf Pharma & Chemicals Sector}: Shows relatively low dependence compared to other sectors, with an average loading of 0.0070. {\bf Primary Manufacturing Sector}: Exhibits moderate dependence, with an average loading of 0.0133. {\bf Machinery Sector}: Demonstrates moderate dependence, with an average loading of 0.0109. {\bf Health Sector}: Indicates relatively low dependence, with an average loading of 0.0061. {\bf Transportation Sector}: Shows moderate dependence, with an average loading of 0.0099. {\bf Trade Sector}: Reflects moderate dependence, with an average loading of 0.0084. {\bf Services Sector}: Exhibits moderate dependence, with an average loading of 0.0111. {\bf Mining Sector}: Indicates relatively low dependence on common risk factors, with an average loading of 0.0079.
Overall, the Finance sector shows the highest average dependence on common risk factors, while sectors like Food, Health, and Electricity exhibit lower average dependence. These variations highlight the differing levels of sensitivity across industry sectors to common economic and market risk factors.
To further identify the systematic factors over time, referred to as time factors, Table S.VII in Sec. D of the Supplement provides a comprehensive list of dates with significant systematic risk factors, the reasons for these factors, and their associated time factors. Each entry explains the specific event and its impact on stock returns, summarized according to the descriptions of the time factors in Table S.VI of the Supplement. These reasons are extracted from the daily reports on {\it CNN Money} (\url{www.money.cnn.com}) after the market closes on each trading day. This website was shut down in late 2018 and its content was merged into the main CNN Business section.
We also plot the frequency charts of the nine factors over the time horizon in Figure (ref). The most frequently mentioned factor is market sentiment, highlighting its dominant role in driving stock price fluctuations. Economic indicators and government policies follow closely, indicating the significant impact of macroeconomic data and policy decisions on market behavior. Company-specific factors also play an important role, affecting individual stock performance and, by extension, broader market trends. Factors related to Europe, such as economic conditions and crises, have a moderate influence, reflecting the interconnectedness of global markets. Oil price fluctuations, China's economic activities, and global events are also notable contributors, underscoring the importance of global economic dynamics. Credit risk, while mentioned less frequently, remains a critical factor during periods of financial instability. Further implications of the results are discussed in Section B.1 of the Supplement for brevity.
Sparse factor modeling aims to extract meaningful insights by identifying a subset of relevant variables, enhancing interpretability while reducing computational complexity. This study presents a novel perspective by assuming that factor processes exhibit sparsity over time rather than imposing sparsity on loadings. This approach is particularly suited to scenarios where systematic co-movements are significant only during specific periods, such as financial crises or policy changes.
This paper introduces approximate factor models and a novel sparse asymptotic PCA (APCA) method, extending the asymptotic PCA framework of connor1986performance,connor1988risk to high-dimensional settings. Theoretical properties of the proposed method, including the cross-sectional cross-validation approach for estimating sparsity, are rigorously established. Through simulation and empirical studies, the method is shown to be effective in delivering interpretable results and linking co-movements to specific events. This work not only advances the statistical methodology of sparse factor modeling but also establishes a meaningful connection between time-specific events and risk factors in economic and financial systems.