Tests for Large-Dimensional Shape Matrices via Tyler’s M Estimators

Abstract

Tyler’s M estimator, as a robust alternative to the sample covariance matrix, has been widely applied in robust statistics. However, classical theory on Tyler’s M estimator is mainly developed in the low-dimensional regime for elliptical populations. It remains largely unknown when the parameter of dimension p grows proportionally to the sample size n for general populations. By using the eigenvalues of Tyler’s M estimator, this article develops tests for the identity and equality of shape matrices in a large-dimensional framework where the dimension-to-sample size ratio p/n has a limit in (0, 1). The proposed tests can be applied to a broad class of multivariate distributions including the family of elliptical distributions (see model (2.1) for details). To analyze both the null and alternative distributions of the proposed tests, we provide a unified theory on the spectrum of a large-dimensional Tyler’s M estimator when the underlying population is general. Simulation results demonstrate good performance and robustness of our tests. An empirical analysis of the Fama-French 49 industrial portfolios is carried out to demonstrate the shape of the portfolios varying. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

KEYWORDS:

• Central limit theorem
• High-dimensional tests
• Linear spectral statistics
• Shape matrix
• Tyler’s M estimator

1 Introduction

Tests for shape (or scatter) matrix for multivariate nonnormal data play the same role as tests for covariance matrix in multivariate normal distributions. This is one of the central parts in multivariate analysis (Anderson 1984). Many applications including, but not limited to, signal processing, wireless communication system, biology, and financial engineering may impose certain structures on shape or covariance matrix. Thus, test of shape or covariance matrix structure is of primary interest in modern multivariate analysis. Motivated by an empirical analysis of financial economic data in Section 3.3, this article aims to develop tests for shape matrix of large-dimensional data without imposing normality assumption.

As a natural extension of multivariate normal distribution, elliptical symmetric distributions (Fang, Kotz, and Ng 1990) are commonly used in the literature of nonparametric multivariate analysis. The shape or scatter matrix of elliptical distribution is proportional to the covariance matrix if its second moment exists and is finite. This motivates researchers to develop statistical inference procedures for the shape matrix rather than the covariance matrix of elliptical distributions. In many real-world applications such as the empirical analysis in Section 3.3, data may exhibit heavy-tailed nature or potentially contaminated by outliers. As it is well known, sample covariance matrix based inference on the shape is sensitive to outliers and heavy tails, and therefore it will lose its efficiency and become less reliable when the underlying populations are elliptically distributed. To cope with this issue, many robust alternatives are used in the literature, including the M-estimators (Huber 1964), the minimum volume ellipsoid and minimum determinant estimators (Rousseeuw 1985) and Tyler’s M estimator (Tyler 1987), etc. Hallin and Paindaeine (2006) and Hallin, Oja, and Paindaeine (2006) proposed a class of rank-based procedures for testing sphericity of elliptical distribution (equivalently testing the shape matrix being an identity matrix), and showed their test for shape is semiparametrically efficient. All aforementioned works assume that the data are finite and fixed-dimensional. The theoretical analyses in these works may not be applicable for large-dimensional data, which are common in many real-world applications nowadays. Thus, it is of great importance in developing new theory that is applicable for establishing asymptotic properties of statistical inference procedures for large-dimensional shape matrix.

Tests for large-dimensional covariance matrix are closely related to tests on shape matrix and have received increasing attention in the literature. Tests for the identity and sphericity of covariance matrix have been studied in the literature Ledoit and Wolf 2002; Bai et al. 2009; Chen, Zhang, and Zhong 2010. Tests for the equality and proportionality of covariance matrices in the two-sample problems have been studied in Schott (2007), Li and Chen (2012), and Zou et al. (2021). A common feature of these conventional tests is that their statistics are constructed from sample covariance matrix or sample correlation matrix, which are both well-behaved for those Gaussian-type and light-tailed populations, but are sensitive to outliers and heavy-tailed distributions. Meanwhile, random matrix theory for sample covariance or correlation matrices have been well developed for Gaussian-type and light-tailed populations, but is less developed for heavy-tailed populations. One of primary goals of this article is to develop new theory for large-dimensional shape matrix estimate with heavy-tailed populations, and then apply the new theory for developing tests on shape matrices.

In this article, we develop tests for shape matrix in a large-dimensional regime where the parameter of dimension p grows to infinity proportionally to the sample size n such that $p/n\to c\in (0,1)$. We focus on two fundamental problems in statistical inference on the shape matrix: test whether the shape matrix of a population is identity and test whether the shape matrices of two populations are equal. The proposed test statistics are based on the eigenvalues of Tyler’s M estimator since among all aforementioned robust estimators, Tyler’s M estimator enjoys certain superiorities: it is affine equivariant and, under ellipticity assumption, it is the “most robust” one that minimizes the maximum asymptotic variance in the finite-dimensional asymptotic framework (Tyler 1987; Magyar and Tyler 2014). The proposed procedures are applicable to a broad class of multivariate distributions defined in (2.1). In order to analyze both the null and alternative distributions of the proposed statistics, we develop new theory on the limiting spectrum of a large-dimensional Tyler’s M estimator. Specifically, under the population model (2.1), we establish both the convergence of the empirical spectral distribution (ESD) and the central limit theorem (CLT) for linear spectral statistics (LSSs) of a large-dimensional Tyler’s M estimator. Such theoretical results are brand new and are of independent interest in random matrix theory, with potential applications in developing new approaches for other related robust inferences in large dimensions. Let’s briefly summarize the major contributions of this work below.

From methodological point of view, the proposed tests have several appealing properties. First, they are all distribution-free. We show that their limiting null distributions depend only on the dimension p and sample size n of a dataset. Second, Tyler’s M estimator is an affine equivariant estimator of the shape matrix in the sense that when we apply the affine transformation $\mathbf{x}\mapsto \mathbf{\text{Ax}}\triangleq \mathbf{z}$ to the data, their corresponding Tyler’s M estimators, denoted as ${\mathbf{C}}_n(\mathbf{x})$ and ${\mathbf{C}}_n(\mathbf{z})$, satisfy ${\mathbf{C}}_n(\mathbf{z})=\mathbf{A}{\mathbf{C}}_n(\mathbf{x})\mathbf{A}\prime$ for any deterministic full-rank square matrix $\mathbf{A}$. Thus, Tyler’s M estimator behaves similarly to the sample covariance matrix under linear transformations. In particular, for tests of the two-sample shape matrix problem, this affine equivariance property of Tyler’s M estimator can eliminate the impact from the unknown shape matrix under the null hypothesis. This ensures that the proposed tests are universal to the underlying populations. At the same time, it allows us to investigate some nontrivial cases, including those where the shape matrix has an unbounded spectral norm, etc. These findings are also conformed by both theoretical analysis and Monte Carlo studies.

From technical aspects, we develop new strategies and techniques to study the limiting spectrum of a large-dimensional Tyler’s M estimator from general populations. As far as we know, existing theories on Tyler’s M estimator are mainly developed in the low-dimensional regime for elliptical populations, but remain largely unknown when the dimension p grows to infinity and the underlying population deviates away from the elliptical ones. In addition, Tyler’s M estimator is defined to be the solution to a fixed point equation and has no explicit expression. Therefore, standard approaches and tools in random matrix theory given in Bai and Silverstein (2010) cannot be directly imitated here. To solve this problem, we propose a new theoretical development strategy by comparing the resolvent of Tyler’s M estimator to that of sample spatial-sign covariance matrix (Locantore et al. 1999) for the first time, which is eventually achieved by employing a series of new techniques, mainly the method of moments. Note that our techniques can be extended similarly to the study of other implicitly defined M-estimators, such as the regularized Tyler’s M estimator (Kammoun et al. 2016), where $\mathrm{\text{lim}}p/n$ is then allowed to be larger than 1. However, this is beyond the scope of this article.

The rest of the article is organized as follows. Section 2 introduces our model assumptions, and develops tests for shape matrices in both the one-sample and two-sample cases. In Section 3, we conduct simulation results, and illustrate the proposed methodology via an empirical analysis of a real dataset. Section 4 is dedicated to the development of new theory regarding the limiting spectrum of a large-dimensional Tyler’s M-estimator. These are the core technical tools used to derive both the null and alternative distributions of our test statistics. All the proofs are provided in the supplementary material of this article.

2 New Test Procedures for Large-Dimensional Shape Matrix

2.1 Model of the Population

Consider a random vector $\mathbf{x}\in {\mathbb{R}}^p$ admitting the following stochastic representation, (2.1) $$\mathbf{x}\stackrel{d}{=}w{\mathbf{\Sigma }}^{\frac{1}{2}}\mathbf{z},$$(2.1) where w > 0 is a scalar random variable, $\mathbf{\Sigma }$ is a p × p deterministic positive definite matrix normalized by $\text{tr}(\mathbf{\Sigma })=p$ (normalization of $\mathbf{\Sigma }$ guarantees the model identifiability between w and $\mathbf{\Sigma }$), referred to as the shape matrix or scatter matrix of the population, $\mathbf{z}=(z_1,\dots ,z_p)\prime \in {\mathbb{R}}^p$ is a random vector with iid continuous coordinates $(z_i)$ having bounded densities and satisfying (2.2) $$\mathbb{E}(z_i)=0,\hspace{1em}\mathbb{E}(z_i^2)=1,\hspace{1em}\mathbb{E}(z_i^4)\triangleq \tau ,\hspace{1em}\mathbb{E}(z_i^6)<\infty .$$(2.2)

Remark 2.1.

It is worth noting that w and $\mathbf{z}$ in (2.1) are not necessarily independent. Therefore, our model is quite general and covers several important classes of multivariate distributions.

1. When the random variable w degenerates to a constant, model (2.1) becomes the well studied independent components model (ICM), which provides the most standard setup for independent component analysis (ICA). See Ilmonen and Paindaeine (2011) for detailed discussion about ICM and ICA. The ICM is also a conventional structure assumption in the study of random matrix theory (Bai and Silverstein 2010). As demonstrated by the real data example in Section 3.3, model (2.1) can be more general than ICM in practice.

2. Model (2.1) includes the elliptical distribution family with zero location vector as a special case. To see this, since if $\mathbf{x}$ follows an elliptical distribution with zero location vector, it has the stochastic representation $\mathbf{x}\stackrel{d}{=}v{\mathbf{\Sigma }}^{\frac{1}{2}}\mathbf{u},$ where v > 0 is a scalar random variable and u is a random vector uniformly distributed on the unit sphere in ${\mathbb{R}}^p$ (Fang, Kotz, and Ng 1990). Further, let $\mathbf{u}\stackrel{d}{=}\mathbf{z}/\left|\right|\mathbf{z}\left|\right|$ and $w=v/\left|\right|\mathbf{z}\left|\right|$, where $\mathbf{z}$ is set to be a standard Gaussian random vector (τ = 3). Thus, we have $\mathbf{x}\stackrel{d}{=}v{\mathbf{\Sigma }}^{\frac{1}{2}}\mathbf{u}\stackrel{d}{=}w{\mathbf{\Sigma }}^{\frac{1}{2}}\mathbf{z},$ which is indeed a special case of our model (2.1).

3. By using the same argument as that in (b), model (2.1) includes multivariate ${\mathcal{l}}_1$-norm symmetric distributions (Fang, Kotz, and Ng 1990) as special cases by taking all elements of $\mathbf{z}$ be iid according to the exponential distribution with mean 1.

4. With a similar argument to that in (b), model (2.1) includes a subclass of multivariate Liouville distributions (Fang, Kotz, and Ng 1990) as special cases by taking all elements of $\mathbf{z}$ be iid according to a Gamma distribution.

Detailed definitions and properties of elliptical distributions, ${\mathcal{l}}_1$-norm symmetric distributions and multivariate Liouville distributions are referred to Fang, Kotz, and Ng (1990).

2.2 Testing the Identity of a Shape Matrix

Let ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\in {\mathbb{R}}^p$ be a sequence of iid observations from the population $\mathbf{x}$ defined in (2.1). Consider testing the following hypothesis for the shape matrix $\mathbf{\Sigma }$, (2.3) $$H_0:\mathbf{\Sigma }={\mathbf{I}}_p\hspace{1em}\text{versus}\hspace{1em}{H}_1:\mathbf{\Sigma }\ne {\mathbf{I}}_p$$(2.3) in the large-dimensional asymptotic framework where $p\to \infty ,n\to \infty ,c_n\triangleq p/n\to c\in (0,1).$

