The Spurious Factor Dilemma: Robust Inference in Heavy-Tailed Elliptical Factor Models

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.

82,110 characters · 16 sections · 45 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractStandard methods for determining the number of factors often overestimate the true number when data exhibit heavy-tailed randomness, misinterpreting noise-induced outliers as genuine factors. This paper addresses this challenge within the framework of Elliptical Factor Models (EFM), which accommodate both heavy tails and potential non-linear dependencies common in real-world data. We demonstrate, both theoretically and empirically, that heavy-tailed noise generates spurious eigenvalues that mimic true factor signals. To distinguish these, we propose a novel methodology based on a fluctuation magnification algorithm. Under mild conditions, we show that, by magnifying perturbations, the eigenvalues associated with real factors exhibit significantly less fluctuation (stabilizing asymptotically) than spurious eigenvalues arising from heavy-tailed effects. We develop a formal testing procedure based on this principle and apply it to the problem of accurately selecting the number of common factors in heavy-tailed EFMs. Simulation studies and real data analysis confirm the effectiveness of our approach, particularly in scenarios with pronounced heavy-tailedness.

KEY WORDS: Elliptical distributions; Factor models; Heavy tails; Spurious factors.

\addtocontents{toc}{\setcounter{tocdepth}{2}}

Introduction

Factor models serve as a cornerstone in the analysis of large-scale datasets across various disciplines, including economics, finance, genetics, and signal processing. Their power lies in the ability to parsimoniously capture complex dependencies and interactions among numerous variables by attributing them to a small number of latent common factors bai2002determining, stock2002forecasting, fan2018large. By effectively modeling these latent structures, factor analysis provides crucial tools for dimension reduction, feature extraction, and understanding the underlying drivers of observed phenomena fan2021recent.

A fundamental challenge in applying factor models is the determination of the correct number of common factors, denoted by $m$. This problem has received considerable attention, yet remains a subject of ongoing research due to its critical impact on subsequent analysis. Underestimating $m$ leads to the omission of significant systematic components, potentially resulting in biased estimates of factor loadings, inconsistent forecasting, and flawed structural interpretations bai2003inferential, baltagi2017identification. Conversely, overestimating $m$ introduces noise by fitting spurious factors, which can inflate estimation variance, reduce model interpretability, and increase computational costs barigozzi2020consistent.

Existing methodologies for estimating $m$ largely fall into two categories. The first, prevalent in econometrics, leverages the connection between factor analysis and Principal Component Analysis (PCA). Methods like the information criteria of bai2002determining and alessi2010improved, the eigenvalue ratio tests of ahn2013eigenvalue and lam2012factor, and the randomization tests trapani2018randomized, kong2020random rely on the assumption that eigenvalues associated with common factors diverge at a faster rate than those corresponding to idiosyncratic noise as dimensions grow. The second category employs Random Matrix Theory (RMT) to provide finer distinctions, particularly in high-dimensional settings where both the number of variables ($p$) and observations ($n$) are large. RMT-based tests, such as those by onatski2009testing, utilize the fact that the largest noise eigenvalues converge to the Tracy–Widom distribution, while factor-related eigenvalues appear as distinct outliers, often asymptotically Gaussian after proper scaling onatski2010determining, cai2020limiting, ke2023estimation.

However, both classes of methods face limitations when faced with data exhibiting heavy-tailed randomness. Heavy-tailed distributions, characterized by a higher probability of extreme events compared to Gaussian distributions, are ubiquitous in financial returns, climate data, and various other fields roy2021empirical,ke2023estimation. PCA-based methods, while relatively robust to certain types of noise dependence, can falter because heavy tails can generate large sample eigenvalues purely from noise, mimicking the signature of true factors. RMT-based methods typically rely on moment conditions or concentration properties that are violated by heavy-tailed noise, leading them to misinterpret large noise eigenvalues as signals. Figure (ref) provides a visual example of this issue, where a large spurious eigenvalue appears far from the bulk, potentially misleading standard selection criteria.

Furthermore, empirical data often exhibit non-linear dependencies alongside heavy tails. Standard factor models typically assume idiosyncratic errors or linear dependencies, potentially missing complex interaction patterns. Elliptical Factor Models (EFM), based on elliptical distributions, offer a flexible framework that naturally incorporates both heavy-tailedness (via the radial component) and non-linear dependencies (via the elliptical structure) chamberlain1982arbitrage,baltagi2017identification. Recently, bao2025signal has demonstrated that even heteroscedastic or cross-correlated noise with independent entries along the time dimension can produce misleading spikes, which traditional singular value methods may incorrectly interpret as signals. As a result, spurious eigenvalues are quite common in practice.

This paper addresses the critical issue of factor number overestimation in the EFM framework, specifically focusing on the confusion caused by heavy-tailed noise. We pose two central questions: Can we reliably detect spurious factors generated by heavy tails in elliptical models? Can we develop a robust factor selection procedure for such data?

Our primary contribution is a novel methodology with a rigorous theoretical basis for distinguishing between “real" factor signals and “spurious" noise-induced signals among the large sample eigenvalues. We introduce a novel algorithm called the fluctuation magnification algorithm, specifically designed for sample covariance matrices derived from EFM data. The core idea is that perturbing the data via an elaborated magnifier affects real and spurious signals differently. We theoretically establish that, under the fluctuation magnification, the (appropriately scaled) eigenvalues corresponding to true common factors exhibit stability, converging to a normal distribution or having small variance relative to their magnitude. In contrast, spurious eigenvalues generated by heavy tails display significantly larger fluctuations under the perturbation. This difference in stability provides a clear mechanism for detection.

Based on these distinct asymptotic behaviors (detailed in Sections (ref) and (ref)), we develop a test statistic that quantifies the fluctuation of each large sample eigenvalue under the magnification algorithm. This allows us to formally test for the presence of spurious factors (Section (ref)). As a key application, we integrate this detection mechanism into a two-step procedure to robustly estimate the number of common factors ($m$) in heavy-tailed EFMs. The first step uses the magnification algorithm to identify potential spurious signals among the leading eigenvalues. The second step employs existing criteria (e.g., from onatski2010determining) as a safeguard, primarily to handle cases where the noise might be light-tailed, ensuring consistency across different tail behaviors.

Our work can be viewed as being generated from high-dimensional resampling techniques lopes2019bootstrapping, han2018gaussian, ding2023extreme, ke2023estimation, yao2021rates, yu2024testing to the challenging setting of elliptical factor models with heavy tails, providing the first procedure, to our knowledge, specifically designed to detect spurious factors in this context. Numerical simulations demonstrate the superior performance of our method compared to established techniques, especially in heavy-tailed scenarios. Application to real financial data yields results consistent with financial theory and highlights the practical relevance of addressing spurious factors.

The remainder of the paper is organized as follows. Section (ref) formally introduces the Elliptical Factor Model and outlines the key assumptions. Section (ref) further illustrates the problem of spurious factors using examples. Section (ref) presents the main asymptotic theory that details the behavior of real and spurious eigenvalues. Section (ref) describes the fluctuation magnifier algorithm and the proposed testing and factor selection procedures. Section (ref) provides simulation results, and Section (ref) discusses the real data application. We provide a sketch for our proof strategy for the theoretical results in Section (ref). The conclusion is offered in Section (ref). All detailed technical proofs are deferred to Appendix (ref).

Conventions

