Asymptotic locations of bounded and unbounded eigenvalues of sample correlation matrices of certain factor models
– application to a components retention rule

In multivariate analysis (Anderson, 2003; Fujikoshi et al., 2011; Muirhead, 2009), covariance matrices are important objects, for example, in principal component analysis (PCA). Motivated by this, research on sample covariance matrices has a long history and an abundance of results have been established, e.g. various Wishart distributions.

Several well-known methods in multivariate analysis, however, become inefficient or even misleading when the data dimension $N$ is large. To deal with such large-dimensional data, a novel approach in asymptotic statistics has been developed where the data dimension $N$ is no longer fixed but tends to infinity together with the sample size $T$. The feature of this limiting regime are discussed in (Yao et al., 2015, Section 1.1) based on real datasets such as portfolio, climate survey, speech analysis, face recognition, microarrays, signal detection, etc. In this paper, we assume a high-dimensional regime in which the ratio $N/T$, of dimension $N$ to sample size $T$, converges to a positive constant $c$.

One of milestone work about the spectral properties of sample covariance matrices in this limiting regime was the discovery of Marčenko-Pastur distribution, and the determination (Yin et al., 1988; Bai & Yin, 1993) of the largest and the smallest eigenvalues of sample covariance matrices $\mathbf{S}$ for i.i.d. data under the finite fourth moment condition. Baik & Silverstein (2006) studied asymptotic locations of the eigenvalues of $\mathbf{S}$ of a spiked eigenvalue model (Johnstone, 2001), but the largest eigenvalues of the population covariance matrix $\bm{\Sigma}$ are bounded. We refer the reader to (Bai & Silverstein, 2010; Paul & Aue, 2014; Yao et al., 2015).

In some cases, for example, when data are measured in different units, it is more appropriate to utilize sample correlation matrices (Jolliffe, 2002; Jolliffe & Cadima, 2016; Johnson & Wichern, 2007) because sample correlation matrices $\mathbf{C}$ are invariant under scaling and shifting. Moreover, the assumption of i.i.d. data is not very acceptable; Practitioners often hope to be in the presence of a curious covariance structure. However, the asymptotic spectral properties of a sample correlation matrix $\mathbf{C}$ in the high-dimensional regime have not been sufficiently investigated, compared to sample covariance matrices $\mathbf{S}$.

We review typical asymptotic results of the spectral properties of sample correlation matrices $\mathbf{C}$ here. For i.i.d. case, under the finite fourth moment condition, Jiang (2004) showed that the extreme eigenvalues of the sample correlation matrix $\mathbf{C}$ converge almost surely to $(1\pm\sqrt{c})^{2}$, and Bai & Zhou (2008) showed that the limiting spectral distribution is the standard Marčenko-Pastur distribution under a weaker assumption.

In a class of spiked models, Morales-Jimenez et al. (2021) derived asymptotic first-order and distributional results for spiked eigenvalues and eigenvectors of sample correlation matrices $\mathbf{C}$. More specifically, they found that the first-order spectral properties of sample correlation matrices $\mathbf{C}$ match those of sample covariance matrices $\mathbf{S}$, whilst their asymptotic distributions can differ significantly.

El Karoui (2009) revealed that the first-order asymptotic behavior of the spectra of $\mathbf{C}$ is similar to that of $\mathbf{S}$, for unit variance data of a general linear model, except that this similarity requires the boundedness of the population correlation matrix $\mathbf{R}$, a condition not met by some factor models.

The problem is that the boundedness of the population covariance matrices $\bm{\Sigma}$ is not always satisfied in econometrics, finance (Chamberlain & Rothschild, 1983; Bai & Ng, 2002), genomics, and stationary long-memory processes. One of the features of the unbounded population covariance matrices $\bm{\Sigma}$ is the consistency of the eigenvectors (Yata & Aoshima, 2013; Koltchinskii & Lounici, 2017; Wang & Fan, 2017), although the asymptotic location and the fluctuation of the largest unbounded eigenvalues of $\mathbf{C}$ are not available yet even for the equi-correlated normal population, which has the simplest unbounded population covariance matrix $\bm{\Sigma}$.

Unbounded covariance/correlation matrices in high-dimensional problems are studied with a spiked model (Yata & Aoshima, 2013), a time-series model (Merlevède et al., 2019) or a factor model (Cai et al., 2020; Wang & Fan, 2017). Following the latter articles for unbounded sample covariance matrices, we study factor models for unbounded sample correlation matrices.