Our procedures are based on the eigenvalues of Tyler’s M estimator (Tyler 1987) of the shape $\mathbf{\Sigma }$, which is defined to be the matrix ${\mathbf{C}}_n$ satisfying (2.4) $${\mathbf{C}}_n=\frac{p}{n}\sum _{j=1}^n\frac{{\mathbf{x}}_j{\mathbf{x}}_j^{\prime }}{{\mathbf{x}}_j^{\prime }{\mathbf{C}}_n^{-1}{\mathbf{x}}_j},\hspace{1em}\text{tr }{\mathbf{C}}_n=p.$$(2.4)

Tyler’s M estimator is affine equivariant and under ellipticity assumption and compared with existing robust estimators, it is the “most robust” one that minimizes the maximum asymptotic variance in the finite-dimensional asymptotic framework (Tyler 1987; Magyar and Tyler 2014). Kent and Tyler (1988) showed that if any linear subspace in ${\mathbb{R}}^p$ of dimension d ($1\le d\le p-1$) contains less than nd/p observations, then there exists a unique solution to the (2.4). In particular, this condition holds with probability one for the random sample $\{{\mathbf{x}}_j,j=1,\dots ,n\}$ from model (2.1).

We consider the following two statistics for testing (2.3), $$T_{1,\mathit{\text{tr}}}=\text{tr}({\mathbf{C}}_n^2)\hspace{1em}\text{and}\hspace{1em}{T}_{1,\mathrm{\text{log}}}=\hspace{0.17em}\mathrm{\text{log det}}({\mathbf{C}}_n).$$

The trace-based statistic $T_{1,\mathit{\text{tr}}}$ is motivated by the tests for covariance matrix from John (1971, 1972), which measures the squared Frobenius distance between the Tyler’s M estimator ${\mathbf{C}}_n$ and the shape matrix $\mathbf{\Sigma }$ under the null. And the log-based statistic $T_{1,\mathrm{\text{log}}}$ is an analogue of likelihood ratio statistic (Anderson 1984), which is based on the entropy loss $L({\mathbf{C}}_n,\mathbf{\Sigma })\triangleq \text{tr}{\mathbf{C}}_n{\mathbf{\Sigma }}^{-1}-\hspace{0.17em}\mathrm{\text{log det}}({\mathbf{C}}_n{\mathbf{\Sigma }}^{-1})-p$ between ${\mathbf{C}}_n$ and $\mathbf{\Sigma }$ under the null (James and Stein 1961). In the definition of these test statistics, we use the normalization $\text{tr}{\mathbf{C}}_n=p$, see (2.4). The null distributions of the two statistics are presented in the following theorem.

Theorem 2.1.

In the large-dimensional asymptotic framework, we have under H₀, $$\begin{array}{c}\frac{T_{1,\mathit{\text{tr}}}-p(1+c_n)+c(c-3)/(1-c)}{2c}\stackrel{D}{\to }N(0,1),\\ \frac{T_{1,\mathrm{\text{log}}}+p-(p-n)\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1-c_n)-2^{-1}\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1-c)+c/(1-c)}{\sqrt{-2\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1-c)-2c}}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\stackrel{D}{\to }N(0,1).\end{array}$$

We prove this theorem by using the new theory on the CLT for LSSs of Tyler’s M estimator, see Theorem 4.2 in Section 4. In particular, the limiting null distributions of the two statistics are only relevant to the population dimension p and sample size n, independent of the fourth moment τ of $(z_i)$ and also the true underlying distribution of the population $\mathbf{x}$.