Let $\mathbb{C}_+$ denote the complex upper half-plane. We use $C > 0$ to represent a generic positive constant whose value may change from line to line. For sequences of positive deterministic values $\{a_n\}$ and $\{b_n\}$, $a_n = \mathrm{O}(b_n)$ means $a_n \leq C b_n$ for some $C > 0$. If $a_n = \mathrm{O}(b_n)$ and $b_n = \mathrm{O}(a_n)$, we write $a_n \asymp b_n$. We write $a_n = \mathrm{o}(b_n)$ if $a_n \leq c_n b_n$ for some positive sequence $c_n \downarrow 0$. For a sequence of random variables $\{x_n\}$ and positive real values $\{a_n\}$, $x_n = \mathrm{O}_{\mathbb{P}}(a_n)$ indicates that $x_n / a_n$ is stochastically bounded. $x_n = \mathrm{o}_{\mathbb{P}}(a_n)$ means $x_n / a_n$ converges to zero in probability. For a sequence of positive random variables $\{y_n\}$, $y_{(k)}$ denotes the $k$-th order statistic, $y_{(1)} \geq y_{(2)} \geq \cdots \geq y_{(n)} > 0$. Vectors are marked in bold.

Elliptical factor model and assumptions

Elliptical distributions provide a versatile class for modeling multivariate data, extending the normal distribution to accommodate heavy tails and capture specific dependence structures like tail dependence, making them particularly relevant in finance and other fields chamberlain1982arbitrage, fama1993common, baltagi2017identification. A $p$-dimensional random vector $\mathbf{y}$ follows a centered elliptical distribution, denoted $\mathbf{y} \sim EC_p(0, \Sigma, \xi)$, if it admits the stochastic representation:

equation[equation omitted — 98 chars of source]

where $\Sigma \in \mathbb{R}^{p \times p}$ is a positive definite matrix representing the population covariance matrix, $\xi \ge 0$ is a non-negative scalar random variable representing the “radius," and $\mathbf{u} \in \mathbb{R}^p$ is a random vector uniformly distributed on the unit sphere $\mathbb{S}^{p-1}$, independent of $\xi$. The variable $\xi$ governs the tail behavior of the distribution.

We integrate this structure with the standard linear factor model. Let $\mathbf{y}_1, \dots, \mathbf{y}_n$ be independent and identically distributed (i.i.d.) random vectors in $\mathbb{R}^p$. The factor model posits:

equation[equation omitted — 114 chars of source]

where $B \in \mathbb{R}^{p \times m}$ is the deterministic factor loading matrix, $\mathbf{f}_t \in \mathbb{R}^m$ is the vector of $m$ latent common factors, and $\mathbf{e}_t \in \mathbb{R}^p$ is the vector of idiosyncratic errors. We assume $m$ is fixed and much smaller than $p$ and $n$. Standard identification conditions often include $\frac{1}{n} \sum_{t=1}^n \mathbf{f}_t \mathbf{f}_t' \overset{\mathbb{P}}{\rightarrow} I_m$ as $p \to \infty$ fan2018large.

To define the Elliptical Factor Model (EFM), we assume that $\mathbf{f}_t$ and $\mathbf{e}_t$ are uncorrelated and the joint distribution of factors and errors follows an elliptical structure. Specifically, if $(\mathbf{f}_t', \mathbf{e}_t')'$ is elliptically distributed with mean zero and covariance matrix $\operatorname{diag}(I_m, \Sigma_{err})$, where $\operatorname{diag}(I_m, \Sigma_{err})$ means that

align*[align* omitted — 137 chars of source]

Then $(\mathbf{f}_t', \mathbf{e}_t')' \sim EC_{m+p}(\mathbf{0}, \operatorname{diag}(I_m, \Sigma_{err}), \xi_t)$. Consequently, the observation vector $\mathbf{y}_t$ also follows an elliptical distribution: \[ \mathbf{y}_t=

pmatrix[pmatrix omitted — 21 chars of source]

pmatrix[pmatrix omitted — 50 chars of source]

\sim EC_p(\mathbf{0}, \Sigma, \xi_t), \] where

align*[align* omitted — 40 chars of source]

is the population covariance matrix of $\mathbf{y}_t$ (assuming $\mathbb{E}[\xi_t^2]$ is finite and normalized appropriately). The stochastic representation for the observations becomes:

equation[equation omitted — 115 chars of source]

where $\{\xi_t\}_{t=1}^n$ are i.i.d. copies of the radius variable $\xi$, and $\{\mathbf{u}_t\}_{t=1}^n$ are i.i.d. uniform on $\mathbb{S}^{p-1}$, independent of $\xi$. The data matrix is $Y = (\mathbf{y}_1, \dots, \mathbf{y}_n)= \Sigma^{1/2} U D$, where $U = (\mathbf{u}_1, \dots, \mathbf{u}_n)$ and $D = \operatorname{diag}(\xi_1, \dots, \xi_n)$.

The key feature allowing for heavy tails is the random radius $\xi_t$. We impose the following assumption on its distribution.

assumption[Distribution of Radius] Let $D^2 = \operatorname{diag} \{ \xi_1^2, \dots, \xi_n^2 \}$. The entries $\xi_t^2 \sim \xi^2$ for $1 \leq t \leq n$ are i.i.d. copies of a non-negative, non-degenerate random variable $\xi^2$ with support on $\mathbb{R}^{+}$. We assume $\mathbb{E}[\xi^2] = 1$ without loss of generality and $\xi^2$ satisfies either of the following tail conditions: (i) Polynomial decay tail: \begin{equation} \mathbb{P}(\xi^2 > x) = x^{-\alpha} L(x), \quad as x \to \infty, \end{equation} for some $\alpha \in (1, \infty)$, where $L(x)$ is a slowly varying function, i.e., $$\lim_{x \to \infty} L(sx)/L(x) = 1$$ for all $s>0$. (ii) Exponential decay tail: For some constant $\beta > 0$ and any fixed constant $s > 0$, \begin{equation} \mathbb{E} e^{s\xi^{2\beta}} < \infty. \end{equation}

remark[Heavy Tails] Assumption (ref) explicitly allows for heavy tails. The polynomial decay case (ref) includes distributions where $\mathbf{y}_t$ has finite variance but potentially infinite higher moments (e.g., multivariate t-distributions with degrees of freedom $\nu > 2$; here $\alpha = \nu/2$). The case $\alpha \in (1, 2]$ corresponds to particularly heavy tails where the fourth moment of $\xi^2$ (and potentially $\mathbf{y}_t$) may not exist. We term this “serious heavy-tailed randomness." The case $\alpha \in (2, \infty)$ or the exponential decay case (ref) covers distributions with finite fourth moments but still heavier tails than Gaussian, termed “typical heavy-tailed randomness." This contrasts with settings assuming $\xi$ is constant or bounded hu2019high. The threshold \begin{equation*} \mathsf{T} := \begin{cases} n ^{1/\alpha}\log n, & if (ref) holds; \\ (\log n)^{1/\beta}, & if (ref) holds. \end{cases} \end{equation*} characterizes the typical maximal order of $\xi_t^2$. As shown later (Lemma (ref)), $\max_{1\le t \le n} \xi_t^2 \lesssim \mathsf{T}$ with high probability. Note the slight adjustment in definition compared to the draft for technical convenience in proofs related to RMT results under heavy tails.

We assume the population covariance matrix $\Sigma$ exhibits a spiked structure, which is common in factor analysis chamberlain1982arbitrage, yu2024testing.

