Large-dimensional Factor Analysis with Weighted PCA

Large-dimensional Factor Analysis with Weighted PCA

Introduction

Large-dimensional factor analysis is an effective tool for extracting latent factors from noisy observations and is widely used in a variety of fields including economics, finance, and genomics among others. More specifically, consider the setting of panel data where observations, denoted by $\bfm X=(x_{it})_{1\le i\le N, 1\le t\le T}$, are obtained from $N$ units across $T$ time points. The factor model expresses $\bfm X$ as

where $\bfm F\in\mathbb{R}^{T\times r}$ and $\bfm L\in\mathbb{R}^{N\times r}$ represent the latent factors and factor loadings, respectively, and $\bfm E\in\mathbb{R}^{N\times T}$ the centered idiosyncratic noise. Here $r$ is much smaller than both $N$ and $T$ so that the signal lies in a low-dimensional space. In general, (ref) is referred to as the approximate factor model chamberlain1982arbitrage. It is also called the strict factor model if $E_{it}$'s are assumed to be uncorrelated anderson1956statistical. It is clear that $\bfm L$ and $\bfm F$ are only identifiable up to scaling and rotation. It is thus customary to assume that $\bfm F^\top\bfm F/T \to_p \bfm I_r$ as $T\to\infty$.

Seminal works such as connor1986performance,connor1988risk,stock2002forecasting, bai2002determining, bai2003inferential have established Principal Component Analysis (PCA) as the most popular approach for factor models (ref). See bai2008large for a comprehensive review. In spite of the popularity and successes of PCA, its effectiveness for large-dimensional factor analysis should not be taken for granted. To fix ideas, let us focus on estimating the column space, denoted by $\bfm U$, of $\bfm L$. PCA estimator, $\widehat{\bfm U}^{\rm PC}$, for $\bfm U$ can be viewed as the eigenspace of $\bfm X\bfm X^\top$. Intuitively, $\widehat{\bfm U}^{\rm PC}$ can only be expected to perform well if the eigenspace, $\bfm U^{\rm PC}$, of $\mathbb{E}[\bfm X\bfm X^\top]=\bfm L\bfm F^\top\bfm F\bfm L^\top + \mathbb{E}[\bfm E\bfm E^\top]$ is identical or at least close to $\bfm U$. The former ($\bfm U^{\rm PC}=\bfm U$) is true, for example, if $E_{it}$'s are independent and identically distributed; whereas the latter ($\bfm U^{\rm PC}\approx\bfm U$) can hold, more generally, if the factors are strong in that the first term $\bfm L\bfm F^\top\bfm F\bfm L^\top$ dominates the second, $\mathbb{E}[\bfm E\bfm E^\top]$. But what happens if the factors are not strong and $E_{it}$'s are cross-sectionally and temporally dependent? Are there situations where we can do better than the usual PCA for factor analysis? These are the questions we aim to address in this paper.

More specifically, we shall consider a weighted version of PCA: estimating $\bfm U$ by the eigenspace, $\widehat{\bfm U}_\bfm Q$, of $\bfm X \bfm Q\bfm X^\top$ for an appropriately chosen weight matrix $\bfm Q$. The usual PCA can be viewed as a special case with $\bfm Q=\bfm I_T$. More generally, the weight matrix offers a mechanism to re-balance the signal $\bfm L\bfm F^\top\bfm Q\bfm F\bfm L^\top$ against the “noise” $\mathbb{E}[\bfm E\bfm Q\bfm E^\top]$. Depending on the nature of the factors and noise, it is plausible that eigenspace, $\bfm U_{\bfm Q}$, of $\mathbb{E}[\bfm X\bfm Q\bfm X^\top]$ can be close to $\bfm U$ while $\bfm U^{\rm PC}$ is not. In particular, taking $\bfm Q=T^{-1}\bfm F\bfm F^\top$ yields $\bfm L\bfm F^\top\bfm Q\bfm F\bfm L^\top = T^{-1}\bfm L(\bfm F^\top\bfm F)^2\bfm L^\top\to T\bfm L\bfm L^\top$. On the other hand, $\mathbb{E}[\bfm E\bfm Q\bfm E^\top]=\textsf{PTr}_{T}\left\{\textsf{Cov}[\textsf{vec}(\bfm E)](\bfm I_N\otimes T^{-1}\bfm F\bfm F^\top)\right\}$ which is of the order $1$ and thus can be dominated by $\bfm L\bfm F^\top\bfm Q\bfm F\bfm L^\top$ even if $\mathbb{E}[\bfm E\bfm E^\top]$ is not of smaller order than $\bfm L\bfm F^\top\bfm F\bfm L^\top$. Here $\textsf{PTr}_{T}$ stands for the partial trace taken over the time index, i.e., for any square matrix $\bfm A\in \mathbb{R}^{pT\times qT}$, $\textsf{PTr}_{T}\left\{\bfm A\right\}:=\sum_{t=1}^T\left(\bfm I_p\otimes{\boldsymbol e}_t^\top\right)\bfm A \left(\bfm I_q\otimes{\boldsymbol e}_t\right)$, where $\left\{{\boldsymbol e}_t\right\}_{t=1}^T$ is the standard basis in $\mathbb{R}^T$ and $\otimes$ denotes the Kronecker product. Of course this specific choice of the weight matrix is too idealistic and infeasible in practice since it requires knowing $\bfm F$ a priori. Nonetheless, there are numerous situations where simple and agnostic choices of $\bfm Q$ can lead to superior performance over PCA.

See (ref) for further discussion about these two examples. To better understand the operating characteristics of weighted PCA, we shall derive non-asymptotic estimation error bounds and inferential theory for the estimated loading matrix and factors under general weighting schemes. These bounds pinpoint the effect of weighting and explain how it can provide an effective tool to address two of the most pressing challenges in large-dimensional factor analysis: weak factors and complex cross-sectional dependence.

Earlier studies for large-dimensional factor analysis have focused on strong factors in that $\bfm L^\top\bfm L/N^\alpha$ tends to a positive definite limit for some $\alpha\ge 1$. See bai2008large for a survey. The shift toward weak factors can be traced back at least to onatski2012asymptotics, who shows that $\alpha>0$ is necessary for the consistency of PCA. Subsequent advances establish consistency and asymptotic normality under weaker regimes uematsu2022estimation, bai2023approximate, choi2024high, fan2024can. In particular, choi2024high established the consistency and asymptotic normality of PC estimator for $\alpha\in(0,1)$, and fan2024can further relaxed the factor strength assumption for normality to logarithmic order. These results indicate that PCA-based estimates remain consistent and are asymptotically normal for weaker factors, but under much more restrictive assumptions about the dependence among the entries of $\bfm E$ that can be relaxed when using weighted PCA. In addition, we shall see that weighted PCA can be consistent under suitable conditions even when $\bfm L^\top\bfm L\rightarrow 0$.