To analyze the power of the proposed tests, we consider a local alternative such that $\mathbf{\Sigma }$ is diagonal or τ = 3, and compare our Tyler’s M based test $T_{1,\mathit{\text{tr}}}$ with John’s test (John 1972), defined as ${\tilde{T}}_{1,\mathit{\text{tr}}}\triangleq \text{tr}({\mathbf{S}}_n^2)/{[\text{tr}({\mathbf{S}}_n)/p]}^2,$ where ${\mathbf{S}}_n=n^{-1}{\sum }_{j=1}^n{\mathbf{x}}_j{\mathbf{x}}_j^{\prime }=\frac{1}{n}{\sum }_{j=1}^n{w}_j^2{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j{\mathbf{z}}_j^{\prime }{\mathbf{\Sigma }}^{\frac{1}{2}}$ is the sample covariance matrix. Under the local alternative, the asymptotical distributions of both $T_{1,\mathit{\text{tr}}}$ and ${\tilde{T}}_{1,\mathit{\text{tr}}}$ are derived in Theorems 2.2 and 2.3, respectively. Using these two theorems, we further derive the asymptotic power functions of these two tests in Proposition 2.1.

Theorem 2.2.

Suppose that $0<\mathrm{\text{lim}}\underset{p}{\mathrm{\text{inf}}}{\lambda }_{\mathrm{\text{min}}}(\mathbf{\Sigma })\le \mathrm{\text{lim}}\underset{p}{\mathrm{\text{sup}}}{\lambda }_{\mathrm{\text{max}}}(\mathbf{\Sigma })<\infty$, and as $p\to \infty$, ${\beta }_{k,p}\triangleq p^{-1}\text{tr}{\mathbf{\Sigma }}^k\to {\beta }_k,k=2,3,4.$ Then, under either of the two conditions: (i) $\mathbf{\Sigma }$ is diagonal or (ii) τ = 3, we have $T_{1,\mathit{\text{tr}}}-\text{tr}({\mathbf{\Sigma }}^2)-p{c}_n\stackrel{D}{\to }N(\mu ,{\sigma }^2),$ where the limiting mean and limiting variance are $\mu =2c(c-2){(c-1)}^{-1}-c{\beta }_2+c(\tau -1)({\beta }_2-1),$ ${\sigma }^2=4c^2{\beta }_2^2+4c(\tau -1)({\beta }_2^3-2{\beta }_2{\beta }_3+{\beta }_4).$

Theorem 2.3.

In addition to the assumptions of Theorem 2.2, suppose that w is independent of $\mathbf{z}$ and has finite eighth moments, that is, ${\beta }_{\mathit{\text{wk}}}\triangleq \mathbb{E}(w^k)<\infty ,k\in \{1,\dots ,8\}.$ Then, we have (2.5) $$\begin{array}{c}\frac{{\tilde{T}}_{1,\mathit{\text{tr}}}-\text{tr}{\mathbf{\Sigma }}^2-c_{n}p{\beta }_{w2}^{-2}{\beta }_{w4}-n^{-1}\text{tr}{\mathbf{\Sigma }}^2\left[(\tau -1){\beta }_{w2}^{-2}{\beta }_{w4}-1\right]}{{\sigma }_{\mathit{\text{nJ}}}}\\ \hspace{1em}\stackrel{D}{\to }N(0,1),\end{array}$$(2.5) where the scaling parameter σ_nJ is defined by $${\sigma }_{\mathit{\text{nJ}}}^2=\{\begin{array}{c}n{c}_n^4\left[4{\beta }_{w4}({\beta }_{w4}^2-{\beta }_{w2}{\beta }_{w6}){\beta }_{w2}^{-2}+{\beta }_{w8}-{\beta }_{w4}^2\right]{\beta }_{w2}^{-4},\hfill \\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{var}(w^2)>0,\hfill \\ 4c^2{\beta }_2^2+4c(\tau -1)({\beta }_2^3-2{\beta }_2{\beta }_3+{\beta }_4),\hfill \\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{var}(w^2)=0.\hfill \end{array}$$

Proposition 2.1.

Suppose that the assumptions of Theorem 2.3 hold. Then, the asymptotic power functions of $T_{1,\mathit{\text{tr}}}$ and ${\tilde{T}}_{1,\mathit{\text{tr}}}$ are given by (2.6) $${\beta }_{T_{1,\mathit{\text{tr}}}}=\Phi \left(\frac{-2c{z}_{\alpha }+\left(\text{tr}{\mathbf{\Sigma }}^2-p\right)\left[1+n^{-1}(\tau -2)\right]}{\sqrt{4c^2{\beta }_2^2+4c(\tau -1)({\beta }_2^3-2{\beta }_2{\beta }_3+{\beta }_4)}}\right),$$(2.6) (2.7) $${\beta }_{{\tilde{T}}_{1,\mathit{\text{tr}}}}=\Phi \left(-z_{\alpha }\frac{{\sigma }_{\mathit{\text{nJ}}0}}{{\sigma }_{\mathit{\text{nJ}}}}+\frac{\left(\text{tr}{\mathbf{\Sigma }}^2-p\right)\left[1+n^{-1}(\tau -1){\beta }_{w2}^{-2}{\beta }_{w4}-n^{-1}\right]}{{\sigma }_{\mathit{\text{nJ}}}}\right),$$(2.7) where $\Phi$ is the standard normal distribution function, $z_{\alpha }={\Phi }^{-1}(1-\alpha )$ is the upper α quantile and (2.8) $${\sigma }_{\mathit{\text{nJ}}0}^2=\{\begin{array}{c}n{c}_n^4\left[4{\beta }_{w4}({\beta }_{w4}^2-{\beta }_{w2}{\beta }_{w6}){\beta }_{w2}^{-2}+{\beta }_{w8}-{\beta }_{w4}^2\right]{\beta }_{w2}^{-4},\hfill \\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{var}(w^2)>0,\hfill \\ 4c^2,\hfill \\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{var}(w^2)=0.\hfill \end{array}$$(2.8)

Moreover, the asymptotic relative efficiency (ARE) of $T_{1,\mathit{\text{tr}}}$ with respect to ${\tilde{T}}_{1,\mathit{\text{tr}}}$ is $$\begin{array}{c}\text{ARE}(T_{1,\mathit{\text{tr}}},{\tilde{T}}_{1,\mathit{\text{tr}}})\doteq \frac{{\sigma }_{\mathit{\text{nJ}}}}{\sqrt{4c^2{\beta }_2^2+4c(\tau -1)({\beta }_2^3-2{\beta }_2{\beta }_3+{\beta }_4)}}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\to \{\begin{array}{cc}\infty ,\hfill & \text{var}(w^2)>0,\hfill \\ 1,\hfill & \text{var}(w^2)=0.\hfill \end{array}\end{array}$$

Remark 2.2.

Theorem 2.3 unveils a phase transition phenomenon for ${\tilde{T}}_{1,\mathit{\text{tr}}}$ such that its convergence rate is of the order $O\left(\sqrt{n}\right)$ for the case of nondegenerate w² such as scale mixture models and multivariate t-distributions, while is of the order O(1) for the degenerate case where w² is a constant, say, the independent components model.

Remark 2.3.

For the independent components model, Proposition 2.1 demonstrates an equivalent role of Tyler’s M estimator to the sample covariance matrix in terms of the asymptotic power functions (the two are exactly the same). However, in contrast to the test $T_{1,\mathit{\text{tr}}}$, we note that by setting $w\equiv \sigma$ and $\mathbf{\Sigma }={\mathbf{I}}_p$ in (2.5), the null distribution of the sample covariance matrix based test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ still depends on the unknown fourth moment parameter τ of the atom random variable z_i, which is not distribution free. Thus, an extra estimation of this parameter should be carried out for the test procedure.

Remark 2.4.

For populations with nondegenerate w², Proposition 2.1 reveals an overwhelming superiority of Tyler’s M estimator to the sample covariance matrix. In particular, for those local alternatives where $\text{tr}{\mathbf{\Sigma }}^2-p=o\left(\sqrt{p}\right)$, Proposition 2.1 shows that the power of the Tyler’s M based test $T_{1,\mathit{\text{tr}}}$ tends to 1 while the sample covariance matrix based test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ has no power. An intuitive explanation for this is that the noise generated by the variable w is stronger than the signal $\text{tr}{\mathbf{\Sigma }}^2-p$. Therefore, removing the irrelevant variable w from the shape estimation is the main achievement of using Tyler’s M estimator.

2.3 Testing the Equality of Two Shape Matrices

Let ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ be two p-dimensional populations with stochastic representations ${\mathbf{x}}_1\stackrel{d}{=}{w}_1{\mathbf{\Sigma }}_1^{\frac{1}{2}}{\mathbf{z}}_1$ and ${\mathbf{x}}_2\stackrel{d}{=}{w}_2{\mathbf{\Sigma }}_2^{\frac{1}{2}}{\mathbf{z}}_2,$ respectively, where both ${\mathbf{z}}_1=(z_{1j})$ and ${\mathbf{z}}_2=(z_{2j})$ satisfy the moment conditions in (2.2) with $\mathbb{E}(z_{\mathit{\text{ij}}}^4)={\tau }_i$ for i = 1, 2. In this section, we consider the following hypothesis for the shape matrices ${\mathbf{\Sigma }}_1$ and ${\mathbf{\Sigma }}_2$, (2.9) $$H_0:{\mathbf{\Sigma }}_1={\mathbf{\Sigma }}_2=\mathbf{\Sigma }\hspace{1em}\text{versus}\hspace{1em}{H}_1:{\mathbf{\Sigma }}_1\ne {\mathbf{\Sigma }}_2,$$(2.9) where $\mathbf{\Sigma }$ is unknown, normalized by $\text{tr}(\mathbf{\Sigma })=p$. Assume ${\mathbf{\Sigma }}_1$ is invertible and denote by D_p the squared Frobenius distance between the matrix ${\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2/[\text{tr}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)/p]$ and ${\mathbf{I}}_p$, (2.10) $$D_p\triangleq \frac{1}{p}\text{tr}{\left[\frac{{\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2}{\text{tr}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)/p}-{\mathbf{I}}_p\right]}^2=\frac{p\text{tr}{({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)}^2}{\stackrel{2}{\text{tr}}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)}-1,$$(2.10) which equals zero if and only if H₀ is true.

Let $\{{\mathbf{x}}_{1j},j=1,\dots ,n_1\}$ and $\{{\mathbf{x}}_{2j},j=1,\dots ,n_2\}$ be two independent sets of observations from the two populations ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$, respectively. Their sample sizes n₁, n₂ and population dimension p are assumed to be large such that $p\to \infty ,n_1\to \infty ,n_2\to \infty ,c_{n_1}\triangleq p/n_1\to c_1\in (0,1),c_{n_2}\triangleq p/n_2\to c_2\in (0,1)$. Denote by ${\mathbf{M}}_{1n_1}$ and ${\mathbf{M}}_{2n_2}$ Tyler’s M estimators of the two shapes. That is, they are the unique solutions to the fixed point equations ${\mathbf{M}}_{i{n}_i}=p{n}_i^{-1}{\sum }_{j=1}^{n_i}{(\mathbf{x}{\prime }_{\mathit{\text{ij}}}{\mathbf{M}}_{i{n}_i}^{-1}{\mathbf{x}}_{\mathit{\text{ij}}})}^{-1}{\mathbf{x}}_{\mathit{\text{ij}}}\mathbf{x}{\prime }_{\mathit{\text{ij}}}$, $\text{tr}{\mathbf{M}}_{i{n}_i}=p$, i = 1, 2. Then, plugging ${\mathbf{M}}_{1n_1}$ and ${\mathbf{M}}_{2n_2}$ into (2.10) yields a naive estimate of D_p, (2.11) $${\widehat{D}}_p\triangleq \frac{p\text{tr}{({\mathbf{M}}_{1n_1}^{-1}{\mathbf{M}}_{2n_2})}^2}{\stackrel{2}{\text{tr}}({\mathbf{M}}_{1n_1}^{-1}{\mathbf{M}}_{2n_2})}-1.$$(2.11)

However, in the large-dimensional asymptotic framework, such ${\widehat{D}}_p$ is not a consistent estimator of D_p as demonstrated in the first conclusion of Theorem 2.5. It has a bias $c_{n_1}/(1-c_{n_1})+c_{n_2}$ which is not negligible asymptotically. By subtracting this bias from ${\widehat{D}}_p$, we finally define the following consistent estimator of D_p for testing (2.9), $$\begin{array}{c}{T}_{2,\mathit{\text{tr}}}\triangleq {\widehat{D}}_p-\frac{c_{n_1}}{1-c_{n_1}}-c_{n_2}\\ =\frac{p\text{tr}{({\mathbf{M}}_{1n_1}^{-1}{\mathbf{M}}_{2n_2})}^2}{\stackrel{2}{\text{tr}}({\mathbf{M}}_{1n_1}^{-1}{\mathbf{M}}_{2n_2})}-\frac{1}{1-c_{n_1}}-c_{n_2},\end{array}$$and a large value of such an estimate constitutes an evidence in favor of H₁. The limiting null distribution of $T_{2,\mathit{\text{tr}}}$ is given in the following theorem, which is also derived from our general theory presented in Section 4.

Theorem 2.4.

In the large-dimensional asymptotic framework, we have under H₀, $p{T}_{2,\mathit{\text{tr}}}\stackrel{D}{\to }N(\mu ,{\sigma }^2),$ where the limiting mean and limiting variance are $$\begin{array}{c}\mu =\frac{-3c_1+c_1^2-3c_2+8c_1{c}_2-3c_1^2{c}_2+c_2^2-c_1{c}_2^2}{(c_2-1){(c_1-1)}^2},\\ {\sigma }^2=\frac{4c_1^2+8c_1^3+8c_1{c}_2-8c_1^3{c}_2+4c_2^2-8c_1{c}_2^2+4c_1^2{c}_2^2}{{(c_1-1)}^4}.\end{array}$$

Remark 2.5.

The matrix ${\mathbf{M}}_{1n_1}^{-1}{\mathbf{M}}_{2n_2}$ in the test statistic $T_{2,\mathit{\text{tr}}}$ is the F-ratio of the two Tyler’s M estimators. By their affine equivariance, this test can eliminate the impact from the unknown shape matrix $\mathbf{\Sigma }$ under the null hypothesis. To see this, we denote ${\mathbf{C}}_{1n_1}$ and ${\mathbf{C}}_{2n_2}$ as two Tyler’s M estimators from the atom random vectors $\{{\mathbf{z}}_{1j},j=1,\dots ,n_1\}$ and $\{{\mathbf{z}}_{2j},j=1,\dots ,n_2\}$, respectively. That is, they satisfy ${\mathbf{C}}_{i{n}_i}=p{n}_i^{-1}{\sum }_{j=1}^{n_i}{(\mathbf{z}{\prime }_{\mathit{\text{ij}}}{\mathbf{C}}_{i{n}_i}^{-1}{\mathbf{z}}_{\mathit{\text{ij}}})}^{-1}{\mathbf{z}}_{\mathit{\text{ij}}}\mathbf{z}{\prime }_{\mathit{\text{ij}}}$, $\text{tr}{\mathbf{C}}_{i{n}_i}=p$, i = 1, 2. From the affine equivariance, we have the relation between ${\mathbf{M}}_{i{n}_i}$ and ${\mathbf{C}}_{i{n}_i}$ as (2.12) $${\mathbf{M}}_{i{n}_i}=\frac{{\mathbf{\Sigma }}_i^{\frac{1}{2}}{\mathbf{C}}_{i{n}_i}{\mathbf{\Sigma }}_i^{\frac{1}{2}}}{\text{tr}\left({\mathbf{\Sigma }}_i^{\frac{1}{2}}{\mathbf{C}}_{i{n}_i}{\mathbf{\Sigma }}_i^{\frac{1}{2}}\right)/p},\hspace{1em}i=1,2.$$(2.12)

Note that in (2.12), we normalize the matrix ${\mathbf{\Sigma }}_i^{\frac{1}{2}}{\mathbf{C}}_{i{n}_i}{\mathbf{\Sigma }}_i^{\frac{1}{2}}$ so that the trace of the term on the right hand side equals p, which meets the condition $\text{tr}{\mathbf{M}}_{i{n}_i}=p$. Then, under $H_0:{\mathbf{\Sigma }}_1={\mathbf{\Sigma }}_2=\mathbf{\Sigma }$, we have from (2.11) and (2.12), $$\begin{array}{c}{\widehat{D}}_p=\frac{p\text{tr}{\left({\mathbf{\Sigma }}_1^{-\frac{1}{2}}{\mathbf{C}}_{1n_1}^{-1}{\mathbf{\Sigma }}_1^{-\frac{1}{2}}{\mathbf{\Sigma }}_2^{\frac{1}{2}}{\mathbf{C}}_{2n_2}{\mathbf{\Sigma }}_2^{\frac{1}{2}}\right)}^2}{\stackrel{2}{\text{tr}}\left({\mathbf{\Sigma }}_1^{-\frac{1}{2}}{\mathbf{C}}_{1n_1}^{-1}{\mathbf{\Sigma }}_1^{-\frac{1}{2}}{\mathbf{\Sigma }}_2^{\frac{1}{2}}{\mathbf{C}}_{2n_2}{\mathbf{\Sigma }}_2^{\frac{1}{2}}\right)}-1\\ =\frac{p\text{tr}{({\mathbf{C}}_{1n_1}^{-1}{\mathbf{C}}_{2n_2})}^2}{\stackrel{2}{\text{tr}}({\mathbf{C}}_{1n_1}^{-1}{\mathbf{C}}_{2n_2})}-1,\end{array}$$which is indeed independent of the common shape matrix $\mathbf{\Sigma }$ under H₀. Therefore, the conclusion of Theorem 2.4 holds true only if the shape matrix $\mathbf{\Sigma }$ is invertible and no other constraints have to be imposed. This point allows us to investigate some nontrivial cases, for example, the one that $\mathbf{\Sigma }$ has an unbounded spectral norm, etc. This case is also studied by our simulations.

Next, we present the convergence of $T_{2,\mathit{\text{tr}}}$ for general ${\mathbf{\Sigma }}_1$ and ${\mathbf{\Sigma }}_2$, from which its power performance can be evaluated theoretically.

Theorem 2.5.

Suppose that $0<\mathrm{\text{lim}}\underset{p}{\mathrm{\text{inf}}}{\lambda }_{\mathrm{\text{min}}}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)\le \mathrm{\text{lim}}\underset{p}{\mathrm{\text{sup}}}{\lambda }_{\mathrm{\text{max}}}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)<\infty$ and as $p\to \infty$, (2.13) $$\frac{1}{p}\text{tr}{\left(\frac{{\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2}{\text{tr}({\mathbf{\Sigma }}_1^{-1}{\mathbf{\Sigma }}_2)/p}\right)}^k\to {\gamma }_k,\hspace{1em}k=2,3,4.$$(2.13)

1. The statistic $T_{2,\mathit{\text{tr}}}$ is a strongly consistent estimator of D_p, that is, $T_{2,\mathit{\text{tr}}}-D_p\to 0$, almost surely.

2. If the moment conditions in (2.2) hold for the two populations with ${\tau }_1={\tau }_2=3$, then we have $$p{T}_{2,\mathit{\text{tr}}}-p{D}_p\stackrel{D}{\to }N\left(\mu ,{\sigma }^2\right).$$

The limiting mean and limiting variance are $$\begin{array}{c}\mu =\frac{2c_1}{{(c_1-1)}^2}-\frac{c_1{c}_2+c_1{\gamma }_2}{c_1-1}-\frac{2c_2}{c_2-1}+c_2{\gamma }_2,\\ {\sigma }^2=4c_2(\frac{c_1^2(6{\gamma }_2+c_2)}{{(c_1-1)}^2}-\frac{2c_1\left({\gamma }_2{c}_2-{\gamma }_2^2+2{\gamma }_3\right)}{c_1-1}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}+2{\gamma }_4+{\gamma }_2(2{\gamma }_2^2+c_2{\gamma }_2-4{\gamma }_3)-\frac{4c_1^3}{{(c_1-1)}^3})\\ \hspace{1em}+\frac{12c_1^4}{{(c_1-1)}^4}-\frac{16c_1^3{\gamma }_2}{{(c_1-1)}^3}+\frac{4c_1^2\left(4{\gamma }_3-3{\gamma }_2^2\right)}{{(c_1-1)}^2}\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}-\frac{8c_1\left({\gamma }_2^3-2{\gamma }_3{\gamma }_2+{\gamma }_4\right)}{c_1-1}.\end{array}$$

Remark 2.6.

The first conclusion of Theorem 2.5 demonstrates the consistency of $T_{2,\mathit{\text{tr}}}$ toward the distance D_p and the second conclusion provides the non-null distribution of $T_{2,\mathit{\text{tr}}}$ for a local alternative with ${\tau }_1={\tau }_2=3$. In particular, the test can achieve full power if the quantity pD_p grows to infinity.

2.4 Adjustments When the True Location is Unknown

In this section, we discuss adjustments to the tests developed in Sections 2.2 and 2.3 when the true location vector $\mathit{\mu }$ of the population $\mathbf{x}$ is unknown, that is, (2.14) $$\mathbf{x}\stackrel{d}{=}\mathit{\mu }+w{\mathbf{\Sigma }}^{\frac{1}{2}}\mathbf{z},\hspace{1em}\text{tr}(\mathbf{\Sigma })=p.$$(2.14)

In addition to the assumptions of the model (2.1), we further assume that the random variable w is independent of $\mathbf{z}$ and they satisfy (2.15) $$\mathbb{E}(w^{-r})<\infty ,\hspace{1em}\mathbb{E}(w^r)<\infty ,\hspace{1em}\mathbb{E}(z_i^{\gamma })<\infty$$(2.15) for some r > 4 and $\gamma >8$.

For the one-sample case, given a sequence of iid observations ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$ following the same distribution as $\mathbf{x}$ in (2.14), Tyler’s M estimator, denoted by ${\overline{\mathbf{C}}}_n$, is defined to be the solution to the equation ${\overline{\mathbf{C}}}_n=p{n}^{-1}{\sum }_{j=1}^n{(({\mathbf{x}}_j-\widehat{\mathit{\mu }})\prime {\overline{\mathbf{C}}}_n^{-1}({\mathbf{x}}_j-\widehat{\mathit{\mu }}))}^{-1}({\mathbf{x}}_j-\widehat{\mathit{\mu }})({\mathbf{x}}_j-\widehat{\mathit{\mu }})\prime ,\text{ tr}({\overline{\mathbf{C}}}_n)=p$, where $\widehat{\mathit{\mu }}$ is an estimate of the true location $\mathit{\mu }$. For robustness, we take into account the location estimate $\widehat{\mathit{\mu }}$ of the form $\widehat{\mathit{\mu }}\triangleq {\sum }_{j=1}^n{\alpha }_j{\mathbf{x}}_j=\mathit{\mu }+{\sum }_{j=1}^n{\alpha }_j{w}_j{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j$, where ${\alpha }_j={({\sum }_{j=1}^n\left|\right|{\mathbf{x}}_j-\overline{\mathbf{x}}|{|}^{-2})}^{-1}\left|\right|{\mathbf{x}}_j-\overline{\mathbf{x}}|{|}^{-2}$ and $\overline{\mathbf{x}}=n^{-1}{\sum }_{i=1}^n{\mathbf{x}}_i$ is the sample mean. Thus, the two test statistics for the hypothesis (2.3) are $${\overline{T}}_{1,\mathit{\text{tr}}}=\text{tr}({\overline{\mathbf{C}}}_n^2)\hspace{1em}\text{and}\hspace{1em}{\overline{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}=\hspace{0.17em}\mathrm{\text{log det}}({\overline{\mathbf{C}}}_n).$$

For the two-sample case, let ${\widehat{\mathit{\mu }}}_i\triangleq {\sum }_{j=1}^{n_i}{\alpha }_{\mathit{\text{ij}}}{\mathbf{x}}_{\mathit{\text{ij}}}={\mathit{\mu }}_i+{\sum }_{j=1}^{n_i}{\alpha }_{\mathit{\text{ij}}}{w}_{\mathit{\text{ij}}}{\mathbf{\Sigma }}_i^{\frac{1}{2}}{\mathbf{z}}_{\mathit{\text{ij}}}$ be the robust estimates of the unknown locations ${\mathit{\mu }}_i$, i = 1, 2, respectively, where ${\alpha }_{\mathit{\text{ij}}}={({\sum }_{j=1}^{n_i}\left|\right|{\mathbf{x}}_{\mathit{\text{ij}}}-{\overline{\mathbf{x}}}_i|{|}^{-2})}^{-1}\left|\right|{\mathbf{x}}_{\mathit{\text{ij}}}-{\overline{\mathbf{x}}}_i|{|}^{-2}$ and ${\overline{\mathbf{x}}}_i={\sum }_{j=1}^{n_i}{\mathbf{x}}_{\mathit{\text{ij}}}/n_i$ being the sample means. Similarly, the two Tyler’s M estimators ${\overline{\mathbf{M}}}_{1n_1}$ and ${\overline{\mathbf{M}}}_{2n_2}$ are the solutions to the equations ${\overline{\mathbf{M}}}_{i{n}_i}=p{n}_i^{-1}{\sum }_{j=1}^{n_i}{(({\mathbf{x}}_{\mathit{\text{ij}}}-{\widehat{\mathit{\mu }}}_{i2})\prime {\overline{\mathbf{M}}}_{i{n}_i}^{-1}({\mathbf{x}}_{\mathit{\text{ij}}}-{\widehat{\mathit{\mu }}}_{i2}))}^{-1}({\mathbf{x}}_{\mathit{\text{ij}}}-{\widehat{\mathit{\mu }}}_{i2})({\mathbf{x}}_{\mathit{\text{ij}}}-{\widehat{\mathit{\mu }}}_{i2})\prime$, $\text{tr}{\overline{\mathbf{M}}}_{i{n}_i}=p$. Our adjusted test statistic for the equality hypothesis in (2.9) is denoted by (2.16) $${\overline{T}}_{2,\mathit{\text{tr}}}=\frac{p\text{tr}{({\overline{\mathbf{M}}}_{1n_1}^{-1}{\overline{\mathbf{M}}}_{2n_2})}^2}{\stackrel{2}{\text{tr}}({\overline{\mathbf{M}}}_{1n_1}^{-1}{\overline{\mathbf{M}}}_{2n_2})}-\frac{n_1-1}{n_1-p-1}-\frac{p}{n_2-1}.$$(2.16)

The asymptotic null distributions of $\{{\overline{T}}_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}},{\overline{T}}_{2,\mathit{\text{tr}}}\}$ are listed below, which are direct conclusions of our newly developed theory on the CLT for LSSs of mean-adjusted Tyler’s M estimator detailed in Theorem 4.3.

Theorem 2.6.

In the large-dimensional asymptotic framework and under the null hypothesis in (2.3), we have $$\begin{array}{c}{\overline{T}}_{1,\mathit{\text{tr}}}-p\left(1+\frac{p}{n-1}\right)\stackrel{D}{\to }N\left(\frac{3c-c^2}{1-c},4c^2\right),\\ {\overline{T}}_{1,\mathrm{\text{log}}}+p-(p-n+1)\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left(1-\frac{p}{n-1}\right)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\stackrel{D}{\to }N\left(\frac{1}{2}\mathrm{\text{log}}\hspace{0.17em}(1-c)-\frac{c}{1-c},-2\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(1-c)-2c\right).\end{array}$$

Theorem 2.7.

In the large-dimensional asymptotic framework and under the null hypothesis in (2.9), we have $p{\overline{T}}_{2,\mathit{\text{tr}}}\stackrel{D}{\to }N(\mu ,{\sigma }^2)$, where the limiting mean and limiting variance are the same as those given in Theorem 2.4.

Comparing Theorems 2.6 and 2.7 with Theorems 2.1 and 2.4, one can see that the three adjusted statistics $\{{\overline{T}}_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}},{\overline{T}}_{2,\mathit{\text{tr}}}\}$ share the same null distributions with $\{T_{1,\mathit{\text{tr}}},T_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}},T_{2,\mathit{\text{tr}}}\}$ except that we have replaced the sample sizes $(n,n_1,n_2)$ with $(n-1,n_1-1,n_2-1)$ in the centering terms. We will discuss the theoretical foundation of such substitution in Section 4.4.