assumption[Spiked Covariance Structure] The population covariance matrix $\Sigma = BB' + \Sigma_{err}$ satisfies: (1) $\max\{\|\Sigma_{err}\|, \|\Sigma_{err}^{-1}\|\}\le c_1$, for some constant $0 < c_1< \infty$. (2) $BB'$ is of rank $m$ for some fixed integer $m>0$. The descending eigenvalues of $BB'$, $(\lambda_i(BB^{\prime}))^{-1}\mathsf{T}=\mathrm{o}(1)$ for $1\le i\le m$, and $1+c_2\le\lambda_i(BB^{\prime})/\lambda_{i+1}(BB^{\prime})\le c_3 $ for any $1\le i\le m-1$ and some $0 < c_2<c_3< \infty$. (3) The dimensions satisfy $p/n \to \phi \in (c_4, c_4^{-1})$ for some constant $0 < c_4< \infty$, as $\min\{p, n\} \to \infty$.

remarkIn Assumption (ref), (1) is natural, as it imposes that the eigenvalues of the error covariance matrix are uniformly bounded from above and below. (2) indicates that $\Sigma$ has a spiked structure. It is easy to see that $\Sigma$ will inherit the key feature of $BB^{\prime}$ and $\Sigma_{err}$ that, there will be $(\Sigma_{ii})^{-1}\mathsf{T}=\mathrm{o}(1)$ for $1\le i\le m$ and $c\le\Sigma_{ii}\le c^{-1}$ for any $i\ge m+1$. In the following, we denote the ordered components in $\Sigma$ as $\sigma_1>\sigma_2>\dots>\sigma_m>\sigma_{m-1}\ge \dots\ge \sigma_p\ge c$. Also, notice that we demand the large spike of $\Sigma$ should exhibit the same order if they are divergent. $(3)$ confirms the structure of high dimensions, and the fixed $m$ is common in the literature, especially when the target is to determine the number of common factors fan2018large,yu2024testing.

remarkWe want to emphasize that although Assumption (ref) imposes the restriction $(\Sigma_{ii})^{-1}\mathsf{T} = \mathrm{o}(1)$ for $1 \le i \le m$, it is generally difficult to distinguish real signals from spurious ones based solely on their asymptotic divergence rates from the observed data, since the presence of heavy-tailed randomness in the factor noise cannot be ruled out.

remarkThe main difference between Assumption (ref) and typical settings in factor analysis is that we specify only the configuration of the common factors in $\Sigma$ without imposing constraints on the random noise. Additionally, from Assumption (ref), it is evident that we allow for a broad range of heavy-tailed randomness in factor noise.

The problem: spurious factors from heavy tails

Given observations $\{\mathbf{y}_t\}_{t=1}^n$ from the EFM (ref), we form the sample covariance matrix $S$ and its companion $\mathcal{S}$:

equation[equation omitted — 152 chars of source]

Note that $S$ and $\mathcal{S}$ share the same non-zero eigenvalues, conventionally denoted $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{\min(p,n)}$. Standard factor analysis methods aim to estimate $m$ by examining these top eigenvalues. The underlying assumption is that $\lambda_1, \dots, \lambda_m$ are primarily influenced by the population spikes $\sigma_1, \dots, \sigma_m$, while $\lambda_{m+1}, \dots$ reflect the noise structure $\Sigma_{err}$.

However, under the EFM with heavy tails (Assumption (ref)), this distinction becomes blurred. The random radii $\xi_t^2$ can take extremely large values. When a large $\xi_t^2$ aligns with certain directions in $U$ and $\Sigma^{1/2}$, it can inflate some sample eigenvalues dramatically, even those notionally corresponding to the noise part $\Sigma_{err}$. These inflated noise eigenvalues can become outliers, exceeding the bulk and potentially mixing with, or even surpassing, the eigenvalues generated by the true factors $B\mathbf{f}_t$. We term these large eigenvalues not directly associated with the $m$ population spikes as “spurious factors".

Figure (ref) provides a toy example. Standard methods relying on gaps or thresholding (e.g., ahn2013eigenvalue, onatski2010determining) would likely identify two factors ($m=2$) instead of the true $m=1$, mistaking the eigenvalue at 5.5 for a real factor signal. This overestimation poses a significant problem in practice. Therefore, the primary goal of this paper is to develop a methodology to detect these spurious factors. Let $\mathsf{f}$ denote the number of spurious factors among the top eigenvalues (i.e., the number of large eigenvalues $\lambda_k$ that do not correspond to the true population spikes $\sigma_1, \dots, \sigma_m$). We aim to test the hypothesis:

equation[equation omitted — 126 chars of source]

Rejecting $\mathbf{H}_0$ implies the presence of spurious factors among the leading eigenvalues. This detection is crucial for accurately estimating the true number of factors $m$. Our proposed approach, detailed in Section (ref), leverages the differential fluctuation of real and spurious signals under magnification perturbations.

figure[figure omitted — 526 chars of source]

Asymptotic behavior of sample eigenvalues

This section provides the theoretical foundation for our detection method by characterizing the asymptotic behavior of both real factor signals and spurious noise components present in the sample eigenvalues $\{\lambda_k\}$. Due to the spiked structure of $\Sigma$, we decompose $\Sigma=\Sigma_1+\Sigma_2$, where $\Sigma_1$ contains the $m$ largest components of $\Sigma$, and $\Sigma_2$ captures the remaining ones (also recall that $\Sigma$ is diagonal from Assumption (ref)). As a consequence, we can rewrite

equation[equation omitted — 102 chars of source]

In the sequel, we focus on the eigenvalues of $\mathcal{S} = Y'Y$. Note that $S$ and $\mathcal{S}$ share the same non-zero eigenvalues.

Asymptotic behavior of real factors

The top $m$ eigenvalues $\lambda_1, \dots, \lambda_m$ are primarily influenced by the population spikes $\sigma_1, \dots, \sigma_m$. However, their exact location is perturbed by the noise component $\Sigma_{err}$ and the heavy-tailed radii $D^2$.

Following standard techniques for spiked models, we first establish a first-order approximation.

theorem[First Order Approximation] Under Assumptions (ref) and (ref), we have: \begin{itemize} • for polynomial decay tail (ref) with $\alpha\in(2,+\infty)$ or exponential decay tail (ref) in Assumption (ref), \begin{equation*} \frac{\lambda_i}{\sigma_i} = 1 + \mathrm{O}_{\mathbb{P}}\left(\frac{\mathsf{T}}{\sigma_i} + \frac{1}{\sqrt{n}}\right), 1 \le i \le m; \end{equation*} • for polynomial decay tail (ref) with $\alpha\in(1,2]$ in Assumption (ref), \begin{equation*} \frac{\lambda_i}{\sigma_i} = 1 + \mathrm{O}_{\mathbb{P}}\left(\frac{\mathsf{T}}{\sigma_i} + \sqrt{\frac{\mathsf{T}}{n}}\right), 1 \le i \le m. \end{equation*} \end{itemize}

remarkTheorem (ref) indicates that the leading eigenvalues of $S$ are very close to the population ones in the sense that their ratios asymptotically converge to one in probability. However, this doesn't mean that $\lambda_i$ is a consistent estimator of $\sigma_i$ since $\sigma_i$ tends to infinity. On the other hand, the convergence rate in Theorem (ref) can be slower than the order of $n^{-1/2}$ due to the term $\mathsf{T}/\sigma_i$ in case $(1)$. Therefore, the results in Theorem (ref) are far from ideal, especially if our goal is to determine the asymptotic distribution of $\lambda_i$.

To obtain a more precise description for the locations of $\lambda_i$, we employ the methodology from yu2024testing and introduce the random quantities $\theta_i,1\le i\le m$ to trace the asymptotic behavior of each $\lambda_i$. We define $\theta_i$ as the unique solution to

equation[equation omitted — 231 chars of source]