Our target model is a $K$-factor model

$$\mathbf{X}=\begin{bmatrix}\bm{\mu}&\dots&\bm{\mu}\end{bmatrix}+\mathbf{L}\mathbf{F}+\mathbf{\Lambda}\mathbf{\Psi}$$

with $\bm{\mu}\in{\mathbb{R}}^{N}$, $\mathbf{L}\in{\mathbb{R}}^{N\times K}$, $\mathbf{\Lambda}\in{\mathbb{R}}^{N\times N}$ being deterministic, $K$ a fixed positive integer, $\mathbf{F}=[f_{it}]\in{\mathbb{R}}^{K\times T}$ and $\mathbf{\Psi}=[\psi_{it}]\in{\mathbb{R}}^{N\times T}$ two random matrices. Here we assume that the entries of $\mathbf{F}$ (factors) are i.i.d. centered random variables with unit variance, and so are the entries of $\mathbf{\Psi}$ (noises). The entries of $\mathbf{F}$ are independent from the entries of $\mathbf{\Psi}$, but $f_{11}$ and $\psi_{11}$ are not necessarily identically distributed. All rows of $\bm{\Gamma}:=\begin{bmatrix}\mathbf{L}&\mathbf{\Lambda}\end{bmatrix}$ are nonzero.

Our fundamental result is Theorem 2.4: for every $K$-factor model, if the entries $\psi_{it}$ of the noise matrix $\mathbf{\Psi}$ have finite fourth moments and $\mathbf{\Lambda}$ is diagonal, or if there exists $\varepsilon>0$ such that $\operatorname{\mathbb{E}}\left(\,|\psi_{it}|^{4}(\log|\psi_{it}|)^{2+2\varepsilon}\right)<\infty$, then the eigenvalues of the sample correlation matrices $\mathbf{C}$ and those of the sample covariance matrices $\mathbf{S}$ with unit-variance data are asymptotically equal, which extends the result of El Karoui (2009).

In PCA, methods for estimating the number of factors in a sample have been studied by (Bai & Ng, 2002; Lam & Yao, 2012; Aït-Sahalia & Xiu, 2017) in econometrics, etc.

Meanwhile, rules for the retention of principal components have been proposed in many lectures (see, e.g., (Jackson, 1993)). Broken-stick (BS) rule (Frontier, 1976; Jackson, 1993) is a peculiar rule among such rules. BS rule compares the spectral distribution of $\mathbf{C}$, with a distribution of the mean lengths of the subintervals obtained by a random partition of [0,  1] (Holst, 1980). Specifically, for the number of factors or significant principal components, the BS rule provides $i-1$ where $i$ is the lowest among $[1,\,N]$ for which the $i$th largest eigenvalue of $\mathbf{C}$ does not exceed the sum $1/i+1/(i+1)+\cdots+1/N$. The BS rule depends on the number and growth rates of unbounded eigenvalues of $\mathbf{C}$. The idea of the BS rule, initially coming from the species occupation model in an ecological system, is not evident in relation to the distribution of eigenvalues of $\mathbf{C}$.

First, we study the asymptotic spectral properties of $\mathbf{C}$, through the fundamental theorem of $K$-factor model (Theorem 2.4) and techniques from the random matrix theory, for two illustrative $K$-factor models generalizing the equi-correlated normal population. One of the two is a $K$-factor model such that $\mathbf{\Lambda}=\sigma\mathbf{I}_{N}$ and the rows of the factor loading matrix $\mathbf{L}$ have a common length. We call this a constant length factor loading model (CLFM). The other is a $K$-factor model such that the rows of $\mathbf{L}$ and the diagonal entries of $\mathbf{\Lambda}$ are convergent. This model was introduced in Akama (2023), and is called an asymptotic convergent factor model (ACFM) here. We establish that under certain mild conditions, the BS rule precisely matches the principal rank of the factor loading matrix $\mathbf{L}$ for a CLFM (see Theorem 2.7), but the essential rank of $\mathbf{L}$ for an ACFM (see Theorem 2.8).

Then we calculate the BS rule and some modern factor number estimators such as the adjusted correlation thresholding (Fan et al., 2022) and Bai-Ng’s rule (Bai & Ng, 2002) based on an information criterion, for financial datasets and a biological dataset (Quadeer et al., 2014). The financial datasets obtained by Fan et al. (2022) from Fama-French 100 portfolio (Fama & French, 1993) by cleaning-up outliers. The biological dataset is a binary multiple sequence alignment (MSA) of HCV genotype 1a (prevalent in North America) publicly available from the Los Alamos National Laboratory database (Kuiken et al., 2004).