As noted before, PCA-based estimates are expected to perform well when $\bfm E$ has independent and identically distributed entries. The potential limitations of PCA-based estimates when this assumption is violated have also come to light in the past few years. In particular, it has been observed that even if the entries of $\bfm E$ are independent but have different variances across units, the PCA-based estimates are suboptimal and can be much improved by the so-called HeteroPCA proposed in zhang2018heteroskedastic that iteratively re-estimates the eigenspace of $\bfm X\bfm X^\top$ by imputing its diagonal entries in order to reduce the bias caused by the cross-sectional heteroskedasticity. See also yan2021inference, agterberg2022entrywise. While HeteroPCA is designed specifically to address the cross-sectional heteroskedasticity, it remains unclear how to account for more general cross-sectional dependence.

Another related work is bai2013statistical who proposes to estimate factors and loadings via a weighted least-squares criterion with the optimal weight $\bfm W=\bfm \Sigma_{\rm C}^{-1}$ under the assumption that $\bfm \Sigma_{\rm C}:=\textsf{Cov}(\bfm E_{\cdot,t})$ for $t\in[T]$. The performance of their estimators hinges upon this structural assumption and other conditions to allow for consistent estimation of $\bfm \Sigma_{\rm C}$. In fact, they show that even an estimate $\widehat \bfm \Sigma_{\rm C}^{-1}$ that obeys $\big\|\widehat \bfm \Sigma_{\rm C}^{-1}- \bfm \Sigma_{\rm C}^{-1}\big\|\overset{p}{\rightarrow} 0$ may not be sufficient for valid inference. In practice, when $\bfm \Sigma_{\rm C}$ cannot be reliably estimated, it is unclear how to proceed with estimation and inference for the factors and loadings within that framework.

Our approach provides a more generally applicable solution to the above challenges by leveraging the potential temporal dependence among the factors. In particular, we shall show that with appropriately chosen weight matrices, weighted PCA can significantly reduce the estimation error of the factors and their loadings and allow for valid statistical inferences for much weaker factors than the usual PCA.

The choice of the weighting matrix is clearly of great practical importance. Motivated by the examples before, we shall consider choosing the weighting matrix from the following broad class of Toeplitz matrices: $$ {\mathscr Q}_T =\{\bfm Q=\textsf{Toeplitz}(\gamma_0,\gamma_1,\ldots, \gamma_{T-1}): \quad \gamma_0,\ldots,\gamma_{T-1}\ge0,\quad {\rm and}\quad \gamma_0+\cdots+\gamma_{T-1}=1\}. $$ Here we use $\textsf{Toeplitz}(\gamma_0,\gamma_1,\ldots, \gamma_{T-1})$ to denote a $T\times T$ symmetric matrix with diagonals being $\gamma_0$ and $t$-th super- and sub-diagonals being $\gamma_t$ for $t\ge 1$. Note that multiples of a weight matrix yield the same weighted PCA and thus the weighting matrices from both examples can be viewed as instances from ${\mathscr Q}_T$. Moreover, weight matrices from ${\mathscr Q}_T$ have an immediate statistical interpretation: $$\bfm X \bfm Q\bfm X^\top = \gamma_0\sum_{t=1}^T \bfm x_t\bfm x_t^\top + \gamma_1\sum_{t=1}^{T-1} (\bfm x_{t-1}\bfm x_t^\top+\bfm x_t\bfm x_{t-1}^\top)+\cdots,$$ where $\bfm x_t$ is the $t$-th column vector of $\bfm X$ proportional to a weighted linear combination of the covariance and auto-covariances of column vectors of $\bfm X$ if its entries are centered. This connects our work with recent developments in multivariate time series literature to utilize higher order auto-covariances for identifying latent factors. See, e.g., lam2011estimation,lam2012factor,chang2018principal, chang2024autocovariance.

Instead of fixing a weighting matrix a priori, we propose a cross-validation (CV) procedure to adaptively choose a weighting matrix from ${\mathscr Q}_T$ with theoretical guarantees. Conceptually, our method shares a similar spirit to common techniques for matrix completion where missing entries are imputed to recover low-rank structure via spectral methods, nuclear-norm penalties, or EM-type algorithms koltchinskii2011nuclear,xia2021statistical,mazumder2010spectral, chen2020noisy, choi2024matrix. In our setting, no entry of $\bfm X$ is truly missing, yet the same logic of “hiding information to gauge model fit” underlies our CV scheme for choosing $\bfm Q$ from ${\mathscr Q}_T$. We temporarily withhold a uniformly chosen subset of entries, fit the weighted PCA on the retained data, and select the Toeplitz weights that best predict the held-out values. We shall refer to this procedure as adaptive weighted PCA (AdaWPCA).

In summary, our contributions are threefold.

The rest of the paper is organized as follows. We shall introduce the general methodology for weighted PCA and study its performance in the next section. (ref) develops estimation bounds and inferential theory for factors and loadings. In (ref), we shall discuss Examples 1.1 and 1.2 in further detail. (ref) focuses on weight matrices from ${\mathscr Q}_T$ and studies how to adaptively choose a weighting matrix by CV. Our methodological and theoretical developments are further complemented by numerical experiments, both simulated and real, in (ref).

\paragraph{Notation} Let $\big\|\cdot\big\|_k$ denote either vector $k$-norm or the matrix norm induced by vector $k$-norm, and we omit the subscript when $k=2$ for brevity. We use $\mathbb{O}_{p,r}$ to denote the collection of $p\times r$ matrices with orthonormal columns, and write $\mathbb{O}_{p}:=\mathbb{O}_{p,p}$. In addition, we use $\left\| \cdot \right\|_{\rm F}$, $\big\|\cdot \big\|_{2,\infty}$ and $\big\|\cdot\big\|_{\sf max}$ to denote the Frobenius norm, $\ell_2\rightarrow\ell_\infty$ norm, and entry-wise max norm, respectively. For two matrices $\bfm A$ and $\bfm B$ of the same dimension, the Hadamard product $\bfm A\circ\bfm B$ is a matrix of the same dimension as the operands, with elements given by $[\bfm A\circ\bfm B]_{i,j}=A_{ij}B_{ij}$. For two non-negative sequences $a_N$ and $b_N$, we write $a_N \lesssim b_N$ ($a_N \gtrsim b_N$) or $a_N=O(b_N)$ if there exists some universal constant $C>0$ independent of $N$ such that $a_N \le Cb_N$ ($b_N \le C a_N$); $a_N \asymp b_N$ if $a_N \lesssim b_N$ and $b_N \lesssim a_N$ hold simultaneously; $a_N = o(b_N )$ or $b_N=\omega(a_N)$ or $a_N\ll b_N$ or $b_N\gg a_N$ if $b_N > 0$ and $a_N /b_N\rightarrow 0$. We use $C^k(D)$ to denote the space of all functions with $k$ continuous derivatives on a domain $D$. In addition, we shall write $n:=N\vee T$ for brevity.