DXYZspiked showed that under the elliptical model, $\theta_i$ is a closer approximation to $\lambda_i$ compared with $\sigma_i$. In our study, to address the effect of $D^2$, we need the second order description $\zeta_i$ to be the unique solution to

equation[equation omitted — 269 chars of source]

The existence and uniqueness of $\theta_i, \zeta_i$ are guaranteed under our assumptions for large $p, n$.

lemma[Properties of $\theta_i, \zeta_i$] Under Assumptions (ref) and (ref), the solutions $\theta_i$ and $\zeta_i$ exist with probability tending to one as $n\rightarrow\infty$. It holds that \begin{equation*} \frac{\theta_i}{\sigma_i}=1+\mathrm{O}(\frac{\operatorname{tr}\Sigma_2}{n\sigma_i}), \end{equation*} and for polynomial decay tail (ref) with $\alpha\in(2,+\infty)$ or exponential decay tail (ref) in Assumption (ref), \begin{equation*} \zeta_i-\frac{\theta_i}{\sigma_i}=\frac{1}{p}\sum_{j=1}^n(\xi^2_j-\phi)+\frac{\operatorname{tr}\Sigma_2}{p\sigma_i\theta_i}\times\mathbb{E}\big[(\xi^2_1-\phi)^2\big]+\mathrm{o}_{\mathbb{P}}\Big(\frac{1}{\sqrt{n}}\Big); \end{equation*} while for polynomial decay tail (ref) with $\alpha\in(1,2]$, \begin{equation*} \zeta_i-\frac{\theta_i}{\sigma_i}=\frac{1}{p}\sum_{j=1}^n(\xi^2_j-\phi)+\frac{\operatorname{tr}\Sigma_2}{p\sigma_i\theta_i}\times\mathbb{E}\big[(\xi^2_1-\phi)^2\big]+\mathrm{o}_{\mathbb{P}}\Big(\sqrt{\frac{\mathsf{T}}{n}}\Big). \end{equation*}

Using these observations, we obtain a more precise characterization of $\lambda_i$. Let $\mathbf{e}_k$ be the $k$-th standard basis vector.

theorem[Second Order Approximation] Under Assumptions (ref) and (ref), we have that: \begin{itemize} • for for polynomial decay tail (ref) with $\alpha\in(2,+\infty)$ or exponential decay tail (ref) in Assumption (ref), \begin{equation*} \frac{\lambda_i}{\theta_i}-1=\mathbf{e}^{\prime}_iUD^2U^{\prime}\mathbf{e}_i-\frac{1}{p}\operatorname{tr}D^2-\frac{\theta_i}{\sigma_i}+\zeta_i+\mathrm{O}_{\mathbb{P}}\big(\frac{\mathsf{T}}{\sqrt{n}\sigma_i}+\frac{1}{n}\big), 1 \le i \le m; \end{equation*} • for polynomial decay tail (ref) with $\alpha\in(1,2]$, \begin{equation*} \frac{\lambda_i}{\theta_i}-1=\mathbf{e}^{\prime}_iUD^2U^{\prime}\mathbf{e}_i-\frac{1}{p}\operatorname{tr}D^2-\frac{\theta_i}{\sigma_i}+\zeta_i+\mathrm{O}_{\mathbb{P}}\big(\frac{\mathsf{T}}{\sqrt{n}\sigma_i}+\frac{\mathsf{T}}{n}\big), 1 \le i \le m. \end{equation*} \end{itemize}

Then, related to the results in Lemma (ref), Theorem (ref) essentially implies the improved asymptotically rate $\lambda_i/\theta_i-1=\mathrm{O}_{\mathbb{P}}\big(\mathsf{T}/\sigma_i^2+(\sqrt{n}\sigma_i)^{-1}+\frac{1}{\sqrt{n}}+\sqrt{\frac{\mathsf{T}}{n}}\times\mathbbm{1}(\alpha\in(1,2])\big)$, which can be seen as an improvement of Theorem (ref) especially for $\alpha\in(2,+\infty)$ and exponential decay tail. On the other hand, we notice that the error rate in Theorem (ref) will get large when the randomness of $\mathbf{y}$ exhibits heavier distributions (in case of the index $\alpha$). Then, it suggests that the extreme heavy-tailness will strongly influence the convergence rates of the spiked eigenvalues.

Specifically, in the case where $\xi^2$ either has a polynomial decay tail with $\alpha\in(2,+\infty)$ or has an exponential decay tail in Assumption (ref), the following central limit theorem holds for real signals from Theorem (ref),

theoremAssume (ref) with $\alpha\in(2,+\infty)$ or (ref) in Assumption (ref), also assume Assumption (ref). Then we have that for $1\le i\le m$ \begin{equation} \sqrt{n}\Big({\lambda_i}/{\theta_i}-(1+\mathrm{o}(1))\Big)\overset{d}{\rightarrow}\mathrm{N}(0,3\mathbb{E}\xi_1^4-1). \end{equation}

remarkIt should be noted that $\theta_i$ is difficult to compute in practice; however, the above theorem establishes the asymptotic normality of the real factors, which serves as a crucial input for our testing procedure.

When $\xi^2$ exhibits serious heavy-tailed decay with $\alpha\in(1,2]$, the asymptotic normality in Theorem (ref) will not hold, but demands a larger scaling. In practice, it is enough to only consider the asymptotic mean and variance for $\lambda_i/\theta_i$. We have the following results for real factors under the serious heavy-tailed scenario.

theoremSuppose (ref) holds with $\alpha\in(1,2]$ in Assumption (ref). Then we have that for $1\le i\le m$ \begin{equation} \mathbb{E}({\lambda_i}/{\theta_i})=1+\mathrm{o}(1),\quad \operatorname{Var}({\lambda_i}/{\theta_i})=\mathrm{O}({3\mathsf{T}}/{n}). \end{equation}

Asymptotic behavior of spurious factors

We first prepare several notations. Recalling the decomposition of $\mathcal{S}$ in (ref), we denote the matching matrices removing the spike structure from $\mathcal{S}$ (or $S$) as

equation[equation omitted — 128 chars of source]

Denote $\lambda_{i,2}$ as the $i$-th largest eigenvalue of $S_2$ or $\mathcal{S}_2$. Denote the quantity $\mathsf{q}$ as

equation[equation omitted — 272 chars of source]

for some sufficiently small constant $\epsilon>0$.

The following lemma indicates that the leading eigenvalues of $S$, despite the first $m$ spiked eigenvalues, are close to the leading eigenvalues of $\mathcal{S}_2$,

lemmaUnder Assumptions (ref) and (ref), we have that for $k\ge1$, \begin{equation} |\lambda_{m+k}-\lambda_{k,2}|=\mathrm{O}_{\mathbb{P}}(n^{-1/2+2\epsilon}\mathsf{q}). \end{equation}

Typically, the largest eigenvalues of $\mathcal{S}_2$ are strongly influenced by the ordered statistics in $D^2$, whose phenomena is captured in the following theorem,

theoremUnder Assumptions (ref) and (ref), we have that for $k\ge1$ \begin{equation} {\lambda_{k,2}}/{\xi_{(k)}^2}=\Bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1), \end{equation} where $\Bar{\sigma}:=p^{-1}\operatorname{tr}(\Sigma_2)$ and $\xi_{(1)}^2\ge \xi_{(2)}^2\ge\cdots\ge \xi_{(n)}^2$ are the ordered statistics from $D^2$.

As a result of Theorem (ref), we may easily find that $\mathbb{E}(\lambda_{k,2})=\xi_{(k)}^2(\Bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1))$, given any realization of $D^2$.