The structure of this paper is as follows. Section 2 formally introduces the model and presents the main theorems, with proofs provided in Section A of the supplementary material. Numerical analysis is performed on actual financial and biological datasets in Section 3. Section 4 is the conclusion. Section B of the supplementary material includes several important lemmas derived from general random matrix theories.

In this paper we consider the following model of data matrix:

We often write $\mathbf{X}$ in a compact form

$$\mathbf{X}=\mathbf{M}+\bm{\Gamma}\mathbf{Z}$$

with

$$\bm{\Gamma}=\begin{bmatrix}\mathbf{L}&\mathbf{\Lambda}\end{bmatrix}\quad\text{and}\quad\mathbf{Z}=\begin{bmatrix}\mathbf{F}\\ \mathbf{\Psi}\end{bmatrix}.$$

Let $\bm{x}_{j}$ be the $j$th column of $\mathbf{X}$. Then $\bm{x}_{1},\dots,\bm{x}_{T}$ are i.i.d. random vectors. We recall that the population covariance and correlation matrices are $[\operatorname{\mathbb{C}ov}(x_{i1},x_{k1})]_{N\times N}$ and $[\operatorname{\mathbb{C}orr}(x_{i1},x_{k1})]_{N\times N}$, respectively, where $x_{ij}$ is the $i$th component of the vector $\bm{x}_{j}$. Note that under Definition 2.1, the population covariance matrix is

$$\bm{\Sigma}:=\bm{\Gamma}\bm{\Gamma}^{\top}=\mathbf{L}\mathbf{L}^{\top}+\mathbf{\Lambda}\mathbf{\Lambda}^{\top},$$

and the population correlation matrix is

$$\mathbf{R}:=\bm{\Delta}^{-\frac{1}{2}}\bm{\Sigma}\bm{\Delta}^{-\frac{1}{2}},$$

where

$\displaystyle\bm{\Delta}=\operatorname{diag}(\operatorname{\mathbb{V}}(x_{11}),\dots,\operatorname{\mathbb{V}}(x_{N1}))$

is a diagonal matrix containing the variance of components of $\bm{x}_{1}$. Note that $\operatorname{\mathbb{V}}(x_{k1})$ is just the square of Euclidean norm of the $k$th row vector $\bm{\gamma}^{k}$ of $\bm{\Gamma}$, and is also the $k$th diagonal element of $\bm{\Sigma}$. So we can write $\bm{\Delta}=\operatorname{diag}(\bm{\Sigma})$, where for a square matrix $\mathbf{A}$, $\operatorname{diag}(\mathbf{A})$ denotes the diagonal matrix with the same diagonal as $\mathbf{A}$.

The population covariance matrix $\bm{\Sigma}$ and correlation matrix $\mathbf{R}$ play important roles in multivariate statistics. However they are not always available. In order to estimate them, we define the theoretically-centered sample covariance matrix ${\tilde{\mathbf{S}}}$ as

$\displaystyle{\tilde{\mathbf{S}}}$

$\displaystyle=\frac{1}{T}(\mathbf{X}-\mathbf{M})(\mathbf{X}-\mathbf{M})^{\top}.$

(1)

Let $\tilde{\mathbf{D}}=\operatorname{diag}({\tilde{\mathbf{S}}})$. Then the theoretically-centered sample correlation matrix $\tilde{\mathbf{C}}$ is defined as (see e.g. El Karoui (2009))

$\displaystyle\tilde{\mathbf{C}}$

$\displaystyle=\tilde{\mathbf{D}}^{-\frac{1}{2}}{\tilde{\mathbf{S}}}\tilde{\mathbf{D}}^{-\frac{1}{2}}.$

As the mean vector $\bm{\mu}$ is not always known either, we sometimes need to replace $\bm{\mu}$ by the sample mean

$\displaystyle\bar{\bm{x}}=\frac{1}{T}\sum_{t=1}^{T}\bm{x}_{t},$

and use it to define the data-centered sample covariance matrix

$\displaystyle\mathbf{S}=\frac{1}{T-1}(\mathbf{X}-\bar{\mathbf{X}})(\mathbf{X}-\bar{\mathbf{X}})^{\top},$

where $\bar{\mathbf{X}}=\begin{bmatrix}\bar{\bm{x}}&\cdots&\bar{\bm{x}}\end{bmatrix}_{N\times T}$, and the data-centered sample correlation matrix