Remark 2.7.

Note that in Kent and Tyler (1991), the authors studied the existence and uniqueness of simultaneous M-estimates of multivariate location and scatter. These estimates are defined to minimize the objective function $L(\mathit{\mu },\Sigma )={\sum }_{i=1}^n\rho \left\{{({\mathbf{x}}_i-\mathit{\mu })}^{\prime }{\Sigma }^{-1}({\mathbf{x}}_i-\mathit{\mu })\right\}+\frac{1}{2}n\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left|\Sigma \right|,$ and thus have the following form (2.17) $$\{\begin{array}{c}{\widehat{\mathit{\mu }}}_M=\sum _{i=1}^n{v}_i{\mathbf{x}}_i/\sum _{i=1}^n{v}_i,\hfill \\ {\widehat{\Sigma }}_M=\frac{1}{n}\sum _{i=1}^n{v}_i\left({\mathbf{x}}_i-{\widehat{\mathit{\mu }}}_M\right){\left({\mathbf{x}}_i-{\widehat{\mathit{\mu }}}_M\right)}^{\prime },\hfill \end{array}$$(2.17) where $v_i=u(s_i),u(s)=2{\rho }^{\prime }(s)$ and $s_i={({\mathbf{x}}_i-{\widehat{\mathit{\mu }}}_M)}^{\prime }{\widehat{\Sigma }}_M^{-1}({\mathbf{x}}_i-{\widehat{\mathit{\mu }}}_M)$. By taking $\rho (s)=\frac{1}{2}\mathrm{\text{log}}\hspace{0.17em}s$, the solution to (2.17) corresponds to Tyler’s M estimator ${\widehat{\Sigma }}_M$ and another estimate of location ${\widehat{\mathit{\mu }}}_M$. However, this choice of $\rho (s)$ does not satisfy the Condition M outlined in Kent and Tyler (1991) for the existence and uniqueness of the solution to (2.17). Although there is potential for further exploration in this direction for simultaneous estimation when the true location is unknown, such theoretical analysis is currently beyond the scope of this article.

Nevertheless, to understand the simultaneous estimates better for future investigations, we have empirically implemented an iterative algorithm to obtain a solution $({\widehat{\mathit{\mu }}}_M,{\widehat{\Sigma }}_M)$ to (2.17), temporarily setting aside the issue of existence and uniqueness of the solution. Detailed simulations are presented at the end of the supplementary material, and the simulation results seem to imply that such ${\widehat{\mathit{\mu }}}_M$ is asymptotically equivalent to our explicitly defined location estimate $\widehat{\mathit{\mu }}$, in terms of the ${\mathcal{l}}_2$-norm.

3 Simulations and Real Data Analysis

In this section, we conduct simulation studies to compare the performance of the proposed tests with existing ones, and illustrate the proposed methodology via a real data example.

3.1 Simulations for Testing (2.3)

We numerically examine the performance of the proposed tests $T_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\mathit{\text{tr}}},T_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$, and ${\overline{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$. For comparison, we also conduct John’s test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ from Theorem 2.3 and the likelihood ratio test ${\tilde{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$ from Bai et al. (2009), both of which are based on the sample covariance matrix. Simulation data are generated from the following two models.

Model (I). The underlying population is multivariate normal $\mathbf{x}\sim \mathcal{N}(\mathbf{0},\mathbf{\Sigma })$ where $\mathbf{\Sigma }={\mathbf{I}}_p+{\mathbf{\Sigma }}_0$ with ${\mathbf{\Sigma }}_0=({\theta }^{|i-j|})$. We set θ = 0 under the null hypothesis and $\theta =0.2$ under the alternative.

Model (II). The underlying population is a scale mixture, that is, $\mathbf{x}=0.5\mathcal{N}(\mathbf{0},\mathbf{\Sigma })+0.5\mathcal{N}(\mathbf{0},4\mathbf{\Sigma })$. The matrix $\mathbf{\Sigma }$ is the same as the one in Model (I).

The population dimension p is set to be p = 100, 200, 400 and the dimension-to-sample size ratio is c_n = 0.25 and 0.5. The nominal level is fixed at 5%. The test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ depends on the mixing parameters of Model (II), which are treated as known in our simulations. All statistics are averaged from 10,000 independent replications.

Empirical sizes of the six tests are summarized in Table 1. The results show that the four Tyler’s M based tests $T_{1,\mathit{\text{tr}}},T_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}},{\overline{T}}_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$, and the test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ all have reasonable sizes under both models. Meanwhile, Tyler’s M based tests have slight size inflation when p and n are small. As the dimensions grow, their sizes reach the nominal level of 0.05. For the test ${\tilde{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$, it has accurate sizes under Model (I) but suffers from serious size distortions under Model (II). This is reasonable since the test ${\tilde{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$ is only justified for the independent components model.

Table 1 Empirical sizes in percentiles of the six tests $T_{1,\mathit{\text{tr}}},T_{1,\mathrm{\text{log}}},{\overline{T}}_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\mathrm{\text{log}}},{\tilde{T}}_{1,\mathit{\text{tr}}}$, and ${\tilde{T}}_{1,\mathrm{\text{log}}}$ at 5% significance.

		Normal distribution						Scale mixture
cn	n	$T_{1,\mathit{\text{tr}}}$	$T_{1,\mathrm{\text{log}}}$	${\overline{T}}_{1,\mathit{\text{tr}}}$	${\overline{T}}_{1,\mathrm{\text{log}}}$	${\tilde{T}}_{1,\mathit{\text{tr}}}$	${\tilde{T}}_{1,\mathrm{\text{log}}}$	$T_{1,\mathit{\text{tr}}}$	$T_{1,\mathrm{\text{log}}}$	${\overline{T}}_{1,\mathit{\text{tr}}}$	${\overline{T}}_{1,\mathrm{\text{log}}}$	${\tilde{T}}_{1,\mathit{\text{tr}}}$	${\tilde{T}}_{1,\mathrm{\text{log}}}$
	200	6.73	6.09	6.70	6.07	5.31	5.06	7.04	6.38	6.78	6.28	8.01	100
0.25	400	5.58	5.47	5.51	5.31	4.89	4.90	6.39	6.15	6.54	5.97	7.23	100
	800	5.23	5.48	5.11	5.33	4.89	5.14	5.55	5.43	5.51	5.29	6.38	100
	200	6.72	5.77	6.27	5.44	4.99	5.19	6.90	5.88	6.60	5.63	6.77	100
0.50	400	5.85	5.51	5.72	5.42	4.84	4.89	5.80	5.55	5.65	5.37	5.93	100
	800	5.45	4.99	5.29	5.01	4.90	4.80	5.46	5.11	5.35	5.11	5.50	100

Table 2 depicts the empirical powers of the six tests. The results of ${\tilde{T}}_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}}$ under the mixture model are excluded due to its size inflation. Under Model (I), the powers of these tests are comparable. That is, the tests based on Tyler’s M estimator keep the same efficiency as those based on the sample covariance matrix, and all achieve full power as the dimensions grow. When the underlying population shifts from the normal distribution to the scale mixture, that is, Model (II), the four Tyler’s M based tests remain well-behaved, but the test ${\tilde{T}}_{1,\mathit{\text{tr}}}$ shows a significant loss of power, as predicted by Theorem 2.3. For instance, when $(c_n,n)=(0.5,400)$, the power of ${\tilde{T}}_{1,\mathit{\text{tr}}}$ declines from 99.08 to 50.40. This demonstrates that Tyler’s M based tests can gain satisfactory robustness against heavy-tailed and tail-dependent data.

Table 2 Empirical powers in percentiles of the six tests $T_{1,\mathit{\text{tr}}},T_{1,\mathrm{\text{log}}},{\overline{T}}_{1,\mathit{\text{tr}}},{\overline{T}}_{1,\mathrm{\text{log}}},{\tilde{T}}_{1,\mathit{\text{tr}}}$, and ${\tilde{T}}_{1,\mathrm{\text{log}}}$ at 5% significance.

		Normal distribution						Scale mixture
cn	n	$T_{1,\mathit{\text{tr}}}$	$T_{1,\mathrm{\text{log}}}$	${\overline{T}}_{1,\mathit{\text{tr}}}$	${\overline{T}}_{1,\mathrm{\text{log}}}$	${\tilde{T}}_{1,\mathit{\text{tr}}}$	${\tilde{T}}_{1,\mathrm{\text{log}}}$	$T_{1,\mathit{\text{tr}}}$	$T_{1,\mathrm{\text{log}}}$	${\overline{T}}_{1,\mathit{\text{tr}}}$	${\overline{T}}_{1,\mathrm{\text{log}}}$	${\tilde{T}}_{1,\mathit{\text{tr}}}$	${\tilde{T}}_{1,\mathrm{\text{log}}}$
	200	62.64	54.85	62.15	53.98	62.34	54.84	62.09	54.91	61.32	54.38	43.09	–
0.25	400	98.55	96.74	98.47	96.62	98.75	96.94	98.63	96.79	98.57	96.68	79.26	–
	800	100	100	100	100	100	100	100	100	100	100	99.00	–
	200	64.79	49.27	63.39	47.67	64.90	48.84	64.34	49.25	63.03	48.32	29.50
0.50	400	98.93	93.90	98.81	93.64	99.08	94.12	99.06	93.98	99.03	93.70	50.40	–
	800	100	100	100	100	100	100	100	100	100	100	78.81	–

NOTE: The number of independent replications is 10,000.

3.2 Simulations for Testing (2.9)

This section evaluates the performance of $T_{2,\mathit{\text{tr}}}$ and ${\overline{T}}_{2,\mathit{\text{tr}}}$. The population dimension is set to be p = 100, 200, 400 and the dimension-to-sample size ratios are $(c_{n1},c_{n2})=(0.25,0.5)$ and $(0.25,0.25)$. Our first experiment is to examine the robustness of the proposed tests by comparing them with the one constructed from sample covariance matrix (Liu et al. 2014), referred to as ${\tilde{T}}_{2,\mathit{\text{tr}}}$. The following two models are considered.

• Model (III). The two populations are multivariate normal such that ${\mathbf{x}}_1\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_p)$ and ${\mathbf{x}}_2\sim \mathcal{N}(\mathbf{0},\Gamma )$, where $\Gamma =\text{diag}({\gamma }_1,\dots ,{\gamma }_p)$ is diagonal with ${\gamma }_j=1-\theta +2\theta (j-1)/(p-1)$. We set θ = 0 ($\Gamma ={\mathbf{I}}_p$) under the null hypothesis and $\theta =0.3$ under the alternative.

• Model (IV). The first population is standard normal ${\mathbf{x}}_1\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_p)$ and the second is a contaminated normal formulated by a scale mixture ${\mathbf{x}}_2\sim (1-\epsilon )\mathcal{N}(\mathbf{0},\Gamma )+\epsilon \mathcal{N}(\mathbf{0},4{\mathbf{I}}_p),$ where $\epsilon \in (0,1)$ represents the proportion of contamination and the matrix $\Gamma$ is defined as in Model (III). We set the proportion to be fixed at $\epsilon =0.05$ and hence, under the null hypothesis ($\Gamma ={\mathbf{I}}_p$), the population ${\mathbf{x}}_2$ generates a dataset containing 5% outliers.

Table 3 depicts empirical sizes and powers of the three tests under Models (III) and (IV). It shows that, despite having small size inflation, the tests $T_{2,\mathit{\text{tr}}}$ and ${\overline{T}}_{2,\mathit{\text{tr}}}$ are both robust against outliers. However, the test ${\tilde{T}}_{2,\mathit{\text{tr}}}$ shows a significant size distortion when outliers are present. This again demonstrates the robustness of Tyler’s M based tests in comparison to sample covariance matrix based tests.

Table 3 Empirical sizes and powers in percentiles of $T_{2,\mathit{\text{tr}}},{\overline{T}}_{2,\mathit{\text{tr}}}$, and ${\tilde{T}}_{2,\mathit{\text{tr}}}$ at 5% significance under Models (III) and (IV).

		Normal distribution						Contaminated normal distribution
$(c_{n1},c_{n2})$	p	Size			Power			Size			Power
		$T_{2,\mathit{\text{tr}}}$	${\overline{T}}_{2,\mathit{\text{tr}}}$	${\tilde{T}}_{2,\mathit{\text{tr}}}$	$T_{2,\mathit{\text{tr}}}$	${\overline{T}}_{2,\mathit{\text{tr}}}$	${\tilde{T}}_{2,\mathit{\text{tr}}}$	$T_{2,\mathit{\text{tr}}}$	${\overline{T}}_{2,\mathit{\text{tr}}}$	${\tilde{T}}_{2,\mathit{\text{tr}}}$	$T_{2,\mathit{\text{tr}}}$	${\overline{T}}_{2,\mathit{\text{tr}}}$	${\tilde{T}}_{2,\mathit{\text{tr}}}$
	100	6.29	6.51	5.41	39.32	39.20	38.18	6.79	6.94	40.21	34.45	34.00	81.27
$(0.25,0.50)$	200	5.98	6.06	5.67	82.32	82.37	82.88	5.80	5.94	87.51	76.67	76.57	99.99
	400	5.54	5.49	5.32	99.98	99.97	99.98	5.57	5.59	100	99.89	99.88	100
	100	6.20	6.33	5.39	57.22	57.09	56.81	6.19	6.33	60.61	50.49	50.33	96.95
$(0.25,0.25)$	200	5.90	5.78	5.60	97.33	97.29	97.47	5.62	5.68	98.68	94.13	94.24	100
	400	5.51	5.47	5.25	100	100	100	5.60	5.59	100	100	100	100

NOTE: The number of independent replications is 10,000.

Our second experiment is designed such that the difference between the two shape matrices only lies in a finite-dimensional subspace. Hence, we adopt the spiked population model (Johnstone 2001) as follows.

• Model (V). ${\mathbf{\Sigma }}_1\propto {\mathbf{I}}_p+\theta \mathbf{\text{uu}}\prime \text{and} {\mathbf{\Sigma }}_2\propto {\mathbf{I}}_p+\theta \mathbf{\text{vv}}\prime ,\theta \in (-1,+\infty ),$ $\mathbf{u}$ and $\mathbf{v}$ are two unit vectors.

For this model, when θ is close to –1, the matrix ${\mathbf{\Sigma }}_j$ (j = 1, 2) has an eigenvalue close to zero, implying the existence of collinearity among the components of the population. When θ is large or even divergent along with the dimension p, a dominant principle component of the population exists. Hence, our two-sample test detects whether the direction of collinearity or the dominant component of the population changes over spatial scales and through time. Moreover, the common base component ${\mathbf{I}}_p$ in ${\mathbf{\Sigma }}_1$ and ${\mathbf{\Sigma }}_2$ can be extended to any positive definite matrix with bounded eigenvalues.

The experimental settings are as follows. The parameter $\theta =-1+1/\sqrt{p}$ for the case of collinearity and $\theta =\sqrt{p}$ for the case of the dominant factor. The vectors are $\mathbf{u}={\mathbf{w}}_1$ and $\mathbf{v}=\hspace{0.17em}\mathrm{\text{cos}}\hspace{0.17em}(\alpha ){\mathbf{w}}_1+\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}(\alpha ){\mathbf{w}}_2$, where ${\mathbf{w}}_1$ and ${\mathbf{w}}_2$ are orthogonal vectors chosen randomly from the unit sphere in ${\mathbb{R}}^p$. The parameter α which represents the angle between $\mathbf{u}$ and $\mathbf{v}$ takes values in $[0,\pi /4]$. Samples are drawn from normal populations with dimensions $(p,n_1,n_2)=(100,200,400)$, (200, 400, 800) and (300, 600, 1200). Empirical powers of $T_{2,\mathit{\text{tr}}}$ and ${\overline{T}}_{2,\mathit{\text{tr}}}$ under both cases at 5% significance from 10,000 independent replications are plotted in Figure 1. It shows that the powers grow to 1 as the angle α and/or the dimensions $(p,n_1,n_2)$ increase, which confirms the capacity of $T_{2,\mathit{\text{tr}}}$ and ${\overline{T}}_{2,\mathit{\text{tr}}}$ to distinguish the two shape matrices.