remarkFrom Lemma (ref), we observe that the spurious noise signals closely approximate the leading eigenvalues of $\mathcal{S}_2$. Furthermore, Theorem (ref) and Lemma (ref) reveal that these leading eigenvalues of $\mathcal{S}_2$ diverge in alignment with the ordered statistics in $D^2$. This theoretical insight explains why true common factors and spurious noise signals become indistinguishable among the top eigenvalues of the data matrices.

Detection via Fluctuation Magnification

In this section, we develop testing strategies for the problem (ref), guided by the theoretical framework established in Sections (ref) and (ref). We begin by rigorously formalizing the fluctuation magnification mechanism from the idea of resampling techniques for $S$, followed by a systematic algorithm to implement the procedure to identify spurious factors.

Fluctuation magnification for $S$

Suppose we observe the data matrix $Y$. We define a sequential magnifier data matrices $\widetilde{Y}_j,1\le j\le K$ from $Y$ for some large fixed integer $K>0$. Let

equation*[equation* omitted — 122 chars of source]

where $w_t\sim w, 1\le t\le n$ are some properly chosen i.i.d. random variables whose details will be specified below. For $1\le j\le K$, we independently construct from $\widetilde{Y}_j$ with

gather*[gather* omitted — 311 chars of source]

Let $\lambda_k^{(j)}$ be the $k$-th largest eigenvalue of $\widetilde{S}^{(j)}$ (or $\widetilde{S}^{(j)}$) and $\lambda_{k,2}^{(j)}$ the $k$-th largest eigenvalue of $\widetilde{\mathcal{S}}_2^{(j)}$. In the sequel, we use the notation $\mathbb{P}^*$ to denote the probability measure conditional on the sample $Y$ by $\mathbb{P}^*(\cdot)=\mathbb{P}(\cdot|\mathcal{F}_Y)$ where $\mathcal{F}_Y=\sigma(Y)$ is the sigma algebra generated from $\{\mathbf{y}_1,\dots,\mathbf{y}_n\}$. Then, we define $\overset{d^*}{\rightarrow}$, $\mathrm{o}_{\mathbb{P}^*}$ and $\mathrm{O}_{\mathbb{P}^*}$ accordingly from $\mathbb{P}^*$.

We first concern ourselves with the real factors. We consider the two situations in Theorem (ref) and Theorem (ref). Applying these two theorems to each $\widetilde{S}^{(j)}$ for $1\le j\le K$, we have the following proposition.

propositionFor all $1\le j\le K$ and $1\le i\le m$, \begin{itemize} • if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(2,+\infty)$ or exponential decay tail as in Theorem (ref), then we have that \begin{equation*} \sqrt{n}\big(\frac{\lambda_i^{(j)}}{\theta_i}-(1+\mathrm{o}(1))\big)\overset{d^*}{\rightarrow}\mathrm{N}(0,3\mathbb{E}[\xi_1^4w_1^2]-1); \end{equation*} • if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(1,2]$ as in Theorem (ref), then we have that \begin{equation*} \mathbb{E}({\lambda_i}/{\theta_i})=1+\mathrm{o}(1),\quad \operatorname{Var}({\lambda_i}/{\theta_i})=\mathrm{O}({3\mathsf{T}}/{n}), \end{equation*} \end{itemize} given in both cases that $\sigma_m\gg\mathsf{T}w_{(1)}$ where $w_{(1)}$ is the largest order statistic of $\{w_t\}_{1\le t\le n}$.

As a consequence, we have the following proposition.

propositionUnder the assumptions of Proposition (ref), for $1\le i\le m$ and sufficiently large $K$, \begin{itemize} • if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(2,+\infty)$ or exponential decay tail as in Theorem (ref), then we have that \begin{equation*} \frac{1}{K}\sum_{j=1}^K\Big({K\lambda_i^{(j)}}/{\sum_{s=1}^K\lambda_i^{(s)}}-1\Big)^2=\mathrm{O}_{\mathbb{P}^*}(1/n); \end{equation*} • if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(1,2]$ as in Theorem (ref), then we have that \begin{equation*} \frac{1}{K}\sum_{j=1}^K\Big({K\lambda_i^{(j)}}/{\sum_{s=1}^K\lambda_i^{(s)}}-1\Big)^2=\mathrm{O}_{\mathbb{P}^*}(\mathsf{T}/n). \end{equation*} \end{itemize}

Next, we turn to the spurious factors perturbed by some magnifier matrices. Lemma (ref) and Theorem (ref) indicate that for $m+1\le i\le \hat{o}$,

equation*[equation* omitted — 269 chars of source]

To simplify the discussion, we restrict ourselves to $i=m+1$ while the other cases can be handled similarly. The key idea here is to choose a suitable $w$ such that the fluctuation of $\lambda^{(j)}_{m+1}$ will not degenerate with $n$. A good candidate is that $w$ is uniformly distributed on the interval $[a,b]$ (say $w\in U[a,b]$) with $(a+b)/2=1$ and $a\ge b\xi^2_{(2)}/\xi^2_{(1)}$, which satisfies the condition in Proposition (ref). Then, it is easy to find that $(\xi^2w)_{(1)}$ is uniformly distributed on $[a\xi^2_{(1)},b\xi^2_{(1)}]$ with

gather*[gather* omitted — 354 chars of source]

Since $\lambda_{m+1}^{(j)}=(\Bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1))\cdot(\xi^2w)_{(1)}$ are i.i.d. across each fluctuation magnification procedure, we can apply CLT to obtain the following results.

propositionUnder the assumptions in Theorem (ref) with the choice of $w\in U[a,b]$, it holds for sufficiently large $K$ that \begin{equation} \frac{1}{K}\sum_{j=1}^K\Big({K\lambda_{m+1}^{(i)}}/{\sum_{s=1}^K\lambda_{m+1}^{(s)}}-1\Big)^2=c+\mathrm{O}_{\mathbb{P}^*}(\mathsf{T}^2K^{-1/2}). \end{equation} Here $c>0$ is a constant only depending on the choice of $w$.

remarkSeveral remarks are in order. First, it is important to choose suitable magnifiers $w$ such that $c$ will not degenerate and therefore there will be a clear distinction with the error terms in {Proposition (ref).} Second, parallel results can also be obtained for $\lambda_{i}^{(j)}, m+2\le i\le o$ where $o$ is some pre-given value that counts the whole number of outliers of $S$. However, in real practice, there is no need to consider all the outliers, since the first spurious signal is enough to establish the borderline between real factors and spurious factors. Third, the distinct behavior of the real factors in {Proposition (ref)} and spurious factors in {Proposition (ref)} gives the opportunity to detect the spurious factors through the fluctuation magnification mechanism, which is the main content in the next section. Fourth, we emphasize that $K$ should be sufficiently large to reduce the error level, indicating that the magnification procedures should be repeated enough times. Finally, the detailed proof of the results in this section will be provided in the Appendix.

Detection algorithm

Inspired by the theoretical observations in Section (ref), we propose the following statistics.

equation[equation omitted — 157 chars of source]

Recap the test problem (ref) in Section (ref). We propose a novel algorithm that leverages a fluctuation magnification mechanism to detect potentially spurious signals and thereby determine the true number of common factors. To achieve this, we conduct a two-step testing procedure. In the first step, we apply the fluctuation magnification mechanism to identify spurious signals from the sample matrix $S$ without distinguishing outliers in advance. As a result, some bulk eigenvalues (non-outliers) may also be flagged as spurious. In the second step, we verify whether any bulk components were mistakenly identified as spurious by employing standard methods such as onatski2010determining, thus refining the selection of common factors.