$\displaystyle\mathbf{C}=\mathbf{D}^{-\frac{1}{2}}\mathbf{S}\mathbf{D}^{-\frac{1}{2}},$

where $\mathbf{D}=\operatorname{diag}(\mathbf{S})$.

In this paper we focus on the sample correlation matrices. Note that if one row of $\bm{\Gamma}$ is identically $0$, then the corresponding row of $\mathbf{X}$ is deterministically equal to the mean, and the correlation between this row and any other row is not defined. But it is easy to recognize such a row in the data and to eliminate it before further treatment. So we can assume

1. A1

   The rows of $\bm{\Gamma}$ are nonzero:

   $$\bm{\gamma}^{i}\neq\bm{0},\quad i=1,\dots,N.$$

We study the limiting locations of eigenvalues of sample correlation matrices $\mathbf{C}$ and $\tilde{\mathbf{C}}$ in the proportional limiting regime

$\displaystyle{N,T\to\infty},\ \frac{N}{T}\to c>0,$

which will be denoted as

$${N,T\to\infty}$$

for simplicity. In this regime, Theorem 1 in El Karoui (2009) related the spectral properties (limiting spectral distributions and limit locations of individual eigenvalues) of $\mathbf{C}$ and $\tilde{\mathbf{C}}$ to those of sample covariance matrices $\bm{\Delta}^{-\frac{1}{2}}\mathbf{S}\bm{\Delta}^{-\frac{1}{2}}$ and $\bm{\Delta}^{-\frac{1}{2}}{\tilde{\mathbf{S}}}\bm{\Delta}^{-\frac{1}{2}}$, in condition that $\left|\!\left|\!\left|\bm{\Delta}^{-\frac{1}{2}}\bm{\Gamma}\right|\!\right|\!\right|$, the spectral norm of $\bm{\Delta}^{-\frac{1}{2}}\bm{\Gamma}$, are uniformly bounded. However, due to the presence of $\mathbf{L}$, this condition is not satisfied by some factor models. For example, let us consider an ENP studied in (Fan & Jiang, 2019; Akama & Husnaqilati, 2022; Akama, 2023) and defined in the following Definition 2.2. Then $\left|\!\left|\!\left|\bm{\Delta}^{-\frac{1}{2}}\bm{\Gamma}\right|\!\right|\!\right|=\sqrt{(L^{2}N+\sigma^{2})/(L^{2}+\sigma^{2})}$ which diverges to $\infty$ as $N\to\infty$.

In this paper, for the $i$th largest singular value $s_{i}\left(\mathbf{A}\right)$ of a matrix $\mathbf{A}$, we prove the following Theorem in Section A:

Thanks to this theorem, we managed to generalize Theorem 1 in El Karoui (2009) to general $K$-factor models with possibly unbounded $\bm{\Delta}^{-\frac{1}{2}}\bm{\Gamma}$.

Before stating our first main result, we add some additional assumptions.

1. A2

   Each random variable $\psi_{it}$ has finite forth moment:

   $$\operatorname{\mathbb{E}}\left(|\psi_{it}|^{4}\right)<\infty.$$

2. A3

   One of the following assumptions holds: (a) There exists $\varepsilon>0$ such that

   $$\operatorname{\mathbb{E}}\left(\,|\psi_{it}|^{4}(\log|\psi_{it}|)^{2+2\varepsilon}\right)<\infty,$$

   or (b) $\mathbf{\Lambda}$ is diagonal.

For two nonnegative sequences $\{a_{n}\}$ and $\{b_{n}\}$, by $a_{n}\sim b_{n}$ we mean that there is a positive sequence $\{\tau_{n}\in(0,1)\}$ such that $\lim_{n\to\infty}\tau_{n}=1$, and for large enough $n$, $\tau_{n}a_{n}\leq b_{n}\leq\tau_{n}^{-1}a_{n}$. For a family of such sequences $\{a_{n}^{(i)}\}_{i}$ and $\{b_{n}^{(i)}\}_{i}$, we say that $a_{n}^{(i)}\sim b_{n}^{(i)}$ uniformly if there is $\{\tau_{n}\in(0,1)\}$ independent of $i$, with $\lim_{n\to\infty}\tau_{n}=1$, such that $\tau_{n}a_{n}^{(i)}\leq b_{n}^{(i)}\leq\tau_{n}^{-1}a_{n}^{(i)}$ holds for all $i$ and for large enough $n$.