Fig. 1 Empirical powers of $T_{2,\mathit{\text{tr}}}$ and ${\overline{T}}_{2,\mathit{\text{tr}}}$ for $\theta =-1+1/\sqrt{p}$ (Left panel) and $\theta =\sqrt{p}$ (Right panel) under Model (V) at 5% significance (red horizontal line) from 10,000 independent replications. The dimensions are $(p,n_1,n_2)=(100,200,400)$, (200, 400, 800) and $(300,600,1200)$.

3.3 A Real Data Example

To illustrate the proposed methodology, we apply the proposed tests to an empirical analysis of Fama-French 49 industrial portfolios. The dataset containing average value-weighted daily returns of the portfolios from Years 2002 to 2019 can be downloaded from $\_$ library.html\#Benchmarks" class="url" >http://mba.tuck.dartmouth.\edu/pages/faculty/ken.french/data{}$\_$ library.html\#Benchmarks. To avoid the weekend effect in the stock market, we consider only data subset consisting of every Wednesday’s returns. Since there was a financial crisis in 2008–2009, we focus on two segments data: one from 2002 to 2007 and the other one from 2010 to 2019. Thus, the sample sizes of the data subsets from 2002–2007 and 2010–2019 are 310 and 515, respectively. We apply the proposed tests to examine whether the shape of data changes over time during these two periods.

To demonstrate the proposed model in (2.14) is more general than the ICM in practice, we conduct an exploratory data analysis to illustrate the ICM does not fit the datasets well. Define $T_1=\mathbb{E}\left|\right|\mathbf{x}-\mathbb{E}(\mathbf{x})|{|}^2/p$ and $T_2=\mathbb{E}\left|\right|\mathbf{x}-\mathbb{E}(\mathbf{x})|{|}^4/p^2.$ Then it can be shown that under assumptions imposed in Section 2.4, if $\mathbf{x}$ has stochastic representation (2.14), $T_2/T_1^2\to \mathbb{E}(w^4)/{\{\mathbb{E}(w^2)\}}^2$ almost surely as $p\to \infty$. Note that the ICM corresponds to a constant w, and therefore the ratio $\mathbb{E}(w^4)/{\{\mathbb{E}(w^2)\}}^2=1$. In general, the ratio $\mathbb{E}(w^4)/{\{\mathbb{E}(w^2)\}}^2$ is greater than 1. This motivates us to examine scatterplots of ($T_1^2,T_2$) based on its sample version.