1. First step detection. We now proceed to the first step. Under Assumption (ref), it suffices to focus on the largest spurious signal under the alternative hypothesis $\mathbf{H}_a$. As noted in Remark (ref), it is crucial to choose an appropriate distribution for $w$ to ensure the validity of {Proposition (ref)}. Typically, $w\in U[a,b]$ satisfies the following conditions. $$\mbox{(1)}~(a+b)/2=1;~~\mbox{(2)}~ a\ge b\xi^2_{(2)}/\xi^2_{(1)};~~\mbox{(3)}~ \sigma_m\gg\mathsf{T}b.$$

However, in practice, $\xi^2_{(1)},\xi^2_{(2)}$ are unobservable, and $\sigma_m$ is difficult to determine. To solve this issue, we approximate the ratio $\xi^2_{(2)}/\xi^2_{(1)}$ using $\lambda_{m+2}/\lambda_{m+1}$, and $\sigma_m/\mathsf{T}$ using $\lambda_{m}/\lambda_{m+1}$. The validity of the first approximation is ensured by Theorem (ref), while Lemma (ref) and Theorem (ref) support the second. The algorithm is summarized in Algorithm (ref), which gives preliminary estimations for the number of spurious factors and the location of the largest potential spurious signal.

algorithm[algorithm omitted — 1,872 chars of source]

remarkSeveral remarks are in order. First, the output $r^*$ in Algorithm (ref) \[ r^*=\min\{\{1\le i\le o:\mathbb{T}_i>\mathsf{L}_i\}\cup\{o+1\}\} \] marks the first location at which a potential spurious signal is detected. And if none is detected, $r^*=o+1$ indicates that no spurious signal appears among the first $o$ indices. Second, based on the theoretical analysis in Section (ref), it holds that \begin{gather} \lim_{n,K\rightarrow\infty}\mathbb{P}(\mathbb{T}_i\le\mathsf{L}_i)=1,\;under $\mathbf{H}_0$, for $1\le i\le m$;\\ \lim_{n,K\rightarrow\infty}\mathbb{P}(\mathbb{T}_i>\mathsf{L}_i)=1,\;under $\mathbf{H}_a$, for $m+1\le i\le o$. \end{gather} In practice, we set $K=n^3$ in accordance with Proposition (ref). Third, the magnifier $w$ can be extended to other distributions, such as a Bernoulli-type random variable satisfying \begin{equation} w= \begin{cases} a & with probability p, \\ b & with probability 1-p, \end{cases} \end{equation} under conditions analogous to those specified in Algorithm (ref). As demonstrated in Step 1, for any choice of magnifier, it is essential to select appropriate regions $[a_i, b_i]$ to ensure that the ordering of ${\lambda_i}$ remains invariant under fluctuation magnification. Fourth, under Assumptions (ref) and (ref), the condition $\mathbb{T}_i\ge \mathsf{L}_i$ for all $m+1\le i\le o$ implies that heavy-tailed noise facilitates the emergence of spurious signals. However, the convolution of the stochastic components $U$ and $D$ may obscure the boundary between outliers (spurious signals) and bulk components (non-signals), potentially leading to the misclassification of the latter. To achieve a more refined estimation of the number of common factors, we introduce an auxiliary algorithm that serves as the centerpiece of the second-step testing procedure. Fifth, in numerical experiments, we adopt a more stringent threshold $\mathsf{L}_i$, particularly in the presence of heavy tails. This conservative choice prioritizes the protection of the null hypothesis, reflecting the principle that overestimation is generally more tolerable than underestimation. Finally, we anticipate that the proposed fluctuation magnification strategy is applicable to spurious signal detection in more general settings, as long as the real signals are significantly large.

2. Second step detection. First step detection aims to identify the most prominent spurious signals. The locations of the detected spurious signals provide useful information for separating potential real signals from noise contamination. When the data do not contain heavy-tailed noise, the first step procedure may incorrectly classify certain bulk eigenvalues as spurious. To mitigate this issue, we subsequently apply a refined factor‑number detection method exclusively to the corresponding potential real signals, which yields more accurate estimates and avoids over‑estimation. In this paper, we adopt the detection procedures of onatski2010determining and Dobriban as representative examples to obtain a more accurate estimate of the number of common factors in elliptical factor models that may exhibit heavy-tailed randomness. Algorithm (ref) proposed below builds on the first-step results together with the robustness properties of established factor selection methods, thereby effectively addressing the challenges posed by heavy-tailed distributions.

algorithm[algorithm omitted — 949 chars of source]

lemmaLet $\widehat{r}$ be obtained by Algorithms (ref) and (ref). Under the assumptions in Propositions (ref) and (ref), we have that \begin{equation} \lim_{n\rightarrow\infty}\mathbb{P}^*(\widehat{r}=m)=1. \end{equation}

remarkSeveral remarks are in order. First, Algorithm (ref) incorporates the methods developed in onatski2010determining and Dobriban as auxiliary inputs. This serves as a correction to our initial detection procedure when the first potential spurious signal is identified at some position $r^*$. In such cases, applying the methods from onatski2010determining or Dobriban to the original sample yields an initial estimate $k$. The final estimate is then obtained by setting $\widehat{r}=k$ if $k \le r^*-1$ (indicating that no spurious signal occurs before or at $k$), and $\widehat{r}=r^*-1$ otherwise (where the first spurious signal at $r^*$ is excluded). The consistency of this overall two-step methodology is established in Lemma (ref). Second, the choice of the specific methods in onatski2010determining and Dobriban for Algorithm (ref) is primarily due to their simplicity and established theoretical properties. They can readily be replaced by other well-known consistent procedures in the literature, such as those in bai2018consistency.

Simulation

In this section, we conduct simulation studies to validate two main aspects: (1) The overestimation of the number of common factors in EFM when noise exhibits heavy-tailed randomness, confirming the existence of spurious signals. (2) The effectiveness of our fluctuation magnification approach (referred to as the detection algorithm in Section 5.2) in distinguishing spurious signals and correctly identifying the true number of common factors in such cases. Consequently, two scenarios of heavy-tailed data $\mathbf{y}$ from the model (ref) are analyzed: (I) A typical heavy-tailed case, where $\mathbf{y}$ follows a distribution with a polynomial decay tail satisfying $\alpha \in (2, +\infty)$. (II) A serious heavy-tailed case, where $\mathbf{y}$ follows a distribution with a polynomial decay tail satisfying $\alpha \in (1, 2]$. We focus on the strength of the heavy tail that arises from polynomial decay tails rather than exponential decay tails, as the former allows easier adjustment of the strength of the heavy tail via the parameter $\alpha$.

Comparison of methodologies in common factor selection: Typical heavy-tailed case

We first establish the settings for the typical heavy-tailed case for $\mathbf{y}$ in (ref), characterized by a polynomial decay tail with $\alpha \in (2, +\infty)$. To represent various scenarios of heavy-tailed data, we employ multivariate t-distributions with degrees of freedom set at 4.3, 4.8, and 5.3 (i.e., $t(4.3)$, $t(4.8)$, and $t(5.3)$, respectively), all of which satisfy the condition $\alpha > 2$. We consider a high-dimensional framework with $p=1000$ (number of variables) and $n=1000$ (number of observations). Our analysis focuses on two population covariance matrices $\Sigma$:

itemize• $\Sigma_{\mathrm{I}} = \operatorname{diag}\{16, 8, 1, \ldots, 1\}$ • $\Sigma_{\mathrm{II}} = \operatorname{diag}\{24, 16, 8, 1, \ldots, 1\}$.

