The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
82,110 characters
The Spurious Factor Dilemma: Robust Inference in Heavy-Tailed Elliptical Factor Models
\maketitle
\begin{abstract}
Standard 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.
\end{abstract}
\textbf{KEY WORDS: Elliptical distributions; Factor models; Heavy tails; Spurious factors.
}
\tableofcontents
\addtocontents{toc}{\protect\setcounter{tocdepth}{2}}
\section{Introduction}\label{sec_intro}
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 \citep{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 \citep{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 \citep{bai2003inferential, baltagi2017identification}. Conversely, overestimating $m$ introduces noise by fitting spurious factors, which can inflate estimation variance, reduce model interpretability, and increase computational costs \citep{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 \citep{bai2002determining} and \citep{alessi2010improved}, the eigenvalue ratio tests of \citep{ahn2013eigenvalue} and \citep{lam2012factor}, and the randomization tests \citep{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 \citep{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 \citep{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 \citep{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{fig_illustration_main} 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) \citep{chamberlain1982arbitrage,baltagi2017identification}. Recently, \cite{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: \textit{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{sec_main_realsignals} and \ref{sec_main_fakesignals}), 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{sec_testing_procedures}). 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 \cite{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 \citep{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{sec_model} formally introduces the Elliptical Factor Model and outlines the key assumptions. Section \ref{sec_motivation} further illustrates the problem of spurious factors using examples. Section \ref{sec_mainresults} presents the main asymptotic theory that details the behavior of real and spurious eigenvalues. Section \ref{sec_testing_procedures} describes the fluctuation magnifier algorithm and the proposed testing and factor selection procedures. Section \ref{sec_simulation} provides simulation results, and Section \ref{sec_realdata} discusses the real data application. We provide a sketch for our proof strategy for the theoretical results in Section \ref{sketch for proof strategy}. The conclusion is offered in Section \ref{sec_conclusion}. All detailed technical proofs are deferred to Appendix \ref{appendix}.
\subsection*{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.
\section{Elliptical factor model and assumptions}\label{sec_model}
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 \citep{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:
\begin{equation}\label{eq_intro_elldistr}
\mathbf{y} \overset{d}{=}\xi \Sigma^{1/2} \mathbf{u},
\end{equation}
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:
\begin{equation}\label{eq_intro_factormodel}
\mathbf{y}_t = B \mathbf{f}_t + \mathbf{e}_t, \quad t=1, \dots, n,
\end{equation}
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$ \citep{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
\begin{align*}
\operatorname{diag}(I_m, \Sigma_{err})=
\begin{pmatrix}
I_m &0\\
0 &\Sigma_{err}
\end{pmatrix}.
\end{align*}
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=
\begin{pmatrix}
B ~~~I
\end{pmatrix}
\begin{pmatrix}
\mathbf{f}_t\\
\mathbf{e}_t
\end{pmatrix}
\sim EC_p(\mathbf{0}, \Sigma, \xi_t), \]
where
\begin{align*}
\Sigma=BB'+\Sigma_{err}\end{align*}
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:
\begin{equation}\label{eq_def_ellfactmodel}
\mathbf{y}_t = \xi_t \Sigma^{1/2} \mathbf{u}_t, \quad t=1, \dots, n,
\end{equation}
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.
\begin{assumption}[Distribution of Radius]\label{ass_xi}
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) \textit{Polynomial decay tail:}
\begin{equation}\label{ass_xi_poly}
\mathbb{P}(\xi^2 > x) = x^{-\alpha} L(x), \quad \text{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) \textit{Exponential decay tail:} For some constant $\beta > 0$ and any fixed constant $s > 0$,
\begin{equation}\label{ass_xi_exp}
\mathbb{E} e^{s\xi^{2\beta}} < \infty.
\end{equation}
\end{assumption}
\begin{remark}[Heavy Tails]
Assumption \ref{ass_xi} explicitly allows for heavy tails.
The polynomial decay case \eqref{ass_xi_poly} 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 \eqref{ass_xi_exp} 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 \citep{hu2019high}. The threshold
\begin{equation*}\label{eq_def_thresholdT}
\mathsf{T} :=
\begin{cases}
n ^{1/\alpha}\log n, & \text{if \eqref{ass_xi_poly} holds}; \\
(\log n)^{1/\beta}, & \text{if \eqref{ass_xi_exp} holds}.
\end{cases}
\end{equation*}
characterizes the typical maximal order of $\xi_t^2$. As shown later (Lemma \ref{lem_goodconfiguration}), $\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.
\end{remark}
We assume the population covariance matrix $\Sigma$ exhibits a spiked structure, which is common in factor analysis \citep{chamberlain1982arbitrage, yu2024testing}.
\begin{assumption}[Spiked Covariance Structure]\label{ass_sigma}
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$.
\end{assumption}
\begin{remark} In Assumption \ref{ass_sigma}, (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 \citep{fan2018large,yu2024testing}.
\end{remark}
\begin{remark}We want to emphasize that
although Assumption \ref{ass_sigma} 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.
\end{remark}
\begin{remark}
The main difference between Assumption \ref{ass_sigma} 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{ass_xi}, it is evident that we allow for a broad range of heavy-tailed randomness in factor noise.
\end{remark}
\section{The problem: spurious factors from heavy tails}\label{sec_motivation}
Given observations $\{\mathbf{y}_t\}_{t=1}^n$ from the EFM \eqref{eq_def_ellfactmodel}, we form the sample covariance matrix $S$ and its companion $\mathcal{S}$:
\begin{equation}\label{eq_def_samplecovariancematrices}
S := YY' = \Sigma^{1/2} U D^2 U' \Sigma^{1/2}, \quad \mathcal{S} := Y'Y = D U' \Sigma U D.
\end{equation}
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{ass_xi}), 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{fig_illustration_main} provides a toy example. Standard methods relying on gaps or thresholding (e.g., \cite{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:
\begin{equation}\label{eq_hypothesistest}
\mathbf{H}_0: \mathsf{f} = 0 \quad \text{vs.} \quad \mathbf{H}_a: \mathsf{f} > 0.
\end{equation}
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{sec_testing_procedures}, leverages the differential fluctuation of real and spurious signals under magnification perturbations.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.5\textwidth]{illl3.pdf}
\caption{Illustration of sample eigenvalues from an EFM with $p=n=1000$. The data follows a multivariate t-distribution with 4 degrees of freedom ($\alpha=2$). The population covariance is $\Sigma = \operatorname{diag}(7, 1, \dots, 1)$ ($m=1$, $\sigma_1=7$). Besides the eigenvalue near 7 (real signal), a spurious eigenvalue appears around 5.5, well-separated from the main bulk sample eigenvalues.}
\label{fig_illustration_main}
\end{figure}
\section{Asymptotic behavior of sample eigenvalues}\label{sec_mainresults}
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{ass_sigma}). As a consequence, we can rewrite
\begin{equation}\label{eq_decomposition_S}
\mathcal{S}=DU^{\prime}\Sigma_1UD+DU^{\prime}\Sigma_2UD.
\end{equation}
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.
\subsection{Asymptotic behavior of real factors}\label{sec_main_realsignals}
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.
\begin{theorem}[First Order Approximation]\label{lem_realsignal_preratio}
Under Assumptions \ref{ass_xi} and \ref{ass_sigma}, we have:
\begin{itemize}
\item[(1)] for polynomial decay tail \eqref{ass_xi_poly} with $\alpha\in(2,+\infty)$ or exponential decay tail \eqref{ass_xi_exp} in Assumption \ref{ass_xi},
\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*}
\item[(2)] for polynomial decay tail \eqref{ass_xi_poly} with $\alpha\in(1,2]$ in Assumption \ref{ass_xi},
\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}
\end{theorem}
\begin{remark}
Theorem \ref{lem_realsignal_preratio} 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{lem_realsignal_preratio} 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{lem_realsignal_preratio} are far from ideal, especially if our goal is to determine the asymptotic distribution of $\lambda_i$.
\end{remark}
To obtain a more precise description for the locations of $\lambda_i$, we employ the methodology from \cite{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
\begin{equation}\label{eq_def_theta/sigma}
\frac{\theta_i}{\sigma_i}=\Big(1-\frac{1}{p\theta_i}\sum_{k=1}^{p-m}\frac{\sigma_{m+k}}{1-\sigma^{-1}_{i}\sigma_{m+k}}\Big)^{-1},\quad \theta_i\in[\sigma_i,2\sigma_i],\quad 1\le i\le m.
\end{equation}
\cite{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
\begin{equation}\label{eq_def_zetai}
\zeta_i=\frac{1}{p}\sum_{j=1}^n\xi_j^2\Big(1-\frac{\xi^2_j}{p\theta_i}\sum_{k=1}^{p-m}\frac{\sigma_{m+k}}{1-\theta_i^{-1}\sigma_{m+k}\zeta_i}\Big)^{-1},~\zeta_i\in[\frac{1}{p}\operatorname{tr}D^2,\frac{2}{p}\operatorname{tr}D^2].
\end{equation}
The existence and uniqueness of $\theta_i, \zeta_i$ are guaranteed under our assumptions for large $p, n$.
\begin{lemma}[Properties of $\theta_i, \zeta_i$]\label{lem_realsignal_secondratio1}
Under Assumptions \ref{ass_xi} and \ref{ass_sigma}, 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 \eqref{ass_xi_poly} with $\alpha\in(2,+\infty)$ or exponential decay tail \eqref{ass_xi_exp} in Assumption \ref{ass_xi},
\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 \eqref{ass_xi_poly} 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*}
\end{lemma}
Using these observations, we obtain a more precise characterization of $\lambda_i$. Let $\mathbf{e}_k$ be the $k$-th standard basis vector.
\begin{theorem}[Second Order Approximation]\label{thm_realsignal_secondratio2}
Under Assumptions \ref{ass_xi} and \ref{ass_sigma}, we have that: \begin{itemize}
\item[(1)] for for polynomial decay tail \eqref{ass_xi_poly} with $\alpha\in(2,+\infty)$ or exponential decay tail \eqref{ass_xi_exp} in Assumption \ref{ass_xi},
\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*}
\item[(2)] for polynomial decay tail \eqref{ass_xi_poly} 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}
\end{theorem}
Then, related to the results in Lemma \ref{lem_realsignal_secondratio1}, Theorem \ref{thm_realsignal_secondratio2} 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{lem_realsignal_preratio} especially for $\alpha\in(2,+\infty)$ and exponential decay tail. On the other hand, we notice that the error rate in Theorem \ref{thm_realsignal_secondratio2} 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{ass_xi}, the following central limit theorem holds for real signals from Theorem \ref{thm_realsignal_secondratio2},
\begin{theorem}\label{thm_realsignal_clt}
Assume \eqref{ass_xi_poly} with $\alpha\in(2,+\infty)$ or \eqref{ass_xi_exp} in Assumption \ref{ass_xi}, also assume Assumption \ref{ass_sigma}. 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}
\end{theorem}
\begin{remark}
It 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.
\end{remark}
When $\xi^2$ exhibits serious heavy-tailed decay with $\alpha\in(1,2]$, the asymptotic normality in Theorem \ref{thm_realsignal_clt} 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.
\begin{theorem}\label{thm_realsignal_extremeheavy}
Suppose \eqref{ass_xi_poly} holds with $\alpha\in(1,2]$ in Assumption \ref{ass_xi}. 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}
\end{theorem}
\subsection{Asymptotic behavior of spurious factors}\label{sec_main_fakesignals}
We first prepare several notations. Recalling the decomposition of $\mathcal{S}$ in \eqref{eq_decomposition_S}, we denote the matching matrices removing the spike structure from $\mathcal{S}$ (or $S$) as
\begin{equation}\label{eq_def_S2}
S_2:=\Sigma^{1/2}_2UD^2U^{\prime}\Sigma_2^{1/2},\quad \mathcal{S}_2:=DU^{\prime}\Sigma_2UD.
\end{equation}
Denote $\lambda_{i,2}$ as the $i$-th largest eigenvalue of $S_2$ or $\mathcal{S}_2$. Denote the quantity $\mathsf{q}$ as
\begin{equation}\label{eq_def_q}
\mathsf{q}:=
\begin{cases}
n^{1/\alpha-\epsilon}\quad &\text{if \eqref{ass_xi_poly} in Assumption \ref{ass_xi} holds};\\
1 \quad &\text{if \eqref{ass_xi_exp} in Assumption
\ref{ass_xi} holds},
\end{cases}
\end{equation}
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$,
\begin{lemma}\label{lem_fakesignal_compare_nonspikes}
Under Assumptions \ref{ass_xi} and \ref{ass_sigma}, 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}
\end{lemma}
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,
\begin{theorem}\label{thm_fakesignal_limit}
Under Assumptions \ref{ass_xi} and \ref{ass_sigma}, 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$.
\end{theorem}
As a result of Theorem \ref{thm_fakesignal_limit}, 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$.
\begin{remark}
From Lemma \ref{lem_fakesignal_compare_nonspikes}, we observe that the spurious noise signals closely approximate the leading eigenvalues of $\mathcal{S}_2$. Furthermore, Theorem \ref{thm_fakesignal_limit} and Lemma \ref{lem_goodconfiguration} 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.
\end{remark}
\section{Detection via Fluctuation Magnification}\label{sec_testing_procedures}
In this section, we develop testing strategies for the problem \eqref{eq_hypothesistest}, guided by the theoretical framework established in Sections \ref{sec_main_realsignals} and \ref{sec_main_fakesignals}. 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.
\subsection{Fluctuation magnification for $S$}\label{sec_testing_bootstrap}
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
\begin{equation*}
\widetilde{Y}_j:=Y\cdot W^{1/2}_j,\quad W_j=\operatorname{diag}(w_t)_{1\le t\le n},
\end{equation*}
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
\begin{gather*}
\widetilde{S}^{(j)}:=\widetilde{Y}_j\widetilde{Y}_j^{\prime}=\Sigma^{1/2}UDW_jD^{\prime}U^{\prime}\Sigma^{1/2}\quad \widetilde{\mathcal{S}}^{(j)}:=W_j^{1/2}DU^{\prime}\Sigma UDW_j^{1/2}\\
\widetilde{\mathcal{S}}_2^{(j)}:=W_j^{1/2}DU^{\prime}\Sigma_2UDW_j^{1/2}.
\end{gather*}
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{thm_realsignal_clt} and Theorem \ref{thm_realsignal_extremeheavy}. Applying these two theorems to each $\widetilde{S}^{(j)}$ for $1\le j\le K$, we have the following proposition.
\begin{proposition}\label{lem_realsignals_bootstrapped}
For all $1\le j\le K$ and $1\le i\le m$,
\begin{itemize}
\item [(1)] if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(2,+\infty)$ or exponential decay tail as in Theorem \ref{thm_realsignal_clt}, then we have that
\begin{equation*}\label{eq_test_realsignal_clt}
\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*}
\item [(2)] if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(1,2]$ as in Theorem \ref{thm_realsignal_extremeheavy}, 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}$.
\end{proposition}
As a consequence, we have the following proposition.
\begin{proposition}\label{col_realsignals_bootstrapped}
Under the assumptions of Proposition \ref{lem_realsignals_bootstrapped}, for $1\le i\le m$ and sufficiently large $K$,
\begin{itemize}
\item [(1)] if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(2,+\infty)$ or exponential decay tail as in Theorem \ref{thm_realsignal_clt}, 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*}
\item [(2)] if $\xi^2$ exhibits polynomial decay tail with $\alpha\in(1,2]$ as in Theorem \ref{thm_realsignal_extremeheavy}, 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}
\end{proposition}
Next, we turn to the spurious factors perturbed by some magnifier matrices. Lemma \ref{lem_fakesignal_compare_nonspikes} and Theorem \ref{thm_fakesignal_limit} indicate that for $m+1\le i\le \hat{o}$,
\begin{equation*}
\big|({\lambda_{i}^{(j)}-\lambda_{i-m,2}^{(j)}})/{(\xi^2w)_{(i-m)}}\big|=\mathrm{o}_{\mathbb{P}}(n^{-1/2+\epsilon}),
\mbox{ with }
\lambda_{i-m,2}^{(j)}=(\bar{\sigma}+\mathrm{o}_{\mathbb{P}}(1))\cdot (\xi^2w)_{(i-m)}.
\end{equation*}
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{lem_realsignals_bootstrapped}. Then, it is easy to find that $(\xi^2w)_{(1)}$ is uniformly distributed on $[a\xi^2_{(1)},b\xi^2_{(1)}]$ with
\begin{gather*}
\mathbb{E}\big[(\xi^2w)_{(1)}\big]=\mathbb{E}[\xi^2_{(1)}w_1]=\frac{a+b}{2}\cdot\xi^2_{(1)},\\
\mathbb{E}\big[(\xi^2w)^2_{(1)}\big]=\mathbb{E}[\xi^4_{(1)}w_1^2]=\frac{\xi^4_{(1)}}{3}\cdot(b^2+ab+a^2),\\
\mathbb{E}[(\xi^2w)^4_{(1)}]=\mathbb{E}[\xi^8_{(1)}w_1^4]=\frac{\xi^8_{(1)}}{5}\cdot (b^4+b^3a+b^2a^2+ba^3+a^4).
\end{gather*}
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.
\begin{proposition}\label{lem_fakesignals_bootstrapped}
Under the assumptions in Theorem \ref{thm_fakesignal_limit} 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$.
\end{proposition}
\begin{remark}\label{rem_bootstrap}
Several 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{lem_fakesignals_bootstrapped}.} 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{col_realsignals_bootstrapped}} and spurious factors in {Proposition \ref{lem_fakesignals_bootstrapped}} 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.
\end{remark}
\subsection{Detection algorithm}\label{sec_testing_algorithm}
Inspired by the theoretical observations in Section \ref{sec_testing_bootstrap}, we propose the following statistics.
\begin{equation}\label{eq_testing_statistics}
\mathbb{T}_i:=\frac{1}{K}\sum_{j=1}^K\Big(\frac{\lambda_i^{(j)}}{\hat{\theta}_i}-1\Big)^2.
\end{equation}
Recap the test problem \eqref{eq_hypothesistest} in Section \ref{sec_motivation}. 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 \cite{onatski2010determining}, thus refining the selection of common factors.
\textbf{1. First step detection.}
We now proceed to the first step. Under Assumption \ref{ass_sigma}, it suffices to focus on the largest spurious signal under the alternative hypothesis $\mathbf{H}_a$. As noted in Remark \ref{rem_bootstrap}, it is crucial to choose an appropriate distribution for $w$ to ensure the validity of {Proposition \ref{lem_fakesignals_bootstrapped}}. 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{thm_fakesignal_limit}, while Lemma \ref{lem_realsignal_secondratio1} and Theorem \ref{thm_realsignal_clt} support the second. The algorithm is summarized in Algorithm \ref{alg_firstround_bootstrap}, which gives preliminary estimations for the number of spurious factors and the location of the largest potential spurious signal.
\begin{algorithm}[htbp]
\caption{Detection for spurious signals}\label{alg_firstround_bootstrap}
\normalsize
\begin{flushleft}
\noindent{\bf Inputs:} Pre-given integer $o$, original data matrix $Y$, sufficiently large $K\equiv K(n)$.
\noindent{\bf Step One:} Calculate the eigenvalues of $S=YY^{\prime}$ as $\lambda_1\ge\lambda_2\ge\dots\ge\lambda_{n}\ge0$. For $i=1,\dots,o$, we generate $w_{i}\sim U[a_i,b_i]$ satisfying: (1). $(a_i+b_i)/2=1$; (2). $b_i/a_i\le\lambda_i/\lambda_{i+1}$; (3). $b_i\log n\times\lambda_{i}/\lambda_{i-1}<1$, where we set $\lambda_{0}=\log^2 n\lambda_1$.
\noindent{\bf Step Two:} For each $1\le i\le o$,
generate $K$ i.i.d. $n\times n$ diagonal matrices $W^{(j)}_i, j=1,\dots,K$ with entries independently sampled from $w_i\sim\mathrm{U}[a_i,b_i]$. Compute the associated matrices after fluctuation magnification $\widetilde{S}_{i}^{(j)}=YW^{(j)}_iY^{\prime}$ and the sequential eigenvalues $\lambda_i^{(j)}$, where $\lambda_i^{(j)}$ is the $i$-th largest eigenvalue of $\widetilde{S}_i^{(j)}$.
\noindent{\bf Step Three:} Compute $\mathbb{T}_i=K^{-1}\sum_{j=1}^K(\lambda_i^{(j)}-K^{-1}\sum_{j=1}^K\lambda_i^{(j)})^2/(K^{-1}\sum_{j=1}^K\lambda_i^{(j)})^2$ for each $1\le i\le o$. Set the thresholds $\{\mathsf{L}_i\}_{1\le i\le o}$ according to $\{\lambda_i\}_{1\le i\le o}$ as $\mathsf{L}_i:=\log^2 n/n\times\mathbbm{1}(\lambda_i\le n^{1/2})+\log^2 n/n^{3/2-\tau_i}\times\mathbbm{1}(\lambda_i>n^{1/2})$ where $\tau_i:=\log\lambda_i/\log n$.
Count preliminary estimators ${\mathsf{f}}^*=\#\{1\le i\le o:\mathbb{T}_i>\mathsf{L}_i\}$ and $r^*:=\min\{\{1\le i\le o:\mathbb{T}_i>\mathsf{L}_i\}\cup\{o+1\}\}.$
\noindent{\bf Output:} ${\mathsf{f}}^*$ and $r^*$.
\end{flushleft}
\end{algorithm}
\begin{remark}
Several remarks are in order. First, the output $r^*$ in Algorithm \ref{alg_firstround_bootstrap}
\[
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{sec_testing_bootstrap}, it holds that
\begin{gather}\label{eq_firstround_consistency}
\lim_{n,K\rightarrow\infty}\mathbb{P}(\mathbb{T}_i\le\mathsf{L}_i)=1,\;\text{under $\mathbf{H}_0$, ~~for $1\le i\le m$;}\\
\lim_{n,K\rightarrow\infty}\mathbb{P}(\mathbb{T}_i>\mathsf{L}_i)=1,\;\text{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{lem_fakesignals_bootstrapped}.
Third, the magnifier $w$ can be extended to other distributions, such as a Bernoulli-type random variable satisfying
\begin{equation}
w=
\begin{cases}
a & \text{with probability } p, \\
b & \text{with probability } 1-p,
\end{cases}
\end{equation}
under conditions analogous to those specified in Algorithm \ref{alg_firstround_bootstrap}. 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{ass_xi} and \ref{ass_sigma}, 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.
\end{remark}
\textbf{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 \cite{onatski2010determining} and \cite{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{alg_secondround_estrealsignal} 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.
\begin{algorithm}[htbp]
\caption{Estimation for the number of common factors}\label{alg_secondround_estrealsignal}
\normalsize
\begin{flushleft}
\noindent\textbf{Inputs:} $r^*$, ${\mathsf{f}}^*$ and the leading $o+1$ eigenvalues of $S$: $\lambda_1\ge\lambda_2\ge\dots\ge\lambda_{o+1}$.
\noindent\textbf{Step One:} Compute an initial estimate $k$ using a consistent factor‑number detection method (e.g., the ED estimator of \cite{onatski2010determining} or the DDPA+ procedure of \cite{Dobriban}) applied to $S$.
\noindent\textbf{Step Two:} Set $\widehat{r}=k$ and $\widehat{\mathsf{f}}=0$ if $k < r^*$ (no spurious signal occurs); otherwise, set $\widehat{r}=r^*-1$ (a spurious signal is detected at $r^* \le k$, so the estimate is truncated to $r^*-1$) and $\widehat{\mathsf{f}}={\mathsf{f}}^*$.
\noindent\textbf{Output:} The estimated number of common factors $\widehat{r}$, and reject $\mathbf{H}_0$ if $\widehat{\mathsf{f}}>0$.
\end{flushleft}
\end{algorithm}
\begin{lemma}\label{lem_secondround_consistency}
Let $\widehat{r}$ be obtained by Algorithms \ref{alg_firstround_bootstrap} and \ref{alg_secondround_estrealsignal}. Under the assumptions in Propositions \ref{col_realsignals_bootstrapped} and \ref{lem_fakesignals_bootstrapped}, we have that
\begin{equation}
\lim_{n\rightarrow\infty}\mathbb{P}^*(\widehat{r}=m)=1.
\end{equation}
\end{lemma}
\begin{remark}
Several remarks are in order. First, Algorithm \ref{alg_secondround_estrealsignal} incorporates the methods developed in \cite{onatski2010determining} and \cite{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 \cite{onatski2010determining} or \cite{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{lem_secondround_consistency}. Second, the choice of the specific methods in \cite{onatski2010determining} and \cite{Dobriban} for Algorithm \ref{alg_secondround_estrealsignal} 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 \cite{bai2018consistency}.
\end{remark}
\section{Simulation}\label{sec_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 \eqref{eq_intro_elldistr} 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$.
\subsection{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 \eqref{eq_intro_elldistr}, 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$:
\begin{itemize}
\item[(i)] $\Sigma_{\mathrm{I}} = \operatorname{diag}\{16, 8, 1, \ldots, 1\}$
\item[(ii)] $\Sigma_{\mathrm{II}} = \operatorname{diag}\{24, 16, 8, 1, \ldots, 1\}$.
\end{itemize}
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 \cite{onatski2010determining} (denoted as ``Onta") and the DDPA+ procedure of \cite{Dobriban} (denoted as ``DDPA+"). Their versions enhanced by the Fluctuation Magnification technique proposed in Section~\ref{sec_testing_algorithm} 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{alg_firstround_bootstrap} and $K=1000$ in Algorithm \ref{alg_secondround_estrealsignal}, which yielded sufficiently good results. Detailed experimental results are presented in Figure \ref{eig4.3}, and Tables \ref{tab:af_comparison1} and \ref{tab:af_comparison2}.
Specifically, Figure~\ref{eig4.3} 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{alg_firstround_bootstrap}; the statistics corresponding to spurious signals are markedly larger than those associated with real signals. These observations highlight the effectiveness of Algorithm~\ref{alg_firstround_bootstrap} in detecting spurious signals.
\begin{figure*}[htbp]
\centering
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{eig4.3.eps}
\centerline{(a)}
\end{minipage}
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{onta4.3.eps}
\centerline{(b)}
\end{minipage}
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{ou4.3.eps}
\centerline{(c)}
\end{minipage}
\caption{Detecting behavior under case (i) with $\mathbf{y}\sim t(4.3)$. Screening and ``On" in (a) and (b) indicate that $k = 3$, suggesting the emergence of spurious signals that lead to an overestimation of the true value, $2$. Algorithm \ref{alg_firstround_bootstrap} accurately identifies the initial spurious signal at $3$ by detecting abnormal variance. }
\label{eig4.3}
\end{figure*}
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{tab:af_comparison1} and~\ref{tab:af_comparison2} 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.
\begin{table}[htbp]
\centering
\small
\begin{tabular}{@{} l l l ccc ccc @{}}
\toprule
\multicolumn{3}{c}{} & \multicolumn{3}{c}{\textbf{Without MF}} & \multicolumn{3}{c}{\textbf{With MF}} \\
\cmidrule(lr){4-6} \cmidrule(l){7-9}
\textbf{Dist.} & \textbf{ESS} & \textbf{Method} & \textbf{Under} & \textbf{Over} & \textbf{Right} & \textbf{Under} & \textbf{Over} & \textbf{Right} \\
\midrule
\multirow{2}{*}{$t(4.3)$} & \multirow{2}{*}{499} & \textbf{Onta} & 1.8\% & 10.2\% & 88.0\% & 3.0\% & 2.4\% & 94.6\% \\
& & \textbf{DDPA+} & 0.0\% & 12.4\% & 87.6\% & 1.6\% & 4.6\% & 93.8\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(4.8)$} & \multirow{2}{*}{499} & \textbf{Onta} & 0.6\% & 2.4\% & 97.0\% & 0.8\% & 0.2\% & 99.0\% \\
& & \textbf{DDPA+} & 0.0\% & 2.8\% & 97.2\% & 0.2\% & 0.8\% & 99.0\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(5.3)$} & \multirow{2}{*}{500} & \textbf{Onta} & 0.0\% & 0.8\% & 99.2\% & 0.0\% & 0.0\% & 100\% \\
& & \textbf{DDPA+} & 0.0\% & 1.0\% & 99.0\% & 0.0\% & 0.0\% & 100\% \\
\bottomrule
\end{tabular}
\caption{Comparison of Onta and DDPA+ methods with and without the MF tool Under Case (i). The table shows the percentage of underestimation (Under), overestimation (Over), and correct estimation (Right) for three $t$-distributions with varying degrees of freedom. (Dist. = Distribution, ESS = Effective Sample Size}
\label{tab:af_comparison1}
\end{table}
\begin{table}[htbp]
\centering
\small
\begin{tabular}{@{} l l l ccc ccc @{}}
\toprule
\multicolumn{3}{c}{} & \multicolumn{3}{c}{\textbf{Without MF}} & \multicolumn{3}{c}{\textbf{With MF}} \\
\cmidrule(lr){4-6} \cmidrule(l){7-9}
\textbf{Dist.} & \textbf{ESS} & \textbf{Method} & \textbf{Under} & \textbf{Over} & \textbf{Right} & \textbf{Under} & \textbf{Over} & \textbf{Right} \\
\midrule
\multirow{2}{*}{$t(4.3)$} & \multirow{2}{*}{498} & \textbf{Onta} & 0.8\% & 9.4\% & 89.8\% & 2.6\% & 9.0\% & 95.4\% \\
& & \textbf{DDPA+} & 0.0\% & 13.1\% & 86.9\% & 2.2\% & 5.0\% & 92.8\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(4.8)$} & \multirow{2}{*}{500} & \textbf{Onta} & 0.4\% & 3.0\% & 96.6\% & 0.4\% & 0.0\% & 99.6\% \\
& & \textbf{DDPA+} & 0.2\% & 4.4\% & 95.4\% & 0.2\% & 1.2\% & 98.6\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(5.3)$} & \multirow{2}{*}{500} & \textbf{Onta} & 0.0\% & 1.2\% & 98.8\% & 0.0\% & 0.0\% & 100\% \\
& & \textbf{DDPA+} & 0.0\% & 1.6\% & 98.4\% & 0.0\% & 0.2\% & 99.8\% \\
\bottomrule
\end{tabular}
\caption{Comparison of Onta and DDPA+ methods with and without the MF tool Under Case (i). The table shows the percentage of underestimation (Under), overestimation (Over), and correct estimation (Right) for three $t$-distributions with varying degrees of freedom. (Dist. = Distribution, ESS = Effective Sample Size}
\label{tab:af_comparison2}
\end{table}
From Figure \ref{eig4.3}, and Tables \ref{tab:af_comparison1} and \ref{tab:af_comparison2}, we can draw the following conclusions.
\begin{itemize}
\item 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.
\item When spurious signals emerge, the fluctuation magnification algorithm (Algorithm \ref{alg_firstround_bootstrap}) 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.
\item 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.
\item It is noteworthy that a higher underestimation rate occurring under low degrees of freedom is, paradoxically, a desirable outcome. From Tables \ref{tab:af_comparison1} and \ref{tab:af_comparison2}, 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.
\end{itemize}
\subsection{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 \eqref{eq_intro_elldistr}, 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:
\begin{itemize}
\item[(iii)] $\Sigma_{\mathrm{III}} = \operatorname{diag}\{240, 120, 1, \ldots, 1\}$
\item[(iv)] $\Sigma_{\mathrm{IV}}= \operatorname{diag}\{360, 240, 120, 1, \ldots, 1\}$
\end{itemize}
Figure~\ref{eig2.5} illustrates the performance of the methods for a representative case (Case (iii)) using the ``Onta" estimator under $t(2.5)$ as an example.
\begin{figure*}[htbp]
\centering
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{eig2.5.eps}
\centerline{(a)}
\end{minipage}
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{onta2.5.eps}
\centerline{(b)}
\end{minipage}
\begin{minipage}{0.32\linewidth}
\vspace{1pt}
\includegraphics[width=\textwidth]{ou2.5.eps}
\centerline{(c)}
\end{minipage}
\caption{Detecting behavior under case (iii) with $\mathbf{y}\sim t(2.5)$. Screening and ``On" in (a) and (b) indicate $\widehat{r} = 3$, suggesting the emergence of spurious signals that lead to an overestimation of the true value, $2$. Algorithm \ref{alg_firstround_bootstrap} accurately identifies the initial spurious signal at $3$ by detecting abnormal variance. }
\label{eig2.5}
\end{figure*}
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.
\begin{table}[htbp]
\centering
\small
\begin{tabular}{@{} l l l ccc ccc @{}}
\toprule
\multicolumn{3}{c}{} & \multicolumn{3}{c}{\textbf{Without MF}} & \multicolumn{3}{c}{\textbf{With MF}} \\
\cmidrule(lr){4-6} \cmidrule(l){7-9}
\textbf{Dist.} & \textbf{ESS} & \textbf{Method} & \textbf{Under} & \textbf{Over} & \textbf{Right} & \textbf{Under} & \textbf{Over} & \textbf{Right} \\
\midrule
\multirow{2}{*}{$t(2.5)$} & \multirow{2}{*}{490} & \textbf{Onta} & 11.4\% & 30.2\% & 58.4\% & 13.3\% & 2.4\% & 84.3\% \\
& & \textbf{DDPA+} & 0.4\% & 98.4\% & 1.2\% & 2.2\% & 11.1\% & 86.7\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(3.0)$} & \multirow{2}{*}{499} & \textbf{Onta} & 1.2\% & 9.1\% & 89.7\% & 1.2\% & 0.2\% & 98.6\% \\
& & \textbf{DDPA+} & 0.0\% & 94.4\% & 5.6\% & 0.2\% & 13.6\% & 86.2\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(3.5)$} & \multirow{2}{*}{500} & \textbf{Onta} & 0.4\% & 2.8\% & 96.8\% & 0.4\% & 0.0\% & 99.6\% \\
& & \textbf{DDPA+} & 0.0\% & 66.8\% & 33.2\% & 0.0\% & 15.8\% & 84.2\% \\
\bottomrule
\end{tabular}
\caption{Comparison of Onta and DDPA+ methods with and without the MF tool Under Case (iii). The table shows the percentage of underestimation (Under), overestimation (Over), and correct estimation (Right) for three $t$-distributions with varying degrees of freedom. (Dist. = Distribution, ESS = Effective Sample Size)}
\label{tab:af_comparison3}
\end{table}
\begin{table}[htbp]
\centering
\small
\begin{tabular}{@{} l l l ccc ccc @{}}
\toprule
\multicolumn{3}{c}{} & \multicolumn{3}{c}{\textbf{Without MF}} & \multicolumn{3}{c}{\textbf{With MF}} \\
\cmidrule(lr){4-6} \cmidrule(l){7-9}
\textbf{Dist.} & \textbf{ESS} & \textbf{Method} & \textbf{Under} & \textbf{Over} & \textbf{Right} & \textbf{Under} & \textbf{Over} & \textbf{Right} \\
\midrule
\multirow{2}{*}{$t(2.5)$} & \multirow{2}{*}{491} & \textbf{Onta} & 14.1\% & 30.8\% & 55.2\% & 15.5\% & 5.1\% & 79.4\% \\
& & \textbf{DDPA+} & 0.2\% & 97.4\% & 2.4\% & 1.4\% & 9.0\% & 89.6\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(3.0)$} & \multirow{2}{*}{499} & \textbf{Onta} & 5.4\% & 6.6\% & 88.0\% & 5.4\% & 0.2\% & 94.4\% \\
& & \textbf{DDPA+} & 0.0\% & 93.6\% & 6.4\% & 0.0\% & 12.4\% & 87.6\% \\
\addlinespace[0.05cm]
\multirow{2}{*}{$t(3.5)$} & \multirow{2}{*}{500} & \textbf{Onta} & 1.4\% & 3.2\% & 95.4\% & 1.4\% & 0.0\% & 98.6\% \\
& & \textbf{DDPA+} & 0.0\% & 64.8\% & 35.2\% & 0.0\% & 13.0\% & 87.0\% \\
\bottomrule
\end{tabular}
\caption{Comparison of Onta and DDPA+ methods with and without the MF tool Under Case (iv). The table shows the percentage of underestimation (Under), overestimation (Over), and correct estimation (Right) for three $t$-distributions with varying degrees of freedom. (Dist. = Distribution, ESS = Effective Sample Size)}
\label{tab:af_comparison4}
\end{table}
From Figure \ref{eig2.5} and Tables \ref{tab:af_comparison3} and \ref{tab:af_comparison4}, 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.
\section{Real data analysis}\label{sec_realdata}
In this section, the proposed methods are applied to the real-world FRED-MD dataset, which was previously studied in \cite{XIA2017235} and \cite{YU2019104543}. This dataset, introduced in \cite{McCracken01102016}, is publicly available on the St. Louis Fed's website: \texttt{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 \cite{XIA2017235} and slightly smaller than the estimate $\widehat{r} = 4$ proposed by \cite{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 \cite{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{First order} shows $\mathbb{G}_i$ (``On" estimation), $i=1,\ldots,10$ of every 48 months data matrices.
\begin{table}[htpb]
\caption{First-order estimation $\mathbb{G}_i$ for $i=1,\ldots,10$}
\centering
\begin{tabular}{lllllllllllll}
\hline
Date&1&2&3&4&5&6&7&8&9&10\\
\hline
1992-1996&1.20 & 53.52 & 1.32 & 1.05 & 2.96 & 2.72 & 2.30 & 0.38 & 2.16 & 1.88 \\
1996-2000&1.95 & 28.46 & 1.34 & 3.95 & 7.38 & 0.72 & 1.79 & 0.54 & 1.01 & 9.45 \\
2000-2004&2.21 & 22.65 & 2.48 & 0.46 & 26.72 & 1.10 & 0.67 & 9.77 & 0.12 & 4.93 \\
2004-2008&2.91 & 54.15 & 0.82 & 6.24 & 0.10 & 13.10 & 5.74 & 0.92 & 0.22 & 2.83 \\
2008-2012&2.64 & 10.07 & 2.94 & 1.33 & 3.08 & 3.69 & 0.92 & 1.31 & 9.23 & 0.08\\
2012-2016&3.00 & 68.32 & 0.47 & 6.57 & 1.28 & 2.52 & 2.06 & 0.59 & 2.71 & 1.05 \\
2016-2020&2.72 & 155.09 & 0.23 & 6.80 & 0.38 & 11.58 & 0.38 & 2.20 & 4.96 & 0.69\\
2020-2004&7.98 & 18.41 & 1.38 & 4.44 & 0.65 & 2.18 & 2.94 & 0.49 & 9.15 & 0.25\\
\hline
\end{tabular}
\label{First order}
\end{table}
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{f2000-2004} 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{f2000-2004}, 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{f2004-2008}, we observe an unusual pattern that diverges from the previous trends.
\begin{figure*}[htbp]
\centering
\subfloat[Variance on 2000-2004]{\label{f2000-2004}
\centering
\includegraphics[width=0.315\textwidth]{2000-2004.pdf}
\captionsetup{font=small}
}
\hfill
\subfloat[Variance on 2004-2008]{\label{f2004-2008}
\centering
\includegraphics[width=0.315\textwidth]{2004-2008.pdf}
\captionsetup{font=small}
}
\hfill
\subfloat[Factors change every 4 years]{\label{fperiod}
\centering
\includegraphics[width=0.315\textwidth]{period.pdf}
\captionsetup{font=small}
}
\caption{
(a) Fluctuation in variance on period ``2000-2004".
(b) Fluctuation in variance on period ``2004-2008".
(c) Change in the number of factors from 1992 to 2024 every 48 months.
}
\label{fig:master}
\end{figure*}
From Figure \ref{f2004-2008}, 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{fperiod} shows the change in the number of factors from 1992 to 2024.
\section{Sketch of proof strategy}\label{sketch for 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.
\begin{itemize}
\item
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{lem_fakesignal_compare_nonspikes} and Theorem \ref{thm_fakesignal_limit}.
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{lem_solutionsystem} and \ref{thm_prf_fakesignals_locallaw}. 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 \eqref{eq_def_S2}. 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{thm_fakesignal_limit}, $\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{thm_prf_fakesignals_eigenvaluerigidity}. Specifically, as we define a determinant function $M(\cdot)$ as in \eqref{eq_defnmlambda}. $\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.
\item We next analyze the leading eigenvalues of $S$ associated with divergent factor signals. Our approach builds on the perturbative framework developed in \cite{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 \eqref{eq_def_zetai}, 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.
\item 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.
\end{itemize}
\section{Conclusion}\label{sec_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.