Weighted PCA and Subspace Estimation

In this section, we introduce our weighted PCA framework for recovering the factors and loadings in the approximate factor model. To fix ideas, we shall focus our discussion on random factors. But it is worth pointing out that our approach and analysis can also be readily applied to deterministic factors. We opt for random factors merely for brevity. Recall that $\bfm F=\left[{\boldsymbol f}_1,\cdots,{\boldsymbol f}_T\right]^\top \in\mathbb{R}^{T\times r}$ , $\bfm L=\left[{\boldsymbol l}_1,\cdots,{\boldsymbol l}_N\right]^\top \in \mathbb{R}^{N\times r}$. Denote by $\bfm M=\bfm L\bfm F^\top$ the low rank component. Notice that $\bfm L$ and $\bfm F$ are only identifiable up to rotation and scale in terms of (ref). Following convention in the literature, we assume that

Under (ref), $\bfm F^\top\bfm F/T\rightarrow_p\bfm I_r$, and the loadings $\bfm L^\top\bfm L$ grow in accordance with the factor strengths $\lambda_i^2$s. Similar assumption is commonly adopted as an identifiability condition fan2024can. Throughout the paper, we allow $\lambda_i^2$s to vary from vanishing order to the order of $N$ which covers most practical scenarios.

Weighted PCA

The weighted PCA proceeds in two steps as outlined in (ref). We first compute the leading-$r$ left singular vectors, $\widetilde\bfm U_\bfm Q$, of $\bfm X\bfm Q\bfm X^\top $, which serve as estimators for eigenvectors of $\bfm L\bfm L^\top$. We then project our data onto the subspace spanned by $\widetilde\bfm U_\bfm Q$, as a denoising step, and obtain the truncated rank-$r$ SVD $(\widehat\bfm U_\bfm Q,\widehat\bfm \Sigma_\bfm Q,\widehat\bfm V_\bfm Q)$ of the projected data. Finally, we can construct the estimator for loadings as $\widehat\bfm L_\bfm Q:=\widehat\bfm U_\bfm Q\widehat\bfm \Sigma_\bfm Q$, and for factors as $\widehat\bfm F_\bfm Q:=\sqrt{T}\widehat\bfm V_\bfm Q$.

To gain intuition behind WPCA, recall that

where $\bfm N:=\bfm M\bfm Q\bfm E^\top+\bfm E\bfm Q\bfm M^\top+\bfm E\bfm Q\bfm E^\top$. The choice of $\bfm Q$ (after proper rescaling) shall meet two goals: (i) preserve the signal strength so that $\bfm L\bfm F^\top\bfm Q\bfm F\bfm L^\top$ remains at the same order as $\bfm L\bfm F^\top\bfm F\bfm L^\top$; (ii) reduce the magnitude of noise $\bfm N$ so that it is as small as possible. We shall restrict attention, without loss of generality, to symmetric weight matrices whose operator norm, i.e., the largest singular value, is $1$ unless otherwise indicated.

Clearly the usual PCA amounts to the choice of $\bfm Q=\bfm I_T$. To better understand the benefit of a different weight matrix $\bfm Q$ in (ref), denote $\left(\bfm U,\bfm \Sigma,\bfm V\right)$ as the rank-$r$ SVD of $T^{-1/2}\bfm M$. Following the explicit representation formula of the empirical spectral projectors from xia2021normal, we obtain

where the quantities $\sigma_r\left(\bfm F^\top\bfm Q\bfm F \right)$ and $\inf_{\zeta\in\mathbb{R}}\big\|\bfm U^\top \left(\bfm N-\zeta\bfm I_N\right)\bfm U_\perp\big\|$ can be interpreted as the signal strength and the noise level, respectively, for estimating the singular subspace $\bfm U$. Note the bound above is sharp in that it cannot be improved in general. Thus the SNR condition (and hence the quality of $\widehat\bfm U_{\bfm Q}$) depends on how well $\bfm Q$ aligns with both the factor matrix $\bfm F$ and the noise structure $\bfm E$. We shall now consider this in more detail.

Error Bounds for Subspace Estimation

We begin with the general covariance structure.

(ref) indicates the estimation error of $\widehat{\bfm U}_{\bfm Q}$ is determined by the signal term $\lambda_{r}^{2}\cdot \sigma_r\left(\bfm F^\top\bfm Q\bfm F \right)$ and $\textsf{GErr}_{U}(\bfm Q)$. More specifically, the first two terms of $\textsf{GErr}_{U}(\bfm Q)$ can be understood as the first and second order terms of the variance, and the third term represents the bias. In particular, if the noise is isotropic, i.e., $\bfm \Sigma_{\rm e}=\sigma^{2}\bfm I_{NT}$, then the bias term vanishes. For general $\bfm \Sigma_{\rm e}$, the bias $\big\|\bfm U^\top\textsf{PTr}_{T}\left\{\bfm \Sigma_{\rm e}(\bfm I_N\otimes \bfm Q)\right\}\bfm U_\perp\big\|$ can be reduced depending on the alignment between $\bfm \Sigma_{\rm e}$ and $\bfm Q$.

For illustrative purposes, consider the case when $\bfm \Sigma_{\rm e}$ is decomposable in that $\textsf{vec}\left(\bfm E\right)\sim {\cal N} \left(0,\bfm \Sigma_{\rm C}\otimes \bfm \Sigma_{\rm T}\right)$ where $\bfm \Sigma_{\rm C}\in\mathbb{R}^{N\times N}$ and $\bfm \Sigma_{\rm T}\in \mathbb{R}^{T\times T}$ are symmetric positive definite matrices. It is not hard to see that, for any $\bfm Q\in\mathbb{R}^{T\times T}$, $$\big\|\bfm U^\top\textsf{PTr}_{T}\left\{\bfm \Sigma_{\rm e}(\bfm I_N\otimes \bfm Q)\right\}\bfm U_\perp\big\|=\left|\textsf{Tr}(\bfm \Sigma_{\rm T}\bfm Q)\right|\big\|\bfm U^\top\bfm \Sigma_{\rm C}\bfm U_\perp\big\|.$$ In particular, the bias vanishes when $\textsf{Tr}(\bfm \Sigma_{\rm T}\bfm Q)=0$ or $\bfm \Sigma_{\rm C}$ is proportional to the identity. Following Theorem (ref), we get

Let $\kappa:=\lambda_1/\lambda_r$. Consider the case when $\big\|\bfm \Sigma_{\rm C}\big\|,\big\|\bfm \Sigma_{\rm T}\big\|, \big\|\bfm U^\top\bfm \Sigma_{\rm C}\bfm U_\perp\big\|, \kappa\asymp 1$. Then the error bound of (ref) can be simplified:

In the usual PCA, $\bfm Q=\bfm I_T$, the last term dominates whenever $T\gg \lambda_r^2N$. On the other hand, with appropriately chosen $\bfm Q$, it is possible to reduce the last term and thus yield improved estimates.

We can also derive similar error bounds for $\widehat{\bfm V}_{\bfm Q}$:

Different from the usual PC estimator, the estimation error of $\widehat\bfm V_\bfm Q$ relies on that of $\widehat\bfm U_{\bfm Q}$ as $\widehat\bfm V_\bfm Q$ is obtained from the projected data $\widetilde\bfm U_\bfm Q\widetilde\bfm U_\bfm Q^\top \bfm X$. In principle, our estimator would perform better than left singular vectors obtained by directly applying SVD on $\bfm X$. However, the bottleneck of estimating $\bfm V$, say, in $\big\|\cdot\big\|$, has a leading term of order $\lambda_r^{-1}\big\|\bfm \Sigma_{\rm T}^{1/2}\big\|\big\|\bfm \Sigma_{\rm C}^{1/2}\big\|$, hence we cannot obtain an order improvement compared to the PC estimator. On the other hand, we also conduct numerical experiments to illustrate the superiority of our method in practice in (ref).

Error Bounds in $\ell_2\to\ell_\infty$ Norm

To facilitate inferences about the factors and their loadings, it is often useful to derive error bounds in terms of $\ell_2\to\ell_\infty$ norm, i.e., $\big\|\widehat\bfm U_\bfm Q\widehat\bfm U_\bfm Q^\top -\bfm U\bfm U^\top \big\|_{2,\infty}$ and $\big\|\widehat\bfm V_\bfm Q\widehat\bfm V_\bfm Q^\top -\bfm V\bfm V^\top \big\|_{2,\infty}$. To this end, write

and $$\bar\rho:=\inf_{\zeta\in\mathbb{R}}\big\|\bfm \Sigma_{\rm C}-\zeta\bfm I_N\big\|_{1}.$$ Then we have

To fix ideas, consider in particular the scenario where $\big\|\bfm U\big\|_{2,\infty}\lesssim\sqrt{r/N}$, $\kappa_1,\kappa_2\asymp 1$, and $\sigma_r\left(\bfm F^\top\bfm Q\bfm F \right)\gtrsim T$. Then the error bound of (ref) becomes

under the signal-to-noise requirement:

Similar to before, the last term of the error bound (ref) is the bias term induced by the heteroskedasticity of $\bfm \Sigma_{\rm C}$. When $T\lesssim N$, it is dominated by the first two terms under the SNR condition (ref). It is also helpful to compare with some of the earlier results.

fan2024can studied the PC estimator for weak factor model with the additional assumption that $\bfm \Sigma_{\rm T}=\bfm I_T$. In particular, under the above setting, their Lemma 4 and Lemma 7 therein entail that, with probability at least $1-O\left(n^{-2}\right)$,

provided that

It is not hard to see that this coincides with our bound (ref) by taking $\bfm Q=\bfm I_T$. There is, however, a subtle difference between their SNR requirement and ours: the two are equivalent only when $T\lesssim N$, and (ref) is more stringent when $T\gg N$.

We can also compare the bound (ref) with those for HeteroPCA where $\bfm \Sigma_{\rm C}$ is assumed to be diagonal. To facilitate comparison, we consider $\bfm \Sigma_{\rm T}=\textsf{diag}\big(\sigma_{{\rm T},1}^2,\cdots, \sigma_{{\rm T},T}^2\big)$ and $\bfm \Sigma_{\rm C}=\textsf{diag}\big(\sigma_{{\rm C},1}^2,\cdots, \sigma_{{\rm C},N}^2\big)$. Denote by $\sigma_{\sf max}^2:=\big\|\bfm \Sigma_{\rm C}^{1/2}\big\|\big\|\bfm \Sigma_{\rm T}^{1/2}\big\|$ and $\widehat\bfm U^{\rm HPC}$ the HeteroPCA estimate of the left singular space. Then the results from yan2021inference and agterberg2022entrywise imply that with high probability,

provided that

Note that our bound (ref) and (ref) coincide under this particular setting. Yet our proposed weighted PCA approach is more generally applicable and allows for both cross-sectional and temporal dependence beyond heteroskedasticity. Indeed, our numerical experiments in (ref) suggest that the HeteroPCA might deteriorate significantly when there is cross-sectional dependence, whereas the weighted PCA remains unaffected.

Factor Analysis: Estimation and Inference

We now turn our attention to the estimation error and asymptotic normality of the estimated factors and loadings.

Error Bounds for Factors and Loadings

To start with, let $\bfm H_{U,\bfm Q}:=\widehat\bfm U^\top _\bfm Q\bfm U$, $\bfm H_{V,\bfm Q}:=\widehat\bfm V^\top _\bfm Q\bfm V$ and $\bfm R_{U,\bfm Q}:=\textsf{sign}\big(\widehat\bfm U^\top _\bfm Q\bfm U\big)$, $\bfm R_{V,\bfm Q}:=\textsf{sign}\big(\widehat\bfm V^\top _\bfm Q\bfm V\big)$. In addition, let $\bfm B\in\mathbb{R}^{r\times r}$ be the invertible matrix obeying $\bfm V=T^{-1/2}\bfm F\bfm B$ and define $\bfm R_{L,\bfm Q}:=\left(\bfm B^{-1}\right)^{\top }\bfm R_{V,\bfm Q}^\top $ and $\bfm R_{F,\bfm Q}:=\bfm B\bfm R^\top_{V,\bfm Q}$. Then we get

As an immediate consequence of (ref), (ref) and (ref), we can obtain the estimation error for $\bfm L$ and $\bfm F$. It is worth noting that $\lambda_r=\omega(1)$ is not necessary for consistency. Nonetheless, to simplify narratives, let us consider the regime $T\gtrsim N$, $\big\|\widehat\bfm U_\bfm Q\widehat\bfm U_\bfm Q^\top -\bfm U\bfm U^\top \big\|=o(1)$, $\lambda_r=\omega(1)$, $\kappa, \big\|\bfm \Sigma_{\rm T}\big\|,\big\|\bfm \Sigma_{\rm C}\big\|\asymp1$. Then under the assumptions of (ref) and (ref) we have

If either $\big\|\bfm U^\top\bfm \Sigma_{\rm C}\bfm U_\perp\big\|=0$ or $\bfm Q$ is chosen so that $\textsf{Tr}(\bfm \Sigma_{\rm T}\bfm Q)=0$, then the bias term in (ref) vanishes. Then, the factors converge at rate $\lambda_r ^{-1}$, which matches the result in fan2024can. For loadings, the convergence rate is $T^{-1/2}$. This is to be compared with the convergence rate of $T^{-1/2}+\lambda_r^{-1}N^{-1/2}$ for the usual PC estimator derived by fan2024can, suggesting a faster convergence rate for the weighted PCA when $T\gg \lambda_r^2N$.