To construct the scatterplots, we generate B = 200 bootstrap samples and calculate the sample version of ($T_1^2,T_2$) for each bootstrap sample. Specifically, for a given bootstrap sample $\{{\mathbf{x}}_j^{*},j=1,\dots ,n\}$. Denote ${\widehat{T}}_2^{*}=n^{-1}{\sum }_{j=1}^n\left|\right|{\mathbf{x}}_j^{*}-{\overline{\mathbf{x}}}^{*}|{|}^4/p^2$, and ${\widehat{T}}_1^{*}=n^{-1}{\sum }_{j=1}^n\left|\right|{\mathbf{x}}_j^{*}-{\overline{\mathbf{x}}}^{*}|{|}^2/p$, where ${\overline{\mathbf{x}}}^{*}$ is the sample mean of the bootstrap sample. We obtain $(T_1^{*2},T_2^{*})$ for each bootstrap sample. The scatterplots of $(T_1^{*2},T_2^{*})$ based on the B = 200 bootstrap samples are depicted in Figure 2. If ICM is valid, one would expect the points $(T_1^{*2},T_2^{*})$ scattering around the straight line y = x since the ratio of $T_2^{*}$ to $T_1^{*2}$ should approximately equal 1 under the ICM model assumption. Figure 2 clearly shows that the ratio of $T_2^{*}$ to $T_1^{*2}$ is dramatically greater than 1.

Fig. 2 Scatterplots of the statistics $\{(T_1^{*2},T_2^{*})\}$ for the data from 2002 to 2007 (left panel) and from 2010 to 2019 (right panel). These statistics are calculated from 200 bootstrap samples.

We now examine whether the shape matrix of the 49 assets would change over Years 2002–2007. To this end, we apply the proposed testing procedure for equality of two shape matrices to examine whether the shape matrix for data from 2002 to 2003 equals the shape matrix for data from 2004 to 2005. The resulting p-value is 0.1745, which implies that the shape does not have dramatically change. We further apply the testing procedure for data from 2004 to 2005 versus data from 2006 to 2007, and the resulting p-value is less than 0.0001. This implies that the shape of data has been changing from 2004–2005 to 2006–2007.

Similarly, we test whether the shape matrix of the 49 assets would change over the period of 2010–2019. We apply the testing procedure for data subsets over every two-year period. The resulting p-values for data from 2010–2011 versus 2012–2013, 2012–2013 versus 2014–2015, 2014–2015 versus 2016–2017, 2016–2017 versus 2018–2019 are 0.0030, 0.0112, 0.0037, and 0.0007, respectively. This implies the shape matrix of the data varies over this period.

4 Theoretical Results on the Limiting Spectrum of Tyler’s M Estimator from a General Population

Recall the definitions of our proposed statistics $T_{1,\mathit{\text{tr}}},T_{1,\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}},T_{2,\mathit{\text{tr}}}$ and their adjusted versions, which are all functionals of the eigenvalues of Tyler’s M estimator (or their F-ratio). To our best knowledge, limiting properties of the eigenvalue distributions of large-dimensional Tyler’s M estimator with general populations have not been explored fully yet and almost all existing works in the literature on M-estimators assume elliptical distributions, see Maronna (1976), Tyler (1987), Couillet, Pascal, and Silverstein (2015), Kammoun et al. (2016), Zhang, Cheng, and Singer (2016), Goes, Lerman, and Nadler (2020), etc. Note that by assuming the parameter of dimension p grows to infinity proportionally to the sample size n (also referred to as Marčenko-Pastur asymptotic regime in random matrix theory), it has been established only recently that for those elliptical distributions, the estimators of Maronna’s M (Couillet, Pascal, and Silverstein 2015), Tyler’s M (Zhang, Cheng, and Singer 2016) and regularized Tyler’s M (Kammoun et al. 2016) all behave similarly and are asymptotically equivalent to the sample covariance matrix from normal distribution. Such equivalence gives the first-order approximation. That is, the ESDs of these robust estimators from elliptical distributions all converge to the renowned Marčenko-Pastur law (Marčenko and Pastur 1967). As for the second-order asymptotics for these robust scatters, existing work is even limited and the only available reference is Couillet, Kammoun, and Pascal (2016), where the authors focused on the fluctuation of bilinear forms of regularized Tyler’s M estimator, also from elliptical populations, and it was again established to be the same as the sample covariance matrix asymptotically. However, in the literature of random matrix theory, the convergence rate of bilinear forms of certain random matrix models is $O_p\left(\sqrt{p}\right)$ (Bai, Miao, and Pan 2007; Couillet, Kammoun, and Pascal 2016), which is in general much slower than the convergence rate of LSSs Bai and Silverstein 2004; Pan and Zhou 2008; Zheng 2012. Therefore, unlike the case for dealing with spectral norms (Couillet, Pascal, and Silverstein 2015; Kammoun et al. 2016; Zhang, Cheng, and Singer 2016) or bilinear forms (Couillet, Kammoun, and Pascal 2016), the LSSs of Tyler’s M estimator may fluctuate differently from those of a sample covariance matrix asymptotically.

Therefore, this section aims to develop new theoretical results on the convergence of ESDs (first-order asymptotic) and the CLT for LSSs (second-order asymptotic) of large-dimensional Tyler’s M estimator from general populations. Note that the technique developed in this article is different from the one in Zhang, Cheng, and Singer (2016). The latter depends on the properties of Gaussian and elliptical distributions, say, some concentration bounds, and thus cannot be directly extended to other distributions. When it comes to the second-order asymptotics, the processing becomes more intricate as it involves obtaining high-precision approximations of certain moments of a resolvent matrix, rather than just comparing two matrices in spectral norms. Our strategy is mainly on the method of moment and a comparison technique. To be concrete, the method of moment helps to get rid of the particular distribution of the underlying population, and the comparison technique allows us to compare Tyler’s M estimator with sample spatial-sign covariance matrix for the first time. We investigate the inner structure of their resolvent matrices and finally show that the two share the same first and second order asymptotics in some sense.

4.1 Notations and Assumptions

First, we introduce some notations that are repeatedly used in this article. Let ${\mathbf{A}}_p$ be a p × p symmetric or Hermitian matrix with eigenvalues ${({\mathcal{l}}_i)}_{1\le i\le p}$. The ESD of ${\mathbf{A}}_p$ is defined as the random measure $F^{{\mathbf{A}}_p}=p^{-1}{\sum }_{i=1}^p{\delta }_{{\mathcal{l}}_i},$ where ${\delta }_{{\mathcal{l}}_i}$ is the Dirac mass at point ${\mathcal{l}}_i$. If the ESD sequence $\left\{F^{{\mathbf{A}}_p}\right\}$ has a limit when $p\to \infty$, this limit is referred to as the limiting spectral distribution (LSD) of the sequence. For a given analytical function f(x) defined on $\mathbb{R}$, we call $\int f(x)d{F}^{{\mathbf{A}}_p}(x)=p^{-1}{\sum }_{i=1}^{p}f({\mathcal{l}}_i)$ an LSS of the random matrix ${\mathbf{A}}_p$ associated with the test function f. The Stieltjes transform of a probability measure G is defined as $m_G(z)=\int {(x-z)}^{-1}\mathit{\text{dG}}(x),z\in {\mathbb{C}}^{+}\triangleq \{z\in \mathbb{C}:\mathfrak{I}(z)>0\}$. This definition can be extended to the whole complex plane except the support of G.