Here, the larger diagonal components represent the distinct common factors. It is clear that $\Sigma_{I}$ encompasses two common factors induced by $\{16, 8\}$, while $\Sigma_{II}$ encompasses three common factors induced by $\{24, 16, 8\}$.

In our comparative study, we consider two baseline factor-number detection methods: the ED estimator of onatski2010determining (denoted as “Onta") and the DDPA+ procedure of Dobriban (denoted as “DDPA+"). Their versions enhanced by the Fluctuation Magnification technique proposed in Section (ref) are referred to as “Onta with MF” and “DDPA+ with MF”, respectively. In our simulations, to reduce computational load, we set the magnifiers $w_{i} \sim U[0.1, 1.9]$ in Algorithm (ref) and $K=1000$ in Algorithm (ref), which yielded sufficiently good results. Detailed experimental results are presented in Figure (ref), and Tables (ref) and (ref).

Specifically, Figure (ref) illustrates the performance of the methods for a representative case (Case (i)) using the “Onta" estimator under $t(4.3)$ as an example. Panel (a) displays the first 50 eigenvalues of the sample covariance matrix. Panel (b) shows the estimated factor number $k$ obtained from the standard “Onta" method. Panel (c) demonstrates the identification of spurious signals via our proposed statistics $\mathbb{T}_i$ from Algorithm (ref); the statistics corresponding to spurious signals are markedly larger than those associated with real signals. These observations highlight the effectiveness of Algorithm (ref) in detecting spurious signals.

figure*[figure* omitted — 784 chars of source]

To compare the performance of Onta and DDPA+ with and without MF technique under Case (i) and Case (ii), we conducted 500 repeated experiments for each setting. The effective sample size (ESS) reported in Tables (ref) and (ref) corresponds to the number of valid replications retained after filtering out cases where the eigenvalue ratio between the last real and first spurious signal is too close. Specifically, for Case (i) we exclude replications with $\lambda_2/\lambda_3 < 1.1$ to avoid the second and third eigenvalues being too close; for Case (ii) we exclude those with $\lambda_3/\lambda_4 < 1.1$ to avoid the third and fourth eigenvalues being too close. This filtering ensures a clear distinction between spurious and real signals, preventing ambiguity that could otherwise distort the comparison. The tables present the underestimation (Under), overestimation (Over), and correct estimation (Right) rates for three $t$-distributions with varying degrees of freedom.

table[table omitted — 1,543 chars of source]

From Figure (ref), and Tables (ref) and (ref), we can draw the following conclusions.

itemize• Spurious signals do exist, as (a) exhibits more spiked eigenvalues compared to $\Sigma$. In such cases, traditional first-order estimation tends to overestimate when spurious signals are present. • When spurious signals emerge, the fluctuation magnification algorithm (Algorithm (ref)) can accurately identify them. This is evidenced by the significant increase in magnified variance in (c) at the locations corresponding to the first spurious signal. • As the degrees of freedom increase (i.e., the tails become lighter), the detection performance of all methods improves. However, the approaches incorporating our MF tool consistently achieve higher accuracy than their original counterparts, with the most substantial gains observed at lower degrees of freedom. • It is noteworthy that a higher underestimation rate occurring under low degrees of freedom is, paradoxically, a desirable outcome. From Tables (ref) and (ref), we observe that when the degrees of freedom are low, the effective sample size falls short of 500, indicating that the last real signal and the first spurious signal are very close. Consequently, cases where the first spurious signal even surpasses the last real signal naturally arise. In such scenarios, the underestimation precisely underscores the sensitivity of our method in detecting spurious signals.

Comparison of methodologies in common factor selection: Serious heavy-tailed case

Next, we evaluate the detection behavior under serious heavy-tailed settings for $\mathbf{y}$ described in (ref), characterized by polynomially decayed tails with $\alpha \in (1, 2]$. The remaining settings align with those in the Typical Heavy-tailed Case; only the differences are described here. We use multivariate t-distributions with degrees of freedom 2.5, 3.0, and 3.5 (denoted as $t(2.5)$, $t(3.0)$, and $t(3.5)$, respectively) to model the polynomial decay tail with $\alpha \in (1, 2]$. Additionally, the covariance matrix $\Sigma$ is configured with larger spiked eigenvalues as follows:

itemize• $\Sigma_{\mathrm{III}} = \operatorname{diag}\{240, 120, 1, \ldots, 1\}$ • $\Sigma_{\mathrm{IV}}= \operatorname{diag}\{360, 240, 120, 1, \ldots, 1\}$

Figure (ref) illustrates the performance of the methods for a representative case (Case (iii)) using the “Onta" estimator under $t(2.5)$ as an example.

figure*[figure* omitted — 792 chars of source]

Using the same evaluation framework as in the typical heavy-tailed case, we assess the underestimation, overestimation, and correct estimation rates for three $t$-distributions with varying degrees of freedom.

table[table omitted — 1,553 chars of source]

table[table omitted — 1,551 chars of source]

From Figure (ref) and Tables (ref) and (ref), we observe that the standalone DDPA+ method essentially fails under the Serious Heavy-tailed Case. Integrating our MF tool successfully restores its performance, yielding relatively accurate estimations. Moreover, the improvement achieved by employing the MF tools is even more pronounced in this case, highlighting their enhanced effectiveness under more challenging heavy-tailed settings.

Real data analysis

In this section, the proposed methods are applied to the real-world FRED-MD dataset, which was previously studied in XIA2017235 and YU2019104543. This dataset, introduced in McCracken01102016, is publicly available on the St. Louis Fed's website: https://www.stlouisfed.org/research/economists/mccracken/fred-databases. It consists of monthly data for 128 macroeconomic variables. For this analysis, we focused on the 786 observations covering the period from March 1959 to June 2024. The raw data is non-stationary and contains missing values. The first step is to transform the series into a stationary one using the code provided on the aforementioned website. After this preprocessing, the first two observations were eliminated due to the application of differencing operators, resulting in a $784 \times 128$ panel. The website also offers code to replace outliers with “reasonable" values, but we chose to omit this step, as extreme observations are inherent in data drawn from heavy-tailed distributions. Missing values were imputed with the sample mean of the corresponding non-missing entries in the same column.

We first analyzed the entire panel to determine the number of common factors. Our detection algorithm estimates $\widehat{r} = 3$, consistent with the result in XIA2017235 and slightly smaller than the estimate $\widehat{r} = 4$ proposed by YU2019104543. We consider $\widehat{r} = 3$ more reasonable than $\widehat{r} = 4$, though both studies retain outliers and recognize heavy-tailed distributions (as noted by YU2019104543). This conclusion stems from temporal dynamics analysis: rolling 4-year windows (Feb 1992-Feb 2024) reveal that spurious factors intermittently emerge in certain subperiods. These spurious factors cause an estimate of $\widehat{r} = 4$ to potentially overestimate the true factor number. Thus $\widehat{r} = 3$ provides a more robust full-sample representation.