Similarly, we can also relate the perturbation bound in (ref) to derive the estimation error in $\ell_2\to\ell_\infty$ norm for factors and loadings.

For simplification, consider the case when $T\gtrsim N$, $\big\| \bfm \Sigma_{\rm C}\big\|, \big\| \bfm \Sigma_{\rm T}\big\|\asymp 1$, $\lambda_r=\omega(1)$ and $\big\|\widehat\bfm U_\bfm Q\widehat\bfm U_\bfm Q^\top -\bfm U\bfm U^\top \big\|_{2,\infty}\lesssim\big\|\widehat\bfm U_\bfm Q\widehat\bfm U_\bfm Q^\top -\bfm U\bfm U^\top \big\|\big\|\bfm L\bfm \Lambda^{-1}\big\|_{2,\infty}$. Then the bounds in (ref) reduce to $$\big\|\widehat\bfm L_\bfm Q-\bfm L\bfm R_{L,\bfm Q}\big\|_{2,\infty}\lesssim \big\|\widehat\bfm U_\bfm Q\widehat\bfm U_\bfm Q^\top -\bfm U\bfm U^\top \big\|\big\|\bfm L\bfm \Lambda^{-1}\big\|_{2,\infty}\lambda_1.$$ In other words, the $\ell_2\to\ell_\infty$ errors of factor and loading estimates are of the same order as the subspace estimates after appropriate scaling. More importantly, (ref) provides a basis for us to develop an inferential theory for factors and loadings.

Asymptotic Normality

We are now in a position to establish the asymptotic normality of our factor and loading estimators. To this end, denote by $\bar\bfm U$ the leading singular vectors of $\bfm M\bfm Q\bfm M^\top $, and $\bar\bfm \Sigma$ the diagonal matrix containing singular values of $\bfm M\bfm Q\bfm M^\top$. Note that there exists some $\bar\bfm O\in\mathbb{O}_{r}$ such that $\bar\bfm U\bar\bfm O=\bfm U$. The following theorem presents the asymptotic normality for our estimate of the loadings.

Note that $\lambda_r\gg 1$ is not required for the consistency of $\widehat\bfm U_\bfm Q$ when $T\gg N$. But a diverging $\lambda_r$ is necessary for the asymptotic normality. Similar observations have been made earlier for PCA onatski2012asymptotics, choi2024high.

(ref) suggests that under suitable conditions, the weighted PCA estimate of $\bfm L$ may enjoy asymptotic normality under weaker SNR conditions than the usual PCA. To fix ideas, consider the setting when $\kappa_1,\kappa_2,\kappa,r,\big\|\bfm \Sigma_{\rm C}^{1/2}\big\|_1,\big\|\bfm \Sigma_{\rm T}^{1/2}\big\|\asymp 1$ and ignore logarithmic terms. Then the SNR condition in (ref) for the PC estimator (i.e., $\bfm Q=\bfm I_T$) becomes

As (ref) indicates, our loading estimate can be asymptotic normal under weaker SNR by choosing a weight matrix $\bfm Q$ so that the bias $\left|\textsf{Tr}(\bfm \Sigma_{\rm T}\bfm Q)\right|$ can be appropriately controlled. However, we shall now argue this condition is essentially optimal for PC estimates when the noise is heteroskedastic.

To see this, consider the rank-one model $\bfm X={\boldsymbol l}{\boldsymbol f}^\top+\bfm E$ with $\bfm \Sigma_{\rm C}=\textsf{diag}\big(\sigma_{{\rm C},1}^2,\cdots,\sigma_{{\rm C},N}^2\big)$, $\bfm \Sigma_{\rm T}=\bfm I_T$. For simplicity, we assume $\sigma_{\rm C}\lesssim \min_{i\in[N]}\sigma_{{\rm C},i}$ with $\sigma_{\rm C}:=\max_{i\in[N]}\sigma_{{\rm C},i}$. In addition, suppose that $\big\|\bfm U_\perp^\top\bfm \Sigma_{\rm C}{\boldsymbol u}\big\|\gtrsim \sigma_{\rm C}^2$ and $\lambda=\omega \big(\sigma_{\rm C}\sqrt{{n}/{T}}\log n\big)$, then we can show that

for some constant $c\in \mathbb{R}$, where ${\boldsymbol g}_L\sim N\big(0,\bfm \Sigma_{\rm C}\big)$ and with probability at least $1-\eta_F-O\left(e^{-c_0N}+n^{-9}\right)$,

Note that the bias ${\boldsymbol b}_L$ is non-vanishing unless $\lambda\gg \sigma_{\rm C}^2\sqrt{{T}/{N}}$. This immediately suggests the necessity of (ref).

To further illustrate the difference in SNR requirements between the usual PCA and weighted PCA, we conducted a simulation study with $N=100$ and $T=800$. (ref) shows the histogram of the estimation error for both estimates. It is clear that the PC estimate remains biased in this setting whereas the weighted PC estimate does not.

Similarly, the following theorem establishes the asymptotic normality of the estimated factors.

Note that the incoherence condition $\big\|\bfm L\bfm \Lambda^{-1}\big\|_{2,\infty}$ is commonly adopted in the literature fan2024can, choi2024high. In addition, (ref) can be viewed as a generalization of the results from fan2024can who considered the usual PCA estimate (i.e., $\bfm Q=\bfm I_T$) under the special case when $\bfm \Sigma_{\rm T}=\bfm I_T$.

Examples

To further illustrate the benefit of weighted PCA and the implication of our theoretical development, we now revisit the two motivating examples introduced in (ref). To fix ideas, we shall assume throughout this subsection that $\kappa,\kappa_1,\kappa_2,\big\|\bfm \Sigma_{\rm C}\big\|\asymp1$ for simplicity.

\paragraph{Multivariate Time Series Model.} Suppose the idiosyncratic noise $\bfm E$ is temporally uncorrelated and $\bfm F$ follows a multivariate autoregressive model that is independent of $\bfm E$:

where $\boldsymbol{\varepsilon}_t\sim {\cal N}(0,\bfm I_r-\bfm A\bfm A^\top )$. In this scenario, a practical choice is to take $\bfm Q=\sum_t(\bfm e_t\bfm e_{t-1}^\top+\bfm e_{t-1}\bfm e_t^\top)$, which leverages the first-order auto-covariance information.

In light of the results above, we get