It should be noted that no bounding condition for $\bm{\Delta}^{-1/2}\mathbf{\Lambda}$ is required. Therefore, when $\mathbf{L}=\bm{0}$, the theorem is applicable to a general linear model $\mathbf{X}=\mathbf{M}+\mathbf{\Lambda}\mathbf{\Psi}$, where $\mathbf{\Lambda}$ can be unbounded. Furthermore, assuming $\mathbf{L}\neq\bm{0}$, we accommodate different distributions for factor and noise components, for example, the factors are allowed to have a heavy-tailed distribution, along with a light-tailed noise.

The general result being stated, it is still complex to determine the asymptotic locations of a sample covariance matrix in general. We content ourselves with some particular cases that extend an ENP.

We will check whether the following real datasets are generated by a CLFM or an ACFM, based on Theorems 2.6 and 2.8:

1. 1.

   the datasets Fan et al. (2022) obtained by cleaning outliers from the datasets of the daily excess returns (Jensen, 1968) of Fama-French 100 portfolios (Fama & French (1993, 2015). See Prof. French’s data library http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/).

2. 2.

   binary multiple sequence alignment (Quadeer et al., 2014) by the courtesy of Prof. Quadeer.

For the correlation matrices $\mathbf{C}$ of these datasets, we will compute the eight quantities:

1. 1.

   the sample size $T$,

2. 2.

   $N/T$,

3. 3.

   $\lambda_{1}\left(\mathbf{C}\right)/N$ (cf. Theorem 2.6 (3)),

4. 4.

   the dimension $N$,

5. 5.

   the broken-stick rule $\operatorname{BS}(\mathbf{C})$,

6. 6.

   Adjusted correlation thresholding (Fan et al., 2022) $\operatorname{ACT}(\mathbf{C})$.

7. 7.

   Bai-Ng’s rule based on information criterion (Bai & Ng, 2002).

In a stock return dataset, $N$ ($T$, resp.) intends the number of companies (trading days$-1$, resp.). For $i,t$ $({1\leq i\leq N},{1\leq t\leq T})$, the $i$th row of a data matrix $\mathbf{X}$ is for the $i$th company and $t$th column is for the $t$th trading day. The factors are those of Fama-French, for instance.

We suppose the proportional limiting regime for factor models. We established a general theorem for the asymptotic equispectrality of $\mathbf{C}$ and the naturally normalized sample covariance matrix for factor models. Then, we derived the following assertions: for the introduced two models (CLFM and ACFM), the limiting largest bounded eigenvalue and the limiting spectral distribution of the sample correlation matrix $\mathbf{C}$, are scaling of those of $\mathbf{C}$ of the i.i.d. case, with the scaling constant $1-\rho$. Here $\rho$ is “a reduced correlation coefficient”. The largest eigenvalue of $\mathbf{C}$ divided by the order $N$ of $\mathbf{C}$ converges almost surely to $\rho$.

Notably, for an ACFM, the BS rule computes not the rank of the factor loading matrix $\mathbf{L}$, but the essential rank $(=1)$ (Theorem 2.8). In other words, the deterministic decreasing vector $\bm{y}=(\sum_{j=k}^{K}(jK)^{-1})_{k=1}^{K}$ of the BS rule examines an intrinsic structure of $\mathbf{L}$, for an ACFM. Moreover, $\bm{y}\in{\mathbb{R}}^{K}$ is the descending list of the limits of the eigenvalues of $\mathbf{L}\mathbf{L}^{\top}\in{\mathbb{R}}^{N\times N}$ in $N\to\infty$, for a random factor loading matrix $\mathbf{L}=\left[B_{ik}\sqrt{y_{ik}}\right]_{{1\leq i\leq N},{1\leq k\leq K}}$ where for each $i$ $({1\leq i\leq N})$,

• •

  $\bm{b}_{i}=(B_{ik})_{k=1}^{K}$ is a random vector of independent Rademacher variables,

• •

  $\bm{y}_{i}=(y_{ik})_{k=1}^{K}$ is such that $y_{ik}$ is the length of the $k$th longest subintervals obtained by $K-1$ independent uniformly distributed separators of $[0,\,1]$;

and $\bm{b}_{1},\ldots,\bm{b}_{N},\ \bm{y}_{1},\ldots,\bm{y}_{N}$ are independent. In this case, $\mathbf{L}$ does not satisfy the condition of an ACFM. A limit theory for the order statistics of the eigenvalues of a sample covariance/correlation matrix is awaited to analyze the BS rule.