The choice to begin this time-varying analysis in February 1992 is necessitated by data availability: Prior to this date, datasets for ACOGNO, AMDMNOx, ANDENOx, and AMDMUOx were completely missing. This extensive data gap before 1992 presents a significant challenge for our 48-month window analysis: Filling such large amounts of missing data within these relatively small observation windows would inevitably introduce substantial inaccuracies. Table (ref) shows $\mathbb{G}_i$ (“On" estimation), $i=1,\ldots,10$ of every 48 months data matrices.

table[table omitted — 868 chars of source]

We set $\mathbb{G}_i>9$ as the threshold for testing the factor, and derive that the periods “1992-1996", “1996-2000", and “2012-2016" possess two common factors. However, the period “2000-2004" may possess either two or five factors, “2004-2008" may possess either two or six factors, “2008-2012" may possess either two or nine factors, “2016-2020" may possess either two or six factors, and “2020-2024" may possess either two or nine factors. Now we conduct a fluctuation magnification algorithm on these uncertain periods. Figure (ref) illustrates the variance after fluctuation magnification for the period 2000-2004. The repetition in the magnification algorithm for variance calculation is set to $K=200$.

From Figure (ref), we observe no significant fluctuations (increases) in variance at $i=6$, indicating that there is no substantial change in the factors (i.e., no shift from spiked to bulk). In comparison, the variance shows a slight increase from $i=2$ to $i=3$. This suggests that the variance remains relatively stable (the variance remains close to 0.02), making the judgment of $i=2$ more reasonable. This is not unexpected, as fluctuation magnification for variance primarily serves as an auxiliary tool to confirm the choice of $i=2$, and the observed stability further supports this decision. However, when we examine the period “2004-2008" in Figure (ref), we observe an unusual pattern that diverges from the previous trends.

figure*[figure* omitted — 702 chars of source]

From Figure (ref), we observe significant fluctuations (from 0.02 to 0.28) in variance at $i=6$, while the fluctuations at other values remain relatively small. This suggests that $i=6$ may represent a spurious signal (noting that spikes larger than the first spurious signal are real signals), and the true number of factors is likely five. This outcome is not surprising, as the large-scale financial crisis of 2008 led to substantial disruptions in the global economy. Such economic shocks can result in changes in the underlying structure of the data, influencing a number of factors. The following Figure (ref) shows the change in the number of factors from 1992 to 2024.

Sketch of proof strategy

In this section, we outline the main ideas of the proof and full technical details are deferred to the Appendix. Our proof proceeds in three steps.

itemize• We first show that, in an elliptical factor model, heavy-tailed randomness in the noise component can generate additional outlying eigenvalues beyond those induced by divergent factor signals. This phenomenon is established in Lemma (ref) and Theorem (ref). The starting point is that the sample covariance matrix $S$ admits a separable structure under the elliptical factor model. When the radii in $D$ exhibit heavy-tailed decay, the associated eigenvalues may diverge. Our analysis relies on a perturbative argument; however, it differs substantially from the classical treatment in spiked covariance models. In the present setting, the spectral distribution is unbounded and the eigen-gaps are not necessarily large, particularly when the heavy-tailed noise is moderate. Consequently, a more delicate local-scale analysis is required. To this end, we employ tools from random matrix theory to derive a self-consistent system and establish local laws for the sample covariance matrix $S_2$ (the component without divergent factors), as developed in Theorems (ref) and (ref). On the other hand, the elliptical structure introduces nonlinear dependence across columns, which necessitates new concentration inequalities beyond the standard i.i.d. framework. Based on these preparations, we modify the perturbation arguments by isolating $\mathbf{y}_{i,2}$ corresponding to the largest noise radius $\xi_{(1)}^2$ from the matrix $Y_2$ as in (ref). We introduce a real auxiliary quantity $\mu_1>0$ to be the largest solution of $1+(\xi^2_{(1)}+\mathsf{q})m_{1n}(\mu_1)=0$ with $\mathsf{q}/\xi^2_{(1)}=\mathrm{o}_{\mathbb{P}}(1)$. As can be seen in the proof of Theorem (ref), $\mu_1$ connects with $\xi^2_{(1)}$ naturally by $\mu_1/(\xi^2_{(1)}+\mathsf{q})=\bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1)$. Moreover, it can establish a connection between $\mu_1$ and $\lambda_{1,2}:=\lambda_1(S_2)$ using a modified perturbation argument as in Theorem (ref). Specifically, as we define a determinant function $M(\cdot)$ as in (ref). $\lambda_{1,2}$ can be uniquely characterized by the equation $M(\lambda_{1,2})=0$. By our established local laws, we can further demonstrate that $1+(\xi^2_{(1)}+\mathsf{q})m_{1n}(\lambda_{1,2})\approx0$. Subsequently, a detailed continuity and stability analysis finally shows that $\lambda_{2,1}/\xi^2_{(1)}=\bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1)$. Thus, the leading eigenvalues of $S_2$ inherit the divergence rate of the largest noise radii, confirming that heavy-tailed noise generates spurious outliers. • We next analyze the leading eigenvalues of $S$ associated with divergent factor signals. Our approach builds on the perturbative framework developed in cai2020limiting, suitably generalized to accommodate elliptical noise with heavy tails. The key observation is that if the signal strengths dominate the heavy-tailed noise, the perturbative argument remains valid. A central role is played by the random quantity $\zeta_1$ defined in (ref), which captures the randomness of $\lambda_1/\theta_1$ induced by the noise radii ${\xi_i^2}$. Then, if the noise radii have polynomial decay with index $\alpha \ge 2$ or exponential decay, asymptotic normality follows from a classical central limit theorem applied to $\zeta_1$. If $\alpha \in (1,2)$, asymptotic normality fails, but the limiting mean and variance can still be characterized. Furthermore, we show that, beyond the first $m$ spikes, the leading eigenvalues of $S$ are asymptotically close to those of $S_2$, which are driven by heavy-tailed noise. This yields a complete description of the limiting behavior of eigenvalues corresponding to both true and spurious signals. • In the final step, we introduce a fluctuation magnification mechanism to distinguish genuine factors from noise-induced outliers. We perturb the data matrix $Y$ using carefully chosen multipliers that preserve the ordering of the largest noise radii. For the magnified matrix, the perturbative analysis for true signals continues to hold, and a law-of-large-numbers-type argument shows that the rescaled variance of genuine signals converges to zero. The proof proceeds conditionally on the original data matrix and relies on a refined fluctuation analysis. In contrast, the spurious eigenvalues induced by heavy-tailed noise exhibit amplified fluctuations under the multipliers. Their limiting variance does not shrink, as it is driven by the variance of the multipliers themselves. This behavior is established by linking the spurious eigenvalues to the order statistics of the magnified radii. Consequently, the asymptotic variance of spurious components remains nondegenerate, while that of true signals vanishes. This variance separation provides the theoretical foundation for our testing procedure, which consistently distinguishes genuine factors from heavy-tail-induced artifacts.

Conclusion

This paper has confronted the critical “spurious factor dilemma" in EFMs, where heavy-tailed randomness, prevalent in economic and financial data, can generate noise-induced eigenvalues that masquerade as real factors. We introduce a novel theory-based fluctuation magnification algorithm, which uniquely leverages the differential stability of real versus spurious factors under targeted perturbations: real factor signals exhibit resilience, while spurious ones betray their noisy origins through amplified volatility. This advancement significantly enhances robust factor analysis, offering a more reliable foundation for modeling and forecasting in high-dimensional, heavy-tailed distribution, thereby improving the fidelity of economic and financial decision-making.

While classical RMT often relies on assumptions of finite higher-order moments, which can be challenged by heavy-tailed distributions, this paper indicates that the conceptual framework and analytical power of RMT remain remarkably insightful and adaptable for understanding complex systems. Indeed, ongoing research continues to extend RMT's reach, developing new results and specialized techniques that successfully characterize the spectral properties of matrices with heavy-tailed entries. These advancements, by providing a deeper understanding of how eigenvalues and eigenvectors behave under non-standard conditions, are proving instrumental in developing robust methodologies capable of distinguishing real signals from noise and making reliable inferences from heavy-tailed data, thereby underscoring RMT's enduring effectiveness and its evolving role in the modern high-dimensional world.