This setting connects our work with a fast growing literature in high dimensional time series analysis. See, e.g, lam2011estimation,gao2022modeling,qiao2025weight. In particular, lam2011estimation and gao2022modeling considered the auto-covariance-based method, whilst more recently qiao2025weight proposed a weight-calibrated estimation method. All their theoretical results for estimating $\bfm U$ are restricted to the moderate-dimensional regime in that there exists some constant $\alpha\in(0,1]$ such that $N^{1-\alpha}\ll T\lesssim N$ with the factor strength $ \lambda_r\asymp N^{\alpha}$. They established the rate of $O_p\big(\sqrt{N^{1-\alpha}/{T}}\big)$, which is the same rate entailed by our result under the same regime. However, our result can go beyond this regime, as we do not require the stringent dimension requirement $N^{1-\alpha}\ll T\lesssim N$. In fact, in terms of (ref) we only require $\lambda_r\gtrsim \sqrt{{N}/{T}}$ and hence $\lambda_r$ can go to $0$ if $T\gg N$. This captures a broader regime than previous work. Moreover, to our best knowledge, there is no asymptotic normality result similar to ours that exists in this line of work.

\paragraph{Multivariate Functional Data.} Let $\bfm F$ be a smooth function so that $F_{t,k}=g_k(t/T)$, where $g_k\in C^1([0,1])$ for $k\in[r]$. For identifiability, we assume that

We consider $\bfm Q\in\mathbb{R}^{T\times T}$ such that $Q_{t,s}=\frac{1}{2B}\mathbb{I}\left(\left|t-s\right|\le B\right)$, where $B$ is the number of sub-diagonals that are nonzero.

As noted in (ref), requiring a divergent SNR for consistency is not intrinsic to factor models but a limitation of the PC estimator. The error rate (ref) and (ref) imply that under this assumption of functional data for $\bfm F$, we shall choose $B\ge \lambda_r^{-2}$ to obtain a sharp rate for the subspace and loadings estimation, when the factor signal is extremely small in that $\lambda_r\ll 1$.

We can also derive properties of the estimated factors under this setting which relates to functional PCA hall2006propertiesa, hall2006propertiesb, shang2014survey, zhou2022theory.

We can compare these results with those for functional PCA in the literature. Consider, for example, the setting from zhou2022theory where $r=O(1)$ and $T\gtrsim N$. In addition, it is assumed that $\lambda_k^2\asymp N$ for $k\in[r]$ and $\max_{k\in[r]}\left\{\big\|g_k\big\|_{\infty}\vee\big\|g_k^\prime\big\|_{\infty}\vee\big\|g_k^{\prime\prime}\big\|_{\infty}\right\}=O(1)$. By choosing the optimal tuning parameter and assuming suitable regularity conditions hold, zhou2022theory (Corollary 2 therein) showed that the functional PCA satisfies $$\big\|\widehat g_k^{\rm fPCA}-g_k\big\|_{L_2}^2=\int_{0}^1\left|\widehat g_k(t)-g_k^{\rm fPCA}(t)\right|^2 dt=O_p\left(\frac{1}{N}\right),\quad \forall k\in[r].$$ To obtain the above bound, a strong eigen-gap condition was imposed to ensure that we can differentiate between the different PCs. Under the same condition, we can also show that the rotation matrix $\bfm R_{F,\bfm Q}$ in (ref) obeys $\bfm R_{F,\bfm Q}=\bfm I_r+\bfm \Delta_F$ with $\big\|\bfm \Delta_F\big\|=O(1/T)$. This implies that $$\frac{1}{T}\left\| \widehat\bfm F_{\bfm Q}-\bfm F \right\|_{\rm F}^2=\sum_{k=1}^r\frac{1}{T}\sum_{t=1}^T\left[\widehat g_k\left(\frac{t}{T}\right)-g\left(\frac{t}{T}\right)\right]^2=O_p\left(\frac{1}{N}\right),$$ for finite $r$. However, different from the functional PCA, this holds without the strong factor assumption and without the boundedness condition on $\left\{g_k^{\prime\prime}\left(\cdot \right)\right\}_{k=1}^r$.

Adaptivity: CV

The choice of the weight matrix $\bfm Q$ is clearly essential to weighted PCA. In this section, we shall discuss how we can adaptively choose from a large class of weight matrices via CV. To this end, we shall assume in the rest of this section that $\left\{{\boldsymbol f}_t,t\in \mathbb{Z}\right\}$ is weakly stationary in the sense that $\mathbb{E} {\boldsymbol f}_t=0$, $\mathbb{E} {\boldsymbol f}_{t}{\boldsymbol f}^\top _{t^\prime}=\mathbb{E} {\boldsymbol f}_{t-t^\prime}{\boldsymbol f}^\top _{0}$ and $\mathbb{E} \big\|{\boldsymbol f}_t\big\|^2<\infty$ for $\forall t, t^\prime\in\mathbb{Z}$. Let

denote the auto-covariance at lag $\tau=1,\cdots,K$. Furthermore, we assume that

Note that (ref) can be viewed as a relaxed version of $\tau$-mixing condition with geometric decay. It holds, in particular, with $\bar\eta_F=o\left(n^{-C}\right)$ for any constant $C>0$ under suitable moment conditions on $\bfm F$. See e.g., han2020moment (Assumptions (A1)-(A4) therein).

In light of (ref), we shall consider the following class of weight matrices:

where

In particular, ${\boldsymbol{\gamma}}=\left(1,0,\cdots,0\right)$ corresponds to classical PCA, and ${\boldsymbol{\gamma}}=\left(0,\gamma_1,\cdots,\gamma_K,0,\cdots 0\right)$ corresponds to using aggregated information from lag-$1$ to lag-$K$ auto-covariance.

The signal strength under this setting can be characterized by the following quantity:

It is easy to see that $\sigma_r\left(\mathbb{E}\bfm M\bfm Q_{\boldsymbol{\gamma}}\bfm M^\top\right)=\sigma_r\left(\bfm L\left(\mathbb{E}\bfm F^\top\bfm Q_{\boldsymbol{\gamma}}\bfm F\right)\bfm L^\top \right)\top }\ge \lambda_r^2T \mu({\boldsymbol{\gamma}})$, which serves as a high probability lower bound for $\textsf{Signal}\left(\bfm Q_{\boldsymbol{\gamma}}\right)$ under (ref). A general lower bound for $\mu({\boldsymbol{\gamma}})$ depends on the structures of $\left\{\bfm \Gamma_\tau\right\}_{\tau=1}^K$. However, there are several scenarios that we can have an explicit lower bound on $\mu \left({\boldsymbol{\gamma}}\right)$. First, $\mu\left({\boldsymbol{\gamma}}\right)\gtrsim 1$ holds when $\gamma_0$ is close to either $0$ or $1$. In addition, when $\bfm \Gamma_1+\bfm \Gamma_1^\top $ has positive eigenvalues, where we can deduce that $\mu\left({\boldsymbol{\gamma}}\right)\gtrsim 1$ for any $\gamma\in[0,1]$.