Next, we recall our assumptions and definitions for easy reading. Let ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\in {\mathbb{R}}^p$ be a random sample having the following stochastic representation, (4.1) $${\mathbf{x}}_j\stackrel{d}{=}{w}_j{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j,\hspace{1em}\text{tr}(\mathbf{\Sigma })=p.$$(4.1)

Our main assumptions on this model are listed below.

Assumption

(1). Both the sample size n and population dimension p tend to infinity such that $p\to \infty ,n\to \infty ,c_n\triangleq p/n\to c\in (0,1)$.

Assumption

(2). The ESD of $\mathbf{\Sigma }$, denoted as H_p, has a bounded support, that is $\text{Supp}(H_p)\subset [a,b]\subset (0,\infty )$, and converges weakly to a probability distribution H.

Assumption

(3). The coordinates $(z_{\mathit{\text{ij}}})$ of ${\mathbf{z}}_j$ are iid continuous random variables having bounded density and satisfying the following moment conditions for some $\Delta >0,\mathbb{E}(z_{\mathit{\text{ij}}})=0,\mathbb{E}(z_{\mathit{\text{ij}}}^2)=1,\mathbb{E}(z_{\mathit{\text{ij}}}^4)=\tau ,\mathbb{E}|z_{\mathit{\text{ij}}}{|}^{4+\Delta }<\infty$.

In the following Sections 4.2 and 4.3, we establish the LSD and the CLT for LSSs of Tyler’s M estimator ${\mathbf{C}}_n$ defined by ${\mathbf{C}}_n=p{n}^{-1}{\sum }_{j=1}^n{({\mathbf{x}}_j^{\prime }{\mathbf{C}}_n^{-1}{\mathbf{x}}_j)}^{-1}{\mathbf{x}}_j{\mathbf{x}}_j^{\prime }$, $\text{tr}{\mathbf{C}}_n=p$, where ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$ is a random sample from (4.1) satisfying Assumptions (1)–(3). Finally, our results are generalized to cover the case of unknown location in Section 4.4.

Remark 4.1.

We discuss here the reason for imposing Assumption (1). Initially, the definition of Tyler’s M estimator implies that the matrix ${\mathbf{C}}_n$ is well-defined only when p < n. As a result, we set the constraint c < 1. It’s worth noting that recent works, such as Zhang, Cheng, and Singer (2016) and Goes, Lerman, and Nadler (2020), have explored the properties of Tyler’s M estimator within this very regime. To bypass the constraint c < 1, one might consider alternative matrix models like the regularized Tyler’s M estimator, as proposed by Ollila and Tyler (2014). While the techniques developed in this article could potentially address non-explicitly defined random matrix models, such discussions are beyond the scope of this article and could be a good topic for future research.

4.2 Limiting Spectral Distribution of ${\mathbf{C}}_n$

The following theorem is on the LSD of Tyler’s M estimator ${\mathbf{C}}_n$.

Theorem 4.1.

Under Assumptions (1)–(3), almost surely, the ESD $F^{{\mathbf{C}}_n}$ of Tyler’s M estimator ${\mathbf{C}}_n$ converges weakly to a probability measure $F^{c,H}$, whose Stieltjes transform m(z) is a solution to (4.2) $$m(z)=\int \frac{1}{t\{1-c-\mathit{\text{czm}}(z)\}-z}\mathrm{d}H(t),\hspace{1em}z\in {\mathbb{C}}^{+}.$$(4.2)

The solution is unique in the set $\{m(z)\in \mathbb{C}:-(1-c)/z$ $+\mathit{\text{cm}}(z)\in {\mathbb{C}}^{+}\}$.

Remark 4.2.

Theorem 4.1 extends the main conclusion in Zhang, Cheng, and Singer (2016), where the same LSD was derived under elliptical distributions. Indeed, according to Remark 2.1(b), elliptical distributions require those atom random vectors $\left\{{\mathbf{z}}_j\right\}$ in (4.1) to be standard multivariate normal $\mathcal{N}(\mathbf{0},{\mathbf{I}}_p)$. In that case, any finite moment of $\left\{z_{\mathit{\text{ij}}}\right\}$ exists. Our Theorem 4.1 is indeed universal. It reaches the same conclusion not only for those standard normal $\left\{{\mathbf{z}}_j\right\}$ but also for any continuous distributions with finite $(4+\Delta )$th moment.

Remark 4.3.

The LSD $F^{c,H}$, characterized by its Stieltjes transform m(z) as defined in (4.2), is a generalized Marčenko-Pastur law, initially introduced in the seminal paper Marčenko and Pastur (1967) to describe the LSD of sample covariance matrices. Hence, Theorem 4.1 illustrates that, in terms of the first-order asymptotics for the ESD, Tyler’s M estimator ${\mathbf{C}}_n$ constructed by samples ${\left\{{\mathbf{x}}_j\right\}}_{j=1}^n$ originating from the model (4.1) with shape matrix $\Sigma$, is asymptotically equivalent to the sample covariance matrix ${\mathbf{S}}_n$ formed by samples drawn from the independent components model with the same population covariance matrix $\Sigma$, that is, ${\mathbf{S}}_n=n^{-1}{\sum }_{i=1}^n{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j\mathbf{z}{\prime }_j{\mathbf{\Sigma }}^{\frac{1}{2}}.$ Indeed, the auxiliary Lemma S1.1, which is crucial for proving our main results, has yielded a stronger conclusion than the current one in Theorem 4.1. Lemma S1.1 in the supplementary material demonstrates that $\left|\right|{\mathbf{C}}_n-{\mathbf{B}}_n\left|\right|\to 0$, almost surely, where ${\mathbf{B}}_n=p{n}^{-1}{\sum }_{j=1}^n\left|\right|{\mathbf{x}}_j|{|}^{-2}{\mathbf{x}}_j{\mathbf{x}}_j^{\prime }$ is the spatial-sign covariance matrix constructed by samples ${\left\{{\mathbf{x}}_j\right\}}_{j=1}^n$. Additionally, it is straightforward to verify that $\left|\right|{\mathbf{B}}_n-{\mathbf{S}}_n\left|\right|\to 0$, almost surely. These facts imply that the three matrices, Tyler’s M estimator ${\mathbf{C}}_n$, the spatial-sign covariance matrix ${\mathbf{B}}_n$, and the sample covariance matrix ${\mathbf{S}}_n$, are asymptotically equivalent in terms of spectral norm.

Remark 4.4.

Note that while Theorem 4.1 is established under Assumption (1), if we allow c approach 0, the result indicates that the eigenvalue’s distribution of Tyler’s M estimator ${\mathbf{C}}_n$ converges to that of the shape matrix $\mathbf{\Sigma }$. This finding aligns with the established result derived in Tyler (1987) for the fixed-p scenario.

Remark 4.5.

Let $\underset{¯}{m}(z)\triangleq \mathit{\text{cm}}(z)-(1-c)/z$, which is the Stieltjes transform of ${\underset{¯}{F}}^{c,H}\triangleq c{F}^{c,H}+(1-c){\delta }_0$. Then, (4.2) can be rewritten using $\underset{¯}{m}(z)$ as (4.3) $$z=-\frac{1}{\underset{¯}{m}(z)}+c\int \mathrm{}\frac{t}{1+t\underset{¯}{m}(z)}\mathit{\text{dH}}(t),\hspace{1em}z\in {\mathbb{C}}^{+}.$$(4.3)

This quantity $\underset{¯}{m}(z)$ will be frequently recalled in the next subsection. For simplicity of notation, we shall sometimes suppress the independent variable z in the functions $\underset{¯}{m}(z)$ and m(z).

4.3 Central Limit Theorem for Linear Spectral Statistics of ${\mathbf{C}}_n$

This section aims to establish the CLT for LSSs of Tyler’s M estimator ${\mathbf{C}}_n$, that is, the statistics of the type $\int f(x)d{F}^{{\mathbf{C}}_n}(x),$ where f(x) is a given test function. To establish our main result, we introduce a deterministic matrix T, which relates to the shape matrix $\mathbf{\Sigma }$ in the following way, (4.4) $$\begin{array}{c}\mathbf{T}\triangleq \mathbf{\Sigma }-\frac{2}{p}{\mathbf{\Sigma }}^2-\frac{\tau -3}{p}{\mathbf{\Sigma }}^{\frac{1}{2}}\text{diag}(\mathbf{\Sigma }){\mathbf{\Sigma }}^{\frac{1}{2}}\\ \hspace{1em}+(\frac{2}{p^2}\text{tr}{\mathbf{\Sigma }}^2+\frac{\tau -3}{p^2}\text{tr}[\mathbf{\Sigma }°\mathbf{\Sigma }\left]\right)\mathbf{\Sigma },\end{array}$$(4.4) where “$°$” denotes the Hadamard product of two matrices. Note that the matrix T is asymptotically equivalent to $\mathbf{\Sigma }$. However, the remaining terms which are of the order $O(1/p)$ in spectral norm will still contribute when we study the CLT for LSSs of ${\mathbf{C}}_n$. Let ${\tilde{H}}_p$ denote the ESD of T. Then, by replacing (c, H) with $(c_n,{\tilde{H}}_p)$ in (4.2), this new equation defines a unique solution, referred as ${\tilde{m}}_0(z)$, and its induced probability measure is referred as $F^{c_n,{\tilde{H}}_p}$. In virtue of this measure, we may normalize the LSS of ${\mathbf{C}}_n$ as (4.5) $$G_n(f)\triangleq p\left(\int f(x)d{F}^{{\mathbf{C}}_n}(x)-\int f(x)d{F}^{c_n,{\tilde{H}}_p}(x)\right).$$(4.5)

In addition, we need the following Assumption (4) for establishing our general CLT, which indeed is required for handling the case when $\mathbb{E}(z_{\mathit{\text{ij}}}^4)=\tau \ne 3$. Note that such an assumption depends not only on the eigenvalues of $\mathbf{\Sigma }$, but also on its eigenvectors.

Assumption

(4). As $p\to \infty$, the limits of the following three auxiliary quantities exist, (4.6) $$\begin{array}{c}{\zeta }_p=\frac{1}{p}\text{tr}[\mathbf{\Sigma }°\mathbf{\Sigma }]\to \zeta ,\\ h_p(z)=\frac{1}{p}\text{tr}\left[\right({\mathbf{\Sigma }}^{\frac{1}{2}}{(\mathbf{\Sigma }-z\mathbf{I})}^{-1}{\mathbf{\Sigma }}^{\frac{1}{2}})°\mathbf{\Sigma }]\to h(z),\\ g_p(z_1,z_2)=\frac{1}{p}\text{tr}\left[\left({\mathbf{\Sigma }}^{\frac{1}{2}}{(\mathbf{\Sigma }-z_1\mathbf{I})}^{-1}{\mathbf{\Sigma }}^{\frac{1}{2}}\right)°\left({\mathbf{\Sigma }}^{\frac{1}{2}}{(\mathbf{\Sigma }-z_2\mathbf{I})}^{-1}{\mathbf{\Sigma }}^{\frac{1}{2}}\right)\right]\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\to g(z_1,z_2),\end{array}$$(4.6) where $z,z_1,z_2\in {\mathbb{C}}^{+}$.

We next show that the linear functional $\{G_n(f)\}$ of the eigenvalues of ${\mathbf{C}}_n$ defined in (4.5) converges weakly to a Gaussian process.

Theorem 4.2.

Suppose that Assumptions (1)–(4) hold with some $\Delta \ge 2$. Let $f_1,\dots ,f_k$ be k functions analytic on an open set that includes the interval $I_c=[\underset{p\to \infty }{\mathrm{\text{lim inf}}}{\lambda }_{{\mathrm{\text{min}}}^{\mathbf{\Sigma }}}{(1-\sqrt{c})}^2$, $\underset{p\to \infty }{\mathrm{\text{lim sup}}}{\lambda }_{{\mathrm{\text{max}}}^{\mathbf{\Sigma }}}{(1+\sqrt{c})}^2].$ Then the random vector $\left\{G_n(f_1),\dots ,G_n(f_k)\right\}$ converges weakly to a Gaussian random vector $\{G_{f_1},\dots ,G_{f_k}\}$. Its mean function is given by (4.7) $$\mathbb{E}(G_f)=-\frac{1}{2\pi \mathrm{i}}{\oint }_{\mathcal{C}}f(z)[{\mu }_1(z)+{\mu }_2(z)]\mathit{\text{dz}},$$(4.7) where (4.8) $$\begin{array}{c}{\mu }_1(z)=\int \frac{c{(\underset{¯}{m}\prime t)}^2\mathit{\text{dH}}(t)}{\underset{¯}{m}{(1+\underset{¯}{m}t)}^3}-\int \frac{2\underset{¯}{m}\prime (1+z\underset{¯}{m})t^2\mathit{\text{dH}}(t)}{{(1+\underset{¯}{m}t)}^2}\\ \hspace{1em}+\int \frac{({\beta }_{2}t-t^2)\mathit{\text{dH}}(t)}{1+\underset{¯}{m}t}\int \frac{2c\underset{¯}{m}\underset{¯}{m}\prime \mathit{\text{tdH}}(t)}{{(1+\underset{¯}{m}t)}^2},\\ \hspace{1em}+(\tau -3)\{\frac{c\underset{¯}{m}\prime }{{\underset{¯}{m}}^2}g\prime \left(\frac{-1}{\underset{¯}{m}},\frac{-1}{\underset{¯}{m}}\right)+\zeta \int \frac{(1+z\underset{¯}{m})t\underset{¯}{m}\prime \mathit{\text{dH}}(t)}{{(1+\underset{¯}{m}t)}^2}\\ \hspace{1em}-\frac{(1+z\underset{¯}{m})\underset{¯}{m}\prime }{{\underset{¯}{m}}^2}h\prime \left(\frac{-1}{\underset{¯}{m}}\right)-\int \frac{c\underset{¯}{m}\prime \mathit{\text{tdH}}(t)}{{(1+\underset{¯}{m}t)}^2}h\left(\frac{-1}{\underset{¯}{m}}\right)\},\end{array}$$(4.8) (4.9) $$\begin{array}{c}{\mu }_2(z)=(z\underset{¯}{m})\prime \{\frac{2}{c(1-c)}z\underset{¯}{m}-\frac{2}{c}{\beta }_{2}z\underset{¯}{m}-\frac{2}{c}z+\frac{2}{c(1-c)}+2\}\\ -\frac{2}{c}(1+z\underset{¯}{m})+(\tau -3)\{\frac{d}{\mathit{\text{dz}}}\left[z\underset{¯}{m}h\right(\frac{-1}{\underset{¯}{m}}\left)\right]\\ -\frac{1}{c}(z\underset{¯}{m})\prime z\underset{¯}{m}\zeta -\frac{1}{c}(1+z\underset{¯}{m})(z\underset{¯}{m})\prime \},\end{array}$$(4.9) in which ${\beta }_2=\int t^2\mathit{\text{dH}}(t)$. Its covariance function is given by (4.10) $$\begin{array}{c}\text{cov}\left(G_f,G_g\right)\\ \hspace{1em}=-\frac{1}{4{\pi }^2}{\oint }_{{\mathcal{C}}_1}{\oint }_{{\mathcal{C}}_2}f(z_1)g(z_2)[{\sigma }_1(z_1,z_2)+(\tau -3){\sigma }_2(z_1,z_2)]d{z}_{1}d{z}_2,\end{array}$$(4.10) where $$\begin{array}{c}{\sigma }_1(z_1,z_2)=\frac{2\underset{¯}{m}\prime (z_1)\underset{¯}{m}\prime (z_2)}{{(\underset{¯}{m}(z_1)-\underset{¯}{m}(z_2))}^2}\\ \hspace{1em}+\frac{2}{c}({\beta }_2+\frac{1}{\underset{¯}{m}(z_1)}+\frac{1}{\underset{¯}{m}(z_2)})(z_1\underset{¯}{m}(z_1))\prime (z_2\underset{¯}{m}(z_2))\prime \\ \hspace{1em}-\frac{2\underset{¯}{m}\prime (z_1)}{c{\underset{¯}{m}}^2(z_1)}(z_2\underset{¯}{m}(z_2))\prime (1+z_1\underset{¯}{m}(z_1))\\ \hspace{1em}-\frac{2\underset{¯}{m}\prime (z_2)}{c{\underset{¯}{m}}^2(z_2)}(z_1\underset{¯}{m}(z_1))\prime (1+z_2\underset{¯}{m}(z_2)),\\ {\sigma }_2(z_1,z_2)=\frac{{\partial }^2}{\partial z_1\partial z_2}\{\mathit{\text{cg}}\left(\frac{-1}{\underset{¯}{m}(z_1)},\frac{-1}{\underset{¯}{m}(z_2)}\right)\\ \hspace{1em}+\frac{\zeta }{c}(1+z_1\underset{¯}{m}(z_1))(1+z_2\underset{¯}{m}(z_2))\\ \hspace{1em}-(1+z_1\underset{¯}{m}(z_1))h\left(\frac{-1}{\underset{¯}{m}(z_2)}\right)\\ \hspace{1em}-(1+z_2\underset{¯}{m}(z_2))h\left(\frac{-1}{\underset{¯}{m}(z_1)}\right)\}.\end{array}$$

Here, the contours in (4.7) and (4.10) (two in (4.10), which may be assumed to be nonoverlapping) are closed, counter-clockwise orientated in the complex plane and enclosing the support of $F^{c,H}$.

Remark 4.6.

Theorem 4.2 indicates that the convergence rate for the LSSs of Tyler’s M estimator is of the order O(n) when the dimension p grows to infinity proportionally with the sample size n, while keeping their ratio away from zero asymptotically. Theoretically, our result does not cover scenarios where p is fixed or p diverges but c = 0. In the fixed-p case, Tyler (1987) demonstrates that under certain conditions, $\sqrt{n}({\mathbf{C}}_n-\Sigma )$ converges in distribution to a matrix normal distribution as $n\to \infty$. This large sample behavior of ${\mathbf{C}}_n$ resembles the classical theory on sample covariance matrix introduced in the book by Anderson (1984). When p diverges but c = 0, random matrix theory suggests that the convergence rate for LSSs becomes the order of $O\left(\sqrt{\mathit{\text{pn}}}\right)$, as observed in the results for sample covariance matrices in Chen and Pan (2015). Particularly, in order to determine the centering term, as well as the limiting mean and limiting variance of the CLT for LSSs, we must identify the dominant terms that will contribute by controlling the remaining negligible terms up to certain orders, which depends heavily on the convergence rates. As far as we know, there is currently no unified theory that can easily accommodate all these asymptotic frameworks including fixed p, $p/n\to 0$ and $p/n\to c>0$. Nevertheless, for practitioners applying the CLT in real data analysis or applications, it’s worth noting that once a dataset is available, one can directly set the value of c equal to p/n. Our CLT result remains valid even if c is small. We have applied this approach in all our simulation experiments and real data analysis provided in Section 3.

Finally, we discuss the results of two special cases covered by Theorem 4.2.

Special case (1): $\mathbb{E}{z}_{\mathit{\text{ij}}}^4=3$. This special case includes the elliptical family of distributions. From Remark 2.1(b), when the sample $\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\}$ comes from the family of elliptical distributions, the atom random vectors ${\mathbf{z}}_1,\dots ,{\mathbf{z}}_n$ in (4.1) reduce to the standard multivariate normal vectors $\mathcal{N}(\mathbf{0},{\mathbf{I}}_p)$, which meets $\mathbb{E}(z_{\mathit{\text{ij}}}^4)=3$. According to Theorem 4.2 and by setting τ = 3, all the impacts from the eigenvectors of $\mathbf{\Sigma }$ disappear. We establish in the following Proposition 4.1 that the centering term $\int f(x)d{F}^{c_n,{\tilde{H}}_p}(x)$ given in (4.16) can also be simplified to a functional of H_p only. In this way, all the quantities including the centering term, limiting mean and limiting variance in Theorem 4.2 are only relevant to the ESD of $\mathbf{\Sigma }$, that is, the information contained in $\left\{H_p\right\}$.

Proposition 4.1.

Suppose that Assumptions (1)–(3) hold with τ = 3. We have $$\begin{array}{c}p\left(\int f(x)d{F}^{c_n,H_p}(x)-\int f(x)d{F}^{c_n,{\tilde{H}}_p}(x)\right)\\ \hspace{1em}\to -\frac{1}{2\pi \mathrm{i}}{\oint }_{\mathcal{C}}f(z){\mu }_3(z)\mathit{\text{dz}},\end{array}$$where $${\mu }_3(z)=2\underset{¯}{m}\prime (z)\int \frac{{\beta }_{2}t-t^2}{{(1+t\underset{¯}{m}(z))}^2}\mathit{\text{dH}}(t)\text{ and }{\beta }_2=\int t^2\mathit{\text{dH}}(t).$$

Special case (2): identity shape matrix. When the shape matrix $\mathbf{\Sigma }={\mathbf{I}}_p$ ($H_p\equiv {\delta }_1$), Assumption (4) holds automatically. Indeed, the three auxiliary quantities defined in (4.6) are exactly equal to their limits given by ζ = 1, $h(z)={(1-z)}^{-1}$ and $g(z_1,z_2)={(1-z_1)}^{-1}{(1-z_2)}^{-1}$. After some algebra by setting $H={\delta }_1$ in Remark 4.5, we have all the limiting quantities in Theorem 4.2 involving the factor $\tau -3$ equal zero, that is, the terms in (4.8), (4.9) and ${\sigma }_2(z_1,z_2)$ all disappear. Moreover, by (4.4), we notice that when $\mathbf{\Sigma }={\mathbf{I}}_p$, the matrix T reduces to the identity matrix ${\mathbf{I}}_p$ which does not depend on the factor $\tau -3$. Thus, for those samples with identity shape matrix, the CLT for LSSs of Tyler’s M estimator ${\mathbf{C}}_n$ is irrelevant to the fourth moment τ. Therefore, it is universal in the sense that it does not depend on the specific underlying distribution of the population.

4.4 Adjustments to the CLT When the Location Vector is Unknown

This section considers the CLT for LSSs of Tyler’s M estimator when the true location vector is unknown. To this end, suppose we have a sequence of iid observations $\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_n\}\in {\mathbb{R}}^p$ having the form (4.11) $${\mathbf{x}}_j\stackrel{d}{=}\mathit{\mu }+w_j{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j,\hspace{1em}\text{tr}(\mathbf{\Sigma })=p,$$(4.11) where $\mathit{\mu }$ is now unknown. We update Assumption (3) of the model (2.1) as follows.

Assumption

(3’). The coordinates $(z_{\mathit{\text{ij}}})$ of ${\mathbf{z}}_j$ are iid continuous random variables having bounded density and satisfying $\mathbb{E}(z_{\mathit{\text{ij}}})=0,\mathbb{E}(z_{\mathit{\text{ij}}}^2)=1,\mathbb{E}(z_{\mathit{\text{ij}}}^4)=\tau ,\mathbb{E}(z_{\mathit{\text{ij}}}^{\gamma })<\infty$ for some $\gamma >8$. The random variables $\left\{w_j\right\}$ are independent of $\left\{{\mathbf{z}}_j\right\}$ and satisfying $\mathbb{E}(w_j^{-r})<\infty ,\mathbb{E}(w_j^r)<\infty$ for some r > 4.

Tyler’s M estimator corresponding to the data sample ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$ in (4.11), denoted as ${\overline{\mathbf{C}}}_n={\overline{\mathbf{C}}}_n(\widehat{\mathit{\mu }})$, is then defined to be the solution of equation (4.12) $${\overline{\mathbf{C}}}_n=\frac{p}{n}\sum _{j=1}^n\frac{({\mathbf{x}}_j-\widehat{\mathit{\mu }})({\mathbf{x}}_j-\widehat{\mathit{\mu }})\prime }{({\mathbf{x}}_j-\widehat{\mathit{\mu }})\prime {\overline{\mathbf{C}}}_n^{-1}({\mathbf{x}}_j-\widehat{\mathit{\mu }})},\hspace{1em}\text{tr}({\overline{\mathbf{C}}}_n)=p,$$(4.12) where $\widehat{\mathit{\mu }}$ is an estimate of the true location $\mathit{\mu }$. We consider a class of location estimates ${\widehat{\mathit{\mu }}}_{\lambda }$ of the form (4.13) $${\widehat{\mathit{\mu }}}_{\lambda }\triangleq \sum _{j=1}^n{\alpha }_{j\lambda }{\mathbf{x}}_j=\mathit{\mu }+\sum _{j=1}^n{\alpha }_{j\lambda }{w}_j{\mathbf{\Sigma }}^{\frac{1}{2}}{\mathbf{z}}_j,$$(4.13) (4.14) $${\alpha }_{j\lambda }\triangleq \frac{\left|\right|{\mathbf{x}}_j-\overline{\mathbf{x}}|{|}^{-\lambda }}{\sum _{j=1}^n\left|\right|{\mathbf{x}}_j-\overline{\mathbf{x}}|{|}^{-\lambda }},$$(4.14) where λ is a pre-determined nonnegative number in $[0,2]$ and $\overline{\mathbf{x}}={\sum }_{j=1}^n{\mathbf{x}}_j/n$ is the sample mean. Notice that the quantitie $\left\{\right||{\mathbf{x}}_j-\overline{\mathbf{x}}|{|}^{-\lambda }\}$ is close to $\left\{{(\sqrt{p}{w}_j)}^{-\lambda }\right\}$ and thus, by the construction, the weight $\left\{{\alpha }_{j\lambda }\right\}$ attempts to approximate the oracle one $w_j^{-\lambda }/{\sum }_{j=1}^n{w}_j^{-\lambda }.$ In particular, ${\widehat{\mathit{\mu }}}_{\lambda }$ reduces to the usual sample mean if we set λ = 0. With such a class of location estimates ${\widehat{\mathit{\mu }}}_{\lambda }$ involved, we establish the relation between the LSSs of ${\overline{\mathbf{C}}}_n({\widehat{\mathit{\mu }}}_{\lambda })$ and those of ${\mathbf{C}}_n$ in the following theorem.

Theorem 4.3.

Suppose that Assumptions (1)-(2)-(3’) hold. Then, we have $$p\int f(x)d{F}^{{\overline{\mathbf{C}}}_n({\widehat{\mathit{\mu }}}_{\lambda })}(x)-p\int f(x)d{F}^{{\mathbf{C}}_n}(x)\to -\frac{1}{2\pi \mathrm{i}}{\oint }_{\mathcal{C}}f(z){\mu }_4(z)\mathit{\text{dz}}$$in probability, where $$\begin{array}{c}{\mu }_4(z)=[z\underset{¯}{m}(z)]\prime \left\{\frac{-1-{\tilde{\alpha }}_{\lambda }-2{\tilde{\alpha }}_{\lambda }z\underset{¯}{m}(z)}{z\underset{¯}{m}(z)+{\tilde{\alpha }}_{\lambda }z\underset{¯}{m}(z)[1+z\underset{¯}{m}(z)]}+{\tilde{\alpha }}_{\lambda }-1\right\},\\ {\tilde{\alpha }}_{\lambda }=\frac{\mathbb{E}{w}^{-2}\mathbb{E}{w}^{2-2\lambda }}{{(\mathbb{E}{w}^{-\lambda })}^2}-1.\end{array}$$

Theorem 4.3

demonstrates that substituting the location estimate ${\widehat{\mathit{\mu }}}_{\lambda }$ for the true location vector $\mathit{\mu }$ in Tyler’s M estimator only introduces a mean shift to the CLT for its LSSs. And the shift is characterized by the function ${\mu }_4(z)$, which depends on certain moments of the random variable w as specified in the parameter ${\tilde{\alpha }}_{\lambda }$. In particular, by setting λ = 2 in Theorem 4.3, we get ${\tilde{\alpha }}_{\lambda }=0$ and the function ${\mu }_4(z)$ then becomes irrelevant to the moments of w and possesses a concise form (4.15) $$\begin{array}{c}{\mu }_4(z)=\frac{(1+z\underset{¯}{m}(z))(\underset{¯}{m}(z)+z\underset{¯}{m}\prime (z))}{-z\underset{¯}{m}(z)}\\ =p\int \frac{1}{x-z}\left[d{F}^{c_{n-1},{\tilde{H}}_p}(x)-d{F}^{c_n,{\tilde{H}}_p}(x)\right]+o(1),\end{array}$$(4.15) where $F^{c_n,{\tilde{H}}_p}$ and $F^{c_{n-1},{\tilde{H}}_p}$ with $c_{n-1}=p/(n-1)$ are two finite-dimensional proxies for the LSD $F^{c,H}$ with (c, H) replaced by $(c_n,{\tilde{H}}_p)$ and $(c_{n-1},{\tilde{H}}_p)$, respectively, see Zheng, Bai, and Yao (2015). Such ${\mu }_4(z)$ in (4.15) can be canceled out by a simple replacement of the sample size n with n – 1 in the centering term of the LSS defined in (4.5). We summarize such substitution result in the following proposition.

Proposition 4.2.

Suppose that Assumptions (1), (2), (3’), and (4) hold. Denote the normalized LSS of ${\overline{\mathbf{C}}}_n({\widehat{\mathit{\mu }}}_2)$ as (4.16) $${\overline{G}}_n(f)\triangleq p\left(\int f(x)d{F}^{{\overline{\mathbf{C}}}_n({\widehat{\mathit{\mu }}}_2)}(x)-\int f(x)d{F}^{c_{n-1},{\tilde{H}}_p}(x)\right).$$(4.16)

Then the random vector $\left\{{\overline{G}}_n(f_1),\dots ,{\overline{G}}_n(f_k)\right\}$ converges weakly to the same Gaussian random vector $\{G_{f_1},\dots ,G_{f_k}\}$ as that in Theorem 4.2.

Very interestingly, we find for the first time that by substituting the true location vector $\mathit{\mu }$ with its estimate ${\widehat{\mathit{\mu }}}_{\lambda }$ by setting λ = 2 in (4.13), the CLT for LSSs of their corresponding Tyler’s M estimators do not change. The only difference lies in the centering terms for normalizing the LSSs. This is a parallel result to the substitution principle for LSSs of the sample covariance matrix when the unknown mean is estimated by the sample mean, see Zheng, Bai, and Yao (2015).

Supplementary Materials

This supplementary material consists of proofs of the main results developed in Sections 2 and 4, namely Theorems 2.1, 2.2, 2.3, 2.4, 2.5, 4.1, 4.2, 4.3, and Proposition 4.1. Additionally, it includes some additional lemmas for proving these results and simulations for the simultaneous estimates mentioned in Remark 2.7.

Acknowledgments

The authors would like to thank the Associate Editor and referees for their constructive comments, which lead to significant improvement of this work.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

R. Li’s research is supported by NIH grants NIH R01AI170249 and R01AI136664, W. Li’s research is supported by NSFC (No. 11971293, 12141107), and Q. Wang’s research is supported by NSFC (No. 12171099) and Natural Science Foundation of Shanghai (23ZR1406200).