Now we shall choose $\bfm Q_{\boldsymbol{\gamma}}\in {\mathscr Q}_{T,K}$, or equivalently, ${\boldsymbol{\gamma}}\in \Delta_{T,K}$ via CV: a random subset of the entries of $\bfm X$ is selected as the validation set and ${\boldsymbol{\gamma}}$ is set to the value that minimizes the validation error. More formally, let $\bfm \Omega$ be a $N$ by $T$ matrix with entries following i.i.d. Ber$(p_*)$ for a prespecified value of $p_\ast$ and $\Omega:=\left\{(i,t)\in [N]\times [T]:\Omega_{i,t}=1\right\}$ be the index set. For each $\Omega$, define the projection operator ${\cal P}_{\Omega}:\mathbb{R}^{N\times T}\rightarrow \mathbb{R}^{N\times T}$ as

In addition, let $\bfm \Omega_{\perp}:=\mathbf{1}_N\mathbf{1}_T^\top-\bfm \Omega $ and similarly define $\Omega_{\perp}$. Let $\widehat\bfm L_{\boldsymbol{\gamma}}^{\Omega}=\left[\widehat{\boldsymbol l}_{{\boldsymbol{\gamma}},1}^{\Omega},\cdots,\widehat{\boldsymbol l}_{{\boldsymbol{\gamma}},N}^{\Omega}\right]^\top \in\mathbb{R}^{N\times r}$ and $\widehat\bfm F_{\boldsymbol{\gamma}}^{\Omega}=\left[\widehat{\boldsymbol f}_{{\boldsymbol{\gamma}},1}^{\Omega},\cdots,\widehat{\boldsymbol f}_{{\boldsymbol{\gamma}},T}^{\Omega}\right]^\top \in\mathbb{R}^{T\times r}$ be our estimator using $p_*^{-1}{\cal P}_{\Omega}\left(\bfm X\right)$. Finally, we choose

where ${\cal G}_{K}=\left\{{\boldsymbol{\gamma}}_1,\cdots,{\boldsymbol{\gamma}}_M\right\}\subset \Delta_{T,K}$ is a grid set of $\Delta_{T,K}$ with cardinality $M$. See (ref) for details. It is not difficult to check that an $\epsilon$-cover of $\Delta_{T,K}$ requires only $M=O\big(\left(\log (1/\epsilon)\right)^K\big)$ grid points.

For technical reasons, we impose the following assumption for theoretical development, though our experience with the numerical experiments suggests that the CV approach works well under more general settings.

For simplicity, consider the case when $\kappa,r\asymp1$ and under the weak factor model with $\lambda_1/\left(\sigma_{\rm C}\sigma_{\rm T}\right) \asymp\log^{3/2} T$ and $T\asymp N\log^{2}T$. (ref) states that we can choose $\widehat{\boldsymbol{\gamma}}$ such that $\widehat\gamma_0=o(1)$ while avoiding those ${\boldsymbol{\gamma}}$ such that $\mu\left({\boldsymbol{\gamma}}\right)=o\big(\left(\log T\right)^{-1}\big)$. In practice, it is oftentimes not known whether the heteroskedasticity of noise or the autocorrelation structure of factors exists. In general, a desirable goal is to consistently choose ${\boldsymbol{\gamma}}$s with $\gamma_0=o(1)$ over those ${\boldsymbol{\gamma}}$s with $\gamma_0=1$ when autocorrelation is strong, and consistently choose ${\boldsymbol{\gamma}}$s with $\gamma_0=1$ over ${\boldsymbol{\gamma}}$s with $\gamma_0=o(1)$ when autocorrelation is weak. Nonetheless, it is unclear whether achieving this stronger consistency demands a new algorithm beyond the current CV scheme or simply more refined analysis, and we leave it for future work.

It is also worth pointing out that similar CV approaches have been adopted by zeng2019double, wei2020determining, jin2021factor. A nuanced difference between our goal and theirs is that they use CV to choose the number of factors $r$, and the inherent low-rankness can benefit analysis for delivering consistency of CV. In contrast, choosing the Toeplitz weights ${\boldsymbol{\gamma}}$ by CV brings out additional technical difficulty, which in turn hinges on a lower bound on singular subspace perturbation in the missing setting.

In practice, we can apply multiple random draws of training set $\left\{\bfm \Omega_{j}\right\}_{j=1}^{K_{\sf cv}}$. For each ${\boldsymbol{\gamma}}\in{\cal G}_K$, we can compute $\textsf{CV}_j\left({\boldsymbol{\gamma}}\right)$ and choose $$\widehat{\boldsymbol{\gamma}}=\operatorname*{arg\,min}_{{\boldsymbol{\gamma}}\in{\cal G}_K}\frac{1}{K_{\sf cv}}\sum_{j=1}^{K_{\sf cv}}\textsf{CV}_j\left({\boldsymbol{\gamma}}\right).$$ (ref) continues to hold for this procedure, and one would expect the resultant $\widehat{\boldsymbol{\gamma}}$ to benefit from the reduction of variance introduced by a single draw of $\bfm \Omega$.

When $r$ is unknown, there has been a line of research on estimating it from the data bai2007determining,lam2012factor,wei2020determining, fan2022estimating. In particular, we can determine the number of factors by the ratio-based estimator as $\widehat r:=\operatorname*{arg\,min}_{j\in[R]}\sigma_{j+1}\left(\bfm X\bfm X^\top\right)/\sigma_{j}\left(\bfm X\bfm X^\top\right)$ in lam2012factor, where $R$ is a reasonable upper bound of $r$ and can be taken as $R=N/2$ in practice. Alternatively, we can use CV with a $(K+1)$-dimensional grid set to choose $\left({\boldsymbol{\gamma}},r\right)\in {\cal G}_{K}\times [R]$ simultaneously.

Numerical Experiments

Simulated-data Analysis

We first consider the factor model in (ref) with rank $r=3$. The factors are assumed to be generated from the VAR model ${\boldsymbol f}_t=\bfm A{\boldsymbol f}_{t-1}+\boldsymbol{\varepsilon}_t$, where $\bfm A=\bfm O_1\bfm D\bfm O_2^\top$ with $\bfm D:=0.9\bfm I_r$ and $\bfm O_1, \bfm O_2$ being two random draws from $\mathbb{O}_{r}$, and $\boldsymbol{\varepsilon}_t\sim {\cal N}(0,\bfm I_r-\bfm A\bfm A^\top )$.

\paragraph{Estimation of Singular Vectors $\bfm U$ and $\bfm V$.} We consider two regimes:

We set $N\in\left\{100,200\right\}$ and $T\in\left\{250, 500, 750,1000\right\}$. For ease of presentation, we consider $K=1$ and ${\boldsymbol{\gamma}}=(\gamma,1-\gamma,0,\cdots,0)\in{\cal G}_1$, for which we use a single parameter $\gamma\in \bar{\cal G}:=\left\{0,1/9,2/9,\cdots,1\right\}$ to represent the Toeplitz weights. Note that all results shown are based on adaptive choice of the Toeplitz weights $\gamma$ by CV with $K_{\sf cv}=10$, and all results are averaged from 100 simulation replicates.

Setting (1) is a realistic and challenging setting where non-zero off-diagonal entries in $\bfm \Sigma_{\rm C}$ exist. As shown in (ref), AdaWPCA outperforms the other two in all different $(N,T)$ settings regarding the estimation of $\bfm U$ and $\bfm V$.

Setting (2) is the preferable setting for HeteroPCA, which is specifically designed for diagonal-only heteroskedasticity. As shown in (ref), AdaWPCA has almost comparable performance to PCA and HeteroPCA in terms of estimation of $\bfm U$, while it outperforms the other two regarding the estimation of $\bfm V$.

\paragraph{Inference of Loadings $\bfm L$ and Factors $\bfm F$.} We consider the regime $\bfm \Sigma_{\rm T}=\bfm I_T$ and $\bfm \Sigma_{\rm C}=\textsf{diag}\big(\left\{\omega_i\right\}_{i=1}^N\big)+0.6\left(\mathbf{1}_N\mathbf{1}_N^\top -\bfm I_N\right)$ with $\omega_i\overset{i.i.d.}{\sim} \textsf{Unif}[1,20]$, and set $\bfm Q=\sum_{t=1}^{T-1}\big({\boldsymbol e}_{t}{\boldsymbol e}^\top _{t+1}+{\boldsymbol e}_{t+1}{\boldsymbol e}^\top _{t}\big)$. In each simulation for loadings, we calculate $\bfm \Sigma_{L,N/2}^{-1/2}\big(\widehat\bfm L_{\bfm Q}-\bfm L\bfm R_{L,\bfm Q}\big)_{N/2,\cdot}^\top $ as in (ref) and plot Q-Q plot and histogram of its first dimension. Each data point is averaged from $500$ simulation replicates. Similarly for factors, we calculate $\bfm \Sigma_{F,T/2}^{-1/2}\big(\widehat\bfm F_{\bfm Q}-\bfm F\bfm R_{F,\bfm Q}\big)_{T/2,\cdot}^\top $ as in (ref) and plot Q-Q plot and histogram of its first dimension. The same setting applies for (ref). From (ref)-(ref), we can see that the standard normal density curve provides good approximations to the normalized histograms, certifying our inferential theories in (ref).

\paragraph{CV for Optimal Toeplitz weight $\gamma$.} We consider two regimes:

Similar to the estimation part, we consider $K=1$ and $\gamma\in \bar{\cal G}=\left\{0,1/9,2/9,\cdots,1\right\}$. In particular, let $R(\widehat\gamma)$ denote the rank of $\widehat\gamma$ according to the CV error for all $\gamma$s in $\bar{\cal G}$. We then calculated the proportions of 100 simulations in which CV choose the top three $\gamma$s in $\bar{\cal G}$ (i.e., $R(\widehat\gamma)\le 3$) and not choosing the bottom three $\gamma$s in $\bar{\cal G}$ (i.e., $R(\widehat\gamma)\ge 8$). The results are summarized in (ref). Notably, the result indicates that CV may not always choose the best $\gamma$ especially when off-diagonal heterogeneity exists, but it can avoid the worst choices of $\gamma$s in most cases. These observations coincide with (ref).

Real-data Analysis

We consider two datasets\footnote{The datasets are publicly available and updated in a timely manner at \url{https://www.stlouisfed.org/research/economists/mccracken/fred-databases.}} from the FRED (Federal Reserve Economic Data): FRED-MD mccracken2016fred and FRED-QD mccracken2020fred, which are monthly and quarterly macroeconomic databases consisting of U.S. indicators. FRED-MD dataset contains $126$ variables spanning from 1959-01-01 to 2025-01-01. FRED-QD dataset contains $246$ variables spanning from 1959-Q1 to 2025-Q1. To obtain a balanced panel, we removed the variables which have more than $5\%$ missing values and kept the time periods without missing values. This leads to a complete data matrix $\bfm X_{\rm MD}\in\mathbb{R}^{N\times T}$ with $N=121$ and $T=777$ for FRED-MD and $\bfm X_{\rm QD}\in\mathbb{R}^{N\times T}$ with $N=207$ and $T=259$\footnote{The preprocessing procedure largely follows the example given by \url{https://github.com/cykbennie/fbi/blob/master/doc/factor_fred.html}, except that we do not perform data transformation.}. In the following, we use $\bfm X$ to denote either $\bfm X_{\rm MD}$ or $\bfm X_{\rm QD}$.

Since we do not have access to true factors and loadings, we adopt the similar idea in CV to showcase the validity of our method. In particular, we randomly mask the entries of $\bfm X$ by $\bfm \Omega^*\in\left\{0,1\right\}^{N\times T}$ with $\Omega_{i,t}^*\overset{i.i.d.}{\sim}\text{Bern}\left(q_{\sf tr}\right)$ with $q_{\sf tr}\in\left\{0.7,0.8,0.9\right\}$. Then we treat $\bfm X\circ\bfm \Omega^*$ as our observed data and use $\bfm X\circ\left(\mathbf{1}\mathbf{1}^\top-\bfm \Omega^*\right)$ as the testing sample. We first determine the rank $r$ by using the ratio-based estimator on the full data matrix $\bfm X$, which gives us $\widehat r=1$ for both FRED-MD and FRED-QD datasets. To incorporate the case of over-specified rank, we apply three methods with rank $r\in \left\{1,2,3\right\}$. The performance is measured in terms of the relative reconstruction error defined as $\sum_{(i,t)\in\Omega_\perp}\big(X_{i,t}-\widehat{\boldsymbol l}_{i}^{\Omega^*\top }\widehat{\boldsymbol f}^{\Omega^*}_{t}\big)^2/\sum_{(i,t)\in\Omega_\perp}X_{i,t}^2$, where $(\widehat\bfm L^{\Omega^* },\widehat\bfm F^{\Omega^*})$ are corresponding estimators based on $\bfm X\circ \bfm \Omega^*$. All reported results shown in (ref) are averaged out of $100$ replicates in terms of the randomness of $\bfm \Omega^*$.

As shown in (ref), AdaWPCA outperforms PCA and HeteroPCA except for the FRED-QD dataset with $q_{\sf tr}=0.9$ and $r=1$, which is still relatively close. As $q_{\sf tr}$ increases, the performance of all methods improves as we have less noisy data. On the other hand, our method is much more robust to over-specified rank, as the performance of PCA and HeteroPCA becomes dramatically worse when $r>1$.

\setcounter{page}{1}