#  Sub-Gaussian High-Dimensional Covariance Matrix Estimation under Elliptical
Factor Model with $2+\varepsilon$th Moment

Yi Ding Faculty of Business Administration, University of Macau, Macau.
Research is supported in part by NSFC 72101226 of China, and
SRG2023-00001-FBA, MYRG-GRG2023-00119-FBA of the University of Macau. Email:
<EMAIL_ADDRESS>Xinghua Zheng Department of ISOM, Hong Kong University of
Science and Technology, Clear Water Bay, Kowloon, Hong Kong. Research is
supported in part by RGC grants GRF-16304521 and GRF 16304019 of the HKSAR.
Email<EMAIL_ADDRESS>
###### Abstract

We study the estimation of high-dimensional covariance matrices under
elliptical factor models with $2+\varepsilon$th moment. For such heavy-tailed
data, robust estimators like the Huber-type estimator in fan2018large can not
achieve sub-Gaussian convergence rate. In this paper, we develop an
idiosyncratic-projected self-normalization (IPSN) method to remove the effect
of heavy-tailed scalar parameter, and propose a robust pilot estimator for the
scatter matrix that achieves the sub-Gaussian rate. We further develop an
estimator of the covariance matrix and show that it achieves a faster
convergence rate than the generic POET estimator in fan2018large.

Keywords: High-dimension, elliptical model, factor model, covariance matrix,
robust estimation

##  1 Introduction

###  1.1 Covariance Matrix Estimation for Heavy-tailed Distributions

There is an extensive interest in robust estimation of the covariance matrices
for heavy-tailed data, which are commonly seen in various fields such as
finance, economics and biology. Huber-type robust estimators (minsker2018sub;
minsker2020robust), and truncation or shrinkage methods (fan2021shrinkage) are
popular approaches in this area. Notably, these studies typically assume that
data have finite 4th moment.

avella2018robust study the high-dimensional covariance matrix estimation for
cases where data have finite 4th or $2+\varepsilon$th moment with
$\varepsilon<2$. Under sparsity assumptions on the covariance matrix or its
inverse matrix, their approach involves applying thresholding techniques
(cai2011adaptive; cai2016estimating) on a robust pilot matrix. The pilot
matrix combines a Huber-type M-estimator of variances and a rank-based robust
correlation matrix estimator such as the marginal Kendall’s tau (liu2012high;
xue2012regularized). Under the high-dimensional setting where the dimension
$p$ and the sample size $n$ both go to infinity, if the data have 4th moment,
avella2018robust show that their estimator achieves the convergence rate of
$\sqrt{(\log p)/n}$, which is the same convergence rate as in the sub-Gaussian
case. However, for the case when only $2+\varepsilon$th moment exists, the
convergence rate reduces to $((\log p)/n)^{\varepsilon/(2+\varepsilon)}$; see
Proposition 6 therein.

###  1.2 Elliptical Factor Models

Due to their flexibility in modeling heavy-tailedness and capturing tail-
dependences, elliptical models are widely used to model financial returns
(ZL11; han2014scale), gene microarrays (avella2018robust) and brain images
(han2018eca). An elliptical model is typically defined as follows:

$\mathbf{y}_{t}=\bm{\mu}+\sqrt{p}\xi_{t}\bm{\Sigma}_{0}^{1/2}\mathbf{u}_{t},$
(1.1)

where $\bm{\mu}=(\mu_{1},...,\mu_{p})^{T}$ is the expectation, $\xi_{t}\geq 0$
is a random scalar variable, $\bm{\Sigma}_{0}$ is a $p\times p$ nonnegative
matrix, and $\mathbf{u}_{t}$ is a length-$p$ vector uniformly distributed on
the unit sphere $\mathcal{S}^{p-1}$ and independent of $\xi_{t}$. The
covariance matrix of $(\mathbf{y}_{t})$ is
$\bm{\Sigma}=E(\xi_{t}^{2})\bm{\Sigma}_{0}$. Usually, $\bm{\Sigma}_{0}$ is
referred to as the scatter matrix (fan2018large). For identification purposes,
$\bm{\Sigma}_{0}$ can be normalized to have a trace of $p$. Popular estimators
of the scatter matrix include the minimum covariance determinant (MCD) and
minimum volume ellipsoid (MVE) estimators (rousseeuw1984least;
rousseeuw1985multivariate), Maronna’s M-estimator (maronna1976robust) and
Tyler’s M-estimator (tyler1987distribution). However, these estimators can not
handle the case when $p>n$; see, for example, the discussions in
van2009minimum, zhang2016robust; zhang2016marvcenko and hubert2018minimum.

High-dimensional data such as financial returns usually exhibit a factor
structure. FFL08; FLM11; fan2013large utilize factor models coupled with
idiosyncratic sparsity to estimate the high-dimensional covariance matrix
under such settings. These studies focus on the case where data are sub-
Gaussian or sub-Exponential.

fan2018large consider elliptical factor models under which only $4$th moment
exists and propose a generic Principal Orthogonal ComplEment Thresholding
(POET) procedure to estimate the high-dimensional covariance matrix. Their
approach relies on three components, $\widehat{\bm{\Sigma}}$,
$\widehat{\bm{\Lambda}}_{K}$ and $\widehat{\bm{\Gamma}}_{K}$, which are robust
pilot estimators of the covariance matrix $\bm{\Sigma}$, its leading
eigenvalues $\bm{\Lambda}_{K}$ and eigenvectors $\bm{\Gamma}_{K}$,
respectively. The theoretical properties of the generic POET procedure
crucially depend on the consistency of the pilot components, which have sub-
Gaussian convergence rate under the $4$th moment condition:

$\displaystyle\|\widehat{\bm{\Sigma}}-\bm{\Sigma}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)},$ (1.2)
$\displaystyle\|\widehat{\bm{\Lambda}}_{K}{\bm{\Lambda}}^{-1}_{K}-\mathbf{I}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)},\quad\text{and}$ (1.3)
$\displaystyle\|\widehat{\bm{\Gamma}}_{K}-{\bm{\Gamma}}_{K}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{np}}\Bigg{)}.$ (1.4)

When only $2+\varepsilon$th moment exists, the convergence rates in equations
(1.2) and (1.3) reduce to $((\log p)/n)^{\varepsilon/(2+\varepsilon)}$.

###  1.3 Our Contributions

Infinite kurtosis and cross-sectional tail dependence are commonly seen in
financial data. Stock returns contain large and dependent jumps (e.g.,
ding2023stock), which have slowly decaying tails that suggest infinite 4th
moment (bollerslev2013jump). For data with only $2+\varepsilon$th moment, the
aforementioned covariance matrix estimation methods can not achieve the sub-
Gaussian rate. An important and challenging question is:

_Is there a covariance matrix estimator that can achieve a higher rate than
$((\log p)/n)^{\varepsilon/(2+\varepsilon)}$ or even the sub-Gaussian rate of
$\sqrt{(\log p)/n}$ when only $2+\varepsilon$th moment exists?_

We address this question in this paper. We study the robust estimation of
high-dimensional covariance matrices under elliptical factor models with only
$2+\varepsilon$th moment. It is worth mentioning that in the elliptical model
(1.1), the scatter matrix $\bm{\Sigma}_{0}$ determines the cross-sectional
dependence and plays a central role in various applications. For example, the
eigenstructure of $\bm{\Sigma}_{0}$ governs the principal component analysis
(PCA) of $(\mathbf{y}_{t})$. In financial applications, the inverse of the
scatter matrix, $\bm{\Sigma}_{0}^{-1}$, is used in constructing the minimum
variance portfolio. The scalar component, $E(\xi_{t}^{2})$, on the other hand,
is not relevant in the above applications. These observations motivate us to
estimate the two components, $E(\xi^{2}_{t})$ and $\bm{\Sigma}_{0}$,
separately.

In this paper, we propose an estimator of the scatter matrix $\bm{\Sigma}_{0}$
that achieves the sub-Gaussian convergence rate under the condition that
$(\mathbf{y}_{t})$ has only $2+\varepsilon$th moment. We summarize the
construction of our estimator below.

First, we estimate $(\xi_{t})$ and remove the effect of $(\xi_{t})$ in
$(\mathbf{y}_{t})$. ZL11 develop a self-normalization approach to remove
$(\xi_{t})$ in elliptical models. Specifically, when $p$ is large, $(\xi_{t})$
can be consistently estimated by
$(\|\mathring{\mathbf{y}}_{t}\|_{2}/\sqrt{p})$, where
$(\mathring{\mathbf{y}}_{t})$ is a properly demeaned $(\mathbf{y}_{t})$. The
consistency of $(\|\mathring{\mathbf{y}}_{t}\|_{2}/\sqrt{p})$ in estimating
$(\xi_{t})$ requires that the data do not admit a strong factor structure.
Under the elliptical factor model, such a consistency no longer holds. In this
paper, we develop an idiosyncratic projected self-normalization (IPSN)
approach that allows the existence of factors.

Next, we use the idiosyncratic-projected self-normalized data to construct the
pilot components for the generic POET estimation of $\bm{\Sigma}_{0}$. We use
the sample covariance matrix of the idiosyncratic-projected self-normalized
data and standardize it to have trace $p$. The resulting estimator is the
pilot covariance matrix estimator, $\widehat{\bm{\Sigma}}_{0}$. The pilot
leading eigenvalues and leading eigenvector estimators,
$\widehat{\bm{\Lambda}}_{0K}$ and $\widehat{\bm{\Gamma}}_{0K}$, are simply the
leading eigenvalues and eigenvectors of $\widehat{\bm{\Sigma}}_{0}$. We show
that our pilot estimators achieve the sub-Gaussian rate under only
2+$\varepsilon$th moment condition. We then apply the generic POET procedure
on the proposed IPSN pilot estimators to obtain the scatter matrix estimator,
denoted by GPOET-IPSN. We show that the GPOET-IPSN estimator achieves the sub-
Gaussian convergence rate under a conditional sparsity assumption.

Finally, to estimate the covariance matrix $\bm{\Sigma}$, we construct a
consistent estimator of $E(\xi_{t}^{2})$ based on
$(\|\mathring{\mathbf{y}}_{t}\|_{2}^{2})$. We show that the estimator achieves
the same converge rate as the oracle estimator based on the unobserved
$(\xi_{t})$ process. Combining our proposed GPOET-IPSN estimator of
$\bm{\Sigma}_{0}$ with the estimator of $E(\xi_{t}^{2})$ yields our estimator
of $\bm{\Sigma}$, which we show achieves a faster convergence rate than the
generic POET estimator in fan2018large.

The favorable properties of our proposed estimator are clearly demonstrated in
the numerical studies. Simulation studies show that the GPOET-IPSN estimator
performs robustly well under various heavy-tailedness settings and achieves a
higher estimation accuracy than the robust generic POET estimator in
fan2018large and the POET estimator in fan2013large. Empirically, we apply our
proposed estimator to construct minimum variance portfolios using S&P 500
Index constituent stocks. We find that our portfolio has a statistically
significantly lower out-of-sample risk than benchmark portfolios.

The rest of this paper is organized as follows. We present our theoretical
results in Section 2. Simulations and empirical studies are presented in
Sections 3 and 4, respectively. We conclude in Section 5. The proofs of the
theorems and propositions are presented in Section A. The proofs of the lemmas
are collected in the Supplementary Material (DZ24_supp).

We use the following notation throughout the paper. For any matrix
$\mathbf{A}=(A_{ij})$, its spectral norm is defined as
$\|\mathbf{A}\|_{2}=\max_{\|\mathbf{x}\|_{2}\leq
1}\sqrt{\mathbf{x}^{T}\mathbf{A}^{T}\mathbf{A}\mathbf{x}}$, where
$\|\mathbf{x}\|_{2}=\sqrt{\sum x_{i}^{2}}$ for any vector
$\mathbf{x}=(x_{i})$; the Frobenius norm is defined as
$\|\mathbf{A}\|_{F}=\sqrt{\sum_{i,j}A_{ij}^{2}}$; and the relative Frobenius
norm is defined as
$\|\mathbf{A}\|_{\bm{\Sigma}}=\|\bm{\Sigma}^{-1/2}\mathbf{A}\bm{\Sigma}^{-1/2}\|_{F}/\sqrt{p}$,
where $\bm{\Sigma}$ is a $p\times p$ positive-definite matrix.

##  2 Main Results

###  2.1 Settings and Assumptions

Suppose that $(\mathbf{y}_{t})$ follows the elliptical model (1.1). We make
the following assumptions.

###### Assumption 1

$(\xi_{t})$ are i.i.d. and independent with $(\mathbf{u}_{t})$. Moreover,
there exist constants $0<c<C<\infty$ such that $E(\xi_{t}^{2+\varepsilon})<C$
for some $0<\varepsilon<2$ and $P(|\xi_{t}|>c)=1$.

###### Assumption 2

There exist constants $0<c<C<\infty$ such that $c<\min_{1\leq 1\leq
p}(\bm{\Sigma}_{0})_{ii}\leq\max_{1\leq 1\leq p}(\bm{\Sigma}_{0})_{ii}<C$.

Write the eigen decomposition of $\bm{\Sigma}_{0}$ as
$\bm{\Sigma}_{0}=\mathbf{\Gamma}\bm{\Lambda}_{0}\mathbf{\Gamma}^{T}$, where
$\bm{\Lambda}_{0}=\operatorname{diag}(\lambda_{0;1},\lambda_{0;2},...,\lambda_{0;p})$,
$\lambda_{0;1}\geq\lambda_{0;2}...\geq\lambda_{0;p}$ are the eigenvalues of
$\bm{\Sigma}_{0}$, and $\mathbf{\Gamma}=(\phi_{1},...,\phi_{p})$ is the
corresponding matrix of eigenvectors. We assume that factors exist so that the
following eigen-gap condition holds.

###### Assumption 3

There exist a $K\in\mathbb{N}$ and $\delta>0$ such that $\delta
p\leq\lambda_{0;i}\leq p/\delta$, and $0<\delta\leq\lambda_{0;i}\leq 1/\delta$
for $i=K+1,...,p$.

Widely used approaches for estimating the number of factors such as those
developed in bai2002determining, onatski2010determining and ahn2013eigenvalue
require the existence of 4th moment. Under the elliptical factor model with
$2+\varepsilon$th moment, the number of factors $K$ can be consistently
estimated using the robust estimator in yu2019robust, which utilizes
eigenvalue ratios of spatial Kendall’s tau matrices. In the following
analysis, we assume that the number of factors is known.

###  2.2 Idiosyncratic-Projected Self-Normalized Estimator of Covariance
Matrix

Observe that under the model (1.1), we can write $\mathbf{y}_{t}$ as follows:

$\mathbf{y}_{t}=\bm{\mu}+\frac{\sqrt{p}\xi_{t}}{\|\mathbf{z}_{t}\|_{2}}\bm{\Sigma}^{1/2}_{0}\mathbf{z}_{t},$
(2.1)

where $\mathbf{z}_{t}\sim N(0,\mathbf{I}_{p})$, and $\xi_{t}$ and
$\mathbf{z}_{t}$ are independent. The heavy-tailedness of $(\mathbf{y}_{t})$
is solely determined by $(\xi_{t})$. Motivated by this observation, we
separately estimate $E(\xi^{2}_{t})$ and $\bm{\Sigma}_{0}$. To do so, we
develop an idiosyncratic-projected self-normalization (IPSN) approach. We
explain the construction of our proposed estimator below.

Step I. Estimate $\bm{\mu}$

Given that $(\mathbf{y}_{t})$ is heavy-tailed, in Step I, we estimate
$\bm{\mu}$ using a Huber estimator. Specifically, the estimator
$\widehat{\bm{\mu}}=(\widehat{\mu}_{i})_{1\leq i\leq p}$ solves

$\sum_{t=1}^{n}\Psi_{H}(y_{i,t}-\widehat{\mu}_{i})=0\quad\text{ for}\quad
i=1,...,p,$ (2.2)

where $\Psi_{H}$ is the Huber function:

$\Psi_{H}(x)=\min\Big{(}H,\max(-H,x)\Big{)}$ (2.3)

with $H$ a tuning parameter, which goes to infinity with $n$ and satisfies
$H=O(\sqrt{n/(\log p)})$. Numerically, the tuning parameter $H$ can be chosen
with cross-validation (sun2020adaptive).

Step II. Construct idiosyncratic-projected self-normalized variables
$(\widehat{\mathbf{X}}_{t})$

In Step II, we normalize observations to remove the effect of $(\xi_{t})$. We
use the $\ell_{2}$ norm of estimated idiosyncratic variables to conduct
normalization. Specifically, we first compute the spatial Kendall’s tau
matrix:

$\widehat{\bm{\Sigma}}_{2}=\frac{2}{n(n-1)}\sum_{i<i^{\prime}}\frac{(\mathbf{y}_{i}-\mathbf{y}_{i^{\prime}})(\mathbf{y}_{i}-\mathbf{y}_{i^{\prime}})^{T}}{\|\mathbf{y}_{i}-\mathbf{y}_{i^{\prime}}\|_{2}^{2}},$

get its leading eigenvectors, denoted by $\widehat{\bm{\Gamma}_{K}}_{ED}$, and
use them as a robust estimator of $\bm{\Gamma}_{K}=(\phi_{1},...,\phi_{K})$.
Then, we define a $(p-K)\times p$ projection matrix $\widehat{\mathbf{P}}_{I}$
as follows:

$\widehat{\mathbf{P}}_{I}^{T}=\text{\rm
Null}(\widehat{\bm{\Gamma}_{K}}_{ED}),$

that is, $\widehat{\mathbf{P}}_{I}^{T}$ is the $p\times(p-K)$ matrix that
spans the null space of $\widehat{\bm{\Gamma}_{K}}_{ED}$ and satisfies
$\widehat{\mathbf{P}}_{I}\widehat{\bm{\Gamma}_{K}}_{ED}=\mathbf{0}$ and
$\widehat{\mathbf{P}}_{I}\widehat{\mathbf{P}}_{I}^{T}=\mathbf{I}_{p-K}$.

The goal of the idiosyncratic projection is to remove the strong cross-
sectional dependence in the data, thereby enabling consistent estimation of
the scalar component $(\xi_{t})$. This step is crucial. As can be seen in
Section 3.3 below, when there is strong cross-sectional dependence, the self-
normalization approach in ZL11 no longer works in estimating the time-varying
scalar component $(\xi_{t})$.

Next, we define the idiosyncratic projected self-normalized variables
$(\widehat{\mathbf{X}}_{t})$:

$\widehat{\mathbf{X}}_{t}=\sqrt{p}\frac{\mathbf{y}_{t}-\widehat{\bm{\mu}}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}},\quad\quad
t=1,...,n.$ (2.4)

The normalization terms
$(\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2})$ are
approximately proportional to $(\xi_{t})$, hence $(\widehat{\mathbf{X}}_{t})$
are approximately proportional to $(\bm{\Sigma}_{0}^{1/2}\mathbf{u}_{t})$.

Step III. Construct pilot estimators of $\bm{\Sigma}_{0}$

In Step III, we estimate $\bm{\Sigma}_{0}$ with the following estimator:

$\widehat{\bm{\Sigma}}_{0}=\frac{\widehat{\eta}}{n}\sum_{t=1}^{n}\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T},$
(2.5)

where

$\widehat{\eta}=\frac{p}{\operatorname{tr}(\sum_{t=1}^{n}\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T}/n)}.$
(2.6)

By doing so, we essentially normalize $\mathbf{y}_{t}-\widehat{\bm{\mu}}$ with

$\widehat{\xi}_{t}=\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}}{\sqrt{\widehat{\eta}}}.$
(2.7)

By the definition of $\widehat{\eta}$ in equation (2.6), we ensure that
$\operatorname{tr}(\widehat{\bm{\Sigma}}_{0})=p$.

Write the eigen-decomposition of $\widehat{\bm{\Sigma}}_{0}$ as
$\widehat{\bm{\Sigma}}_{0}=\widehat{\mathbf{\Gamma}}\widehat{\bm{\Lambda}}_{0}\widehat{\mathbf{\Gamma}}^{T}$,
where
$\widehat{\bm{\Lambda}}_{0}=\operatorname{diag}(\widehat{\lambda}_{0;1},...,\widehat{\lambda}_{0;p})$
is the matrix of eigenvalues with
$\widehat{\lambda}_{0;1}\geq\widehat{\lambda}_{0;2}\geq...\geq\widehat{\lambda}_{0;p}$,
and $\widehat{\mathbf{\Gamma}}=(\widehat{\phi}_{1},...,\widehat{\phi}_{p})$ is
the matrix of eigenvectors. We then use the leading eigenvectors and leading
eigenvalues of $\widehat{\bm{\Sigma}}_{0}$ as the corresponding pilot
estimators:

$\widehat{\mathbf{\Gamma}}_{K}=(\widehat{\phi}_{1},...,\widehat{\phi}_{K}),\quad\text{
and
}\quad\widehat{\bm{\Lambda}}_{0K}=\operatorname{diag}(\widehat{\lambda}_{0;1},...,\widehat{\lambda}_{0;K}).$

The estimators $\widehat{\bm{\Sigma}}_{0}$, $\widehat{\bm{\Gamma}}_{K}$ and
$\widehat{\bm{\Lambda}}_{0K}$ will be used in the generic POET procedure for
estimating $\bm{\Sigma}_{0}$.

In the generic POET procedure proposed in fan2018large, the leading
eigenvectors are estimated using the spatial Kendall’s tau, and pilot
components of the covariance matrix and its leading eigenvalues are from a
Huber estimator of the covariance matrix. As these estimators are from
different sources, the resulting pilot idiosyncratic covariance matrix
estimator is not guaranteed to be semi-positive definite. Our approach, on the
other hand, does not have such an issue because our pilot estimators of the
leading eigenvalues and leading eigenvectors are from the eigen decomposition
of a same covariance matrix estimator.

Step IV. Estimate $E(\xi_{t}^{2})$

In order to estimate $\bm{\Sigma}$, we also need to estimate $E(\xi_{t}^{2})$.
We propose the following robust estimator $\widehat{E(\xi_{t}^{2})}$ that
solves

$\sum_{i=1}^{n}\Psi_{H}\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-\widehat{E(\xi_{t}^{2})}\Bigg{)}=0,$
(2.8)

where $H\asymp n^{\min(1/(1+\varepsilon/2),1/2)}$. Again, the tuning parameter
$H$ can be chosen using cross-validation.

Step V. Generic POET estimators of $\bm{\Sigma}_{0}$ and $\bm{\Sigma}$

Finally, we estimate the scatter matrix $\bm{\Sigma}_{0}$ and the covariance
matrix $\bm{\Sigma}$. To estimate $\bm{\Sigma}_{0}$, we apply the generic POET
procedure on our proposed idiosyncratic-projected self-normalized pilot
estimators, $\widehat{\bm{\Sigma}}_{0}$, $\widehat{\bm{\Lambda}}_{0K}$ and
$\widehat{\bm{\Gamma}}_{K}$ from Step III. Specifically, the estimator of
$\bm{\Sigma}_{0}$ is

$\widehat{\bm{\Sigma}}^{\tau}_{0}=\widehat{\bm{\Gamma}}_{K}\widehat{\bm{\Lambda}}_{0K}\widehat{\bm{\Gamma}}_{K}^{T}+\widehat{\bm{\Sigma}}_{0u}^{\tau},$
(2.9)

where $\widehat{\bm{\Sigma}}_{0u}^{\tau}$ is obtained by applying the adaptive
thresholding method (cai2011adaptive) to
$\widehat{\bm{\Sigma}}_{0}-\widehat{\bm{\Gamma}}_{K}\widehat{\bm{\Lambda}}_{0K}\widehat{\bm{\Gamma}}_{K}^{T}$.

To estimate $\bm{\Sigma}$, we combine $\widehat{E(\xi_{t}^{2})}$ from Step IV
with $\widehat{\bm{\Sigma}}^{\tau}_{0}$ in equation (2.9):

$\widehat{\bm{\Sigma}}^{\tau}=\widehat{E(\xi_{t}^{2})}\,\widehat{\bm{\Sigma}}^{\tau}_{0}.$
(2.10)

The corresponding idiosyncratic covariance matrix estimator is
$\widehat{\bm{\Sigma}}_{u}^{\tau}=\widehat{E(\xi_{t}^{2})}\,\widehat{\bm{\Sigma}}_{0u}^{\tau}$.
It can be seen that $\widehat{\bm{\Sigma}}^{\tau}$ is equivalent to the
generic POET estimator using the following pilot estimators:
$\widehat{\bm{\Sigma}}=:\widehat{E(\xi_{t}^{2})}\,\widehat{\bm{\Sigma}}_{0}$,
$\widehat{\bm{\Lambda}}_{K}=:\widehat{E(\xi_{t}^{2})}\,\widehat{\bm{\Lambda}}_{0K}$,
and $\widehat{\bm{\Gamma}}_{K}$.

We summarize our IPSN approach in the following algorithm.

Algorithm: Idiosyncratic-Projected Self-Normalization (IPSN) for Covariance
Matrix Estimation Input: $(\mathbf{y}_{t})_{t\leq n}$, $K$. Output:
$\widehat{\bm{\Sigma}}^{\tau}_{0}$, $\widehat{\bm{\Sigma}}^{\tau}$. Step I.
Compute $\widehat{\bm{\mu}}$ via (2.2). Step II. Compute
$(\widehat{\mathbf{X}}_{t})_{t\leq n}$ via (2.2). Step III. Compute
$\widehat{\eta}$ and $(\widehat{\xi}_{t})_{t\leq n}$ via (2.6) and (2.7),
respectively. Compute $\widehat{\bm{\Sigma}}_{0}$ via (2.5), its first $K$
eigenvalues $\widehat{\bm{\Lambda}}_{0K}$, and the corresponding eigenvectors
$\widehat{\mathbf{\Gamma}}_{K}$. Step IV. Compute $\widehat{E(\xi_{t}^{2})}$
via (2.8). Step V. Compute $\widehat{\bm{\Sigma}}^{\tau}_{0}$ via (2.9), and
$\widehat{\bm{\Sigma}}^{\tau}$ via (2.10).

###  2.3 Theoretical Properties

In this subsection, we present the asymptotic properties of our proposed
estimators.

####  2.3.1 Convergence Rates of Idiosyncratic-Projected Self-Normalized Pilot
Estimators

We start with the properties of the idiosyncratic-projected self-
normalization. The key to the success of our approach lies in approximating
$(\xi_{t})$ with $(\widehat{\xi}_{t})$ defined in equation (2.7). Intuitively,
if $(\widehat{\xi}_{t})$ is close to $(\xi_{t})$, then
$((\mathbf{y}_{t}-\widehat{\bm{\mu}})/\widehat{\xi}_{t})$ is close to
$\sqrt{p}(\bm{\Sigma}_{0}^{1/2}\mathbf{u}_{t})$, which is absent of the heavy-
tailed scalar process $(\xi_{t})$. The next proposition gives the properties
of $~{}(\widehat{\xi}_{t})$.

###### Proposition 1

Under Assumptions 1–3, if in addition, $p,n\to\infty$ and satisfy
$(\log(p))^{2+\gamma}=o(n)$ for some $\gamma>0$ and $\log(n)=O(\log(p))$, then
there exists a constant $c_{1}>0$ such that

$P\Bigg{(}\min_{1\leq t\leq n}\widehat{\xi}_{t}>c_{1}\Bigg{)}\to 1,$ (2.11)

and

$\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\xi_{t}^{2}}{\widehat{\xi}_{t}^{2}}-1\Big{)}^{2}=O_{p}\Big{(}\frac{\log
p}{p}+\frac{\log p}{n}\Big{)}.$ (2.12)

Equation (2.11) guarantees that, with probability approaching one,
$(\widehat{\xi}_{t})$ is bounded away from zero, hence the normalization based
on $(\widehat{\xi}_{t})$ is well-behaved. Furthermore, equation (2.12)
guarantees the consistency of $(\widehat{\xi}_{t})$ in estimating
$({\xi}_{t})$. We need both $p,n\to\infty$ to estimate $(\xi_{t})$
consistently, and $p$ can be much larger than $n$.

The next theorem gives the convergence rate of the IPSN pilot estimators,
$\widehat{\bm{\Sigma}}_{0}$, $\widehat{\bm{\Lambda}}_{0K}$, and
$\widehat{\bm{\Gamma}}_{K}$.

###### Theorem 1

Under Assumptions 1–3, if in addition, $p,n\to\infty$ and satisfy
$(\log(p))^{2+\gamma}=o(n)$ for some $\gamma>0$ and $\log(n)=O(\log(p))$, then

$\displaystyle\|\widehat{\bm{\Sigma}}_{0}-\bm{\Sigma}_{0}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p}}+\sqrt{\frac{\log p}{n}}\Bigg{)},$ (2.13)
$\displaystyle\|\widehat{\bm{\Lambda}}_{0K}\bm{\Lambda}_{0K}^{-1}-\mathbf{I}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p}}+\sqrt{\frac{\log p}{n}}\Bigg{)},\quad\text{and}$
$\displaystyle\|\widehat{\bm{\Gamma}}_{K}-\bm{\Gamma}_{K}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p^{2}}}+\sqrt{\frac{\log p}{np}}\Bigg{)}.$

Theorem 1 states that the pilot estimators of the scatter matrix achieve the
sub-Gaussian convergence rates under $2+\varepsilon$th moment assumption. By
contrast, when only $2+\varepsilon$th moment exists, the pilot estimator
$\widehat{\bm{\Sigma}}$ in fan2018large has a convergence rate of only
$O_{p}(((\log p)/n)^{\varepsilon/(2+\varepsilon)})$.

####  2.3.2 Convergence Rates of Generic POET Estimators based on IPSN

Under the factor model specified in Assumption 3, we can write
$\bm{\Sigma}_{0}=\mathbf{B}_{0}\mathbf{B}_{0}^{T}+\bm{\Sigma}_{0u}$, where
$\mathbf{B}_{0}$ is a $p\times K$ matrix, and $\bm{\Sigma}_{0u}$ is the
idiosyncratic component of the scatter matrix $\bm{\Sigma}_{0}$.

We consider factor models with conditional sparsity. Specifically, we assume
that $\bm{\Sigma}_{0u}$ belongs to the following class of sparse matrices: for
some $q\in[0,1]$, $c>0$, and $s_{0}(p)<\infty$,

$\displaystyle\mathcal{U}_{q}\big{(}s_{0}(p)\big{)}$ (2.14) $\displaystyle=$
$\displaystyle\Big{\\{}\bm{\Sigma}:\bm{\Sigma}\text{ is positive semi-
definite},\,\max_{i}\sigma_{ii}\leq c,\text{and
}\max_{i}\sum_{j=1}\sigma_{ij}^{q}\leq s_{0}(p)\Big{\\}}.$

###### Assumption 4

$\bm{\Sigma}_{0u}\in\mathcal{U}_{q}\big{(}s_{0}(p)\big{)}$ for some
$q\in[0,1]$.

Denote $\omega_{n}=\sqrt{(\log p)/p}+\sqrt{(\log p)/n}$. The next theorem
gives the convergence rate of our proposed generic POET estimator of the
scatter matrix.

###### Theorem 2

Under the assumptions of Theorem 1 and Assumption 4, if in addition,
$s_{0}(p)\omega_{n}^{1-q}=o(1)$, then

$\displaystyle\|(\widehat{\bm{\Sigma}}_{0}^{\tau})^{-1}-(\bm{\Sigma}_{0})^{-1}\|_{2}=O_{p}\Big{(}s_{0}(p)\omega_{n}^{1-q}\Big{)},$
(2.15)
$\displaystyle\|\widehat{\bm{\Sigma}}_{0}^{\tau}-\bm{\Sigma}_{0}\|_{\bm{\Sigma}_{0}}=O_{p}\big{(}\sqrt{p}\omega_{n}^{2}+s_{0}(p)\omega_{n}^{1-q}\Big{)},\quad\text{and}$
$\displaystyle\|\widehat{\bm{\Sigma}}_{0u}^{\tau}-\bm{\Sigma}_{0u}\|_{2}=O_{p}\Big{(}s_{0}(p)\omega_{n}^{1-q}\Big{)}=\|(\widehat{\bm{\Sigma}}_{0u}^{\tau})^{-1}-\bm{\Sigma}_{0u}^{-1}\|_{2}.$

Theorem 2 asserts that in estimating of the scatter matrix, its idiosyncratic
components and their inverse matrices, under only $2+\varepsilon$th moment
condition, our generic POET estimator achieves the same convergence rate as
the sub-Gaussian case; see Theorems 3.1 and 3.2 in fan2013large.

Finally, about the estimation of the covariance matrix, recall that we combine
the estimators $\widehat{E(\xi_{t}^{2})}$ and
$\widehat{\bm{\Sigma}}_{0}^{\tau}$ to estimate $\bm{\Sigma}$, and combine
$\widehat{E(\xi_{t}^{2})}$ and $\widehat{\bm{\Sigma}}_{0u}^{\tau}$ to estimate
$\bm{\Sigma}_{u}$. The next proposition gives the convergence rate of
$\widehat{E(\xi_{t}^{2})}$ defined in equation (2.8) in estimating
$E(\xi_{t}^{2})$.

###### Proposition 2

Under Assumptions 1–3, if in addition, $n\to\infty$ and $\log(p)=O(n)$, then

$|\widehat{E(\xi_{t}^{2})}-E(\xi_{t}^{2})|=O_{p}\Bigg{(}\frac{1}{n^{\varepsilon/(2+\varepsilon)}}+\sqrt{\frac{\log
p}{n}}\Bigg{)}.$

The error term $1/n^{\varepsilon/(2+\varepsilon)}$ is from robust Huber
estimation based on the infeasible series $(\xi_{t})$, and the error term
$\sqrt{\log p/n}$ comes from the estimation error of $\widehat{\bm{\mu}}$.
When $\log p=o(n^{(2-\varepsilon)/(2+\varepsilon)})$, the convergence rate of
$\widehat{E(\xi_{t}^{2})}$ is the same as the robust Huber estimator applied
to the unobserved series $(\xi_{t})$.

The next theorem gives the convergence rates of $\widehat{\bm{\Sigma}}^{\tau}$
and $\widehat{\bm{\Sigma}}_{u}^{\tau}$ in estimating the covariance matrix and
the idiosyncratic covariance matrix, respectively.

###### Theorem 3

Under the assumptions of Theorem 1 and Assumption 4, if in addition,
$s_{0}(p)\omega_{n}^{1-q}=o(1)$, then

$\displaystyle\|(\widehat{\bm{\Sigma}}^{\tau})^{-1}-\bm{\Sigma}^{-1}\|_{2}=O_{p}\Bigg{(}s_{0}(p)\omega_{n}^{1-q}+\frac{1}{n^{\varepsilon/(2+\varepsilon)}}\Bigg{)},$
(2.16)
$\displaystyle\|\widehat{\bm{\Sigma}}^{\tau}-\bm{\Sigma}\|_{\bm{\Sigma}}=O_{p}\Bigg{(}\sqrt{p}\omega_{n}^{2}+s_{0}(p)\omega_{n}^{1-q}+\frac{1}{n^{\varepsilon/(2+\varepsilon)}}\Bigg{)},\quad\text{and}$
$\displaystyle\|\widehat{\bm{\Sigma}}_{u}^{\tau}-\bm{\Sigma}_{u}\|_{2}=O_{p}\Bigg{(}s_{0}(p)\omega_{n}^{1-q}+\frac{1}{n^{\varepsilon/(2+\varepsilon)}}\Bigg{)}=\|(\widehat{\bm{\Sigma}}_{u}^{\tau})^{-1}-\bm{\Sigma}_{u}^{-1}\|_{2}.$

One can show that when only $2+\varepsilon$th moment exists, the generic POET
estimator proposed in fan2018large has a lower convergence rate. Specifically,
it has a convergence rate with the term $\omega^{2}_{n}$ replaced by
$\widetilde{\omega}^{2}_{n}$ and the term $\omega_{n}^{1-q}$ replaced by
$\widetilde{\omega}_{n}^{2\varepsilon(1-q)/(2+\varepsilon)}$ for
$\widetilde{\omega}_{n}=\sqrt{\log p/n}+1/\sqrt{p}$.

##  3 Simulation Studies

###  3.1 Simulation Setting

We generate data from the following model. We first simulate
$\mathbf{B}=(b_{ik})_{1\leq i\leq N,1\leq k\leq K}$, where
$b_{ik}\underset{\text{i.i.d.}}{\sim}N(0,s_{i})$, $K=3$, $s_{1}=1$,
$s_{2}=0.75^{2}$, and $s_{3}=0.5^{2}$. We then set
$\bm{\Sigma}_{0}=p(\mathbf{B}\mathbf{B}^{T}+\mathbf{I}_{p})/\operatorname{tr}((\mathbf{B}\mathbf{B}^{T}+\mathbf{I}_{p}))$.

We generate the scalar process $(\xi_{t})$ from a Pareto distribution with
$P(\xi_{t}>x)=(x_{m}/x)^{\alpha}$ for $x>x_{m}$, where $x_{m}$ is chosen such
that $E(\xi_{t}^{2})=1$. We set $\alpha=2+\varepsilon$ with $\varepsilon=2$ or
$0.2$. When $\varepsilon=2$, $(\xi_{t})$ has finite moments with order below
4, and when $\varepsilon=0.2$, $(\xi_{t})$ does not have a finite third
moment.

We then generate $(\mathbf{y}_{t})$ from
$\mathbf{y}_{t}=\xi_{t}\sqrt{p}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}/\|\mathbf{z}_{t}\|_{2}$,
$t=1,...,n$, where
$\mathbf{z}_{t}\underset{\text{i.i.d.}}{\sim}N(0,\mathbf{I})$. Note that the
estimators of the scatter matrices and covariance matrices are location
invariant, hence without lost of generality, we set $\bm{\mu}=\mathbf{0}$. The
dimension and the sample size $p$ and $n$ are set to be $p=500$ and $n=250$,
respectively.

###  3.2 Covariance Matrix Estimators and Evaluation Metrics

We evaluate the performance of the following covariance matrix estimators: the
generic POET estimators based on our idiosyncratic-projected self-
normalization approach, denoted by GPOET-IPSN, the robust generic POET
estimators in fan2018large, denoted by GPOET-FLW, and the original POET
estimator based on the sample covariance matrix (fan2013large), denoted by
POET-S.

We use the same evaluation metrics as in fan2018large. Specifically, errors in
estimating $\bm{\Sigma}_{0}$ and $\bm{\Sigma}$ are measured by the difference
in the relative Frobenius norm:

$\displaystyle\|\widehat{\bm{\Sigma}}^{\tau}_{0}-{\bm{\Sigma}}_{0}\|_{{\bm{\Sigma}}_{0}}=\frac{1}{\sqrt{p}}\|{\bm{\Sigma}}_{0}^{-1/2}(\widehat{\bm{\Sigma}}^{\tau}_{0}-{\bm{\Sigma}}_{0}){\bm{\Sigma}}_{0}^{-1/2}\|_{F},$
$\displaystyle\|\widehat{\bm{\Sigma}}^{\tau}-{\bm{\Sigma}}\|_{{\bm{\Sigma}}}=\frac{1}{\sqrt{p}}\|{\bm{\Sigma}}^{-1/2}(\widehat{\bm{\Sigma}}^{\tau}-{\bm{\Sigma}})\widehat{\bm{\Sigma}}^{-1/2}\|_{F}.$

Errors in estimating the matrices $\bm{\Sigma}_{0u}$, $\bm{\Sigma}_{u}$,
$\bm{\Sigma}_{0}^{-1}$, and $\bm{\Sigma}^{-1}$ are measured by the difference
in spectral norm.

###  3.3 Simulation Results

First, to illustrate how well our approach works in estimating the scalar
process $(\xi_{t})$, we compute the estimator $(\widehat{\xi}_{t})$ from our
IPSN approach in equation (2.7) and obtain the ratio
$(\widehat{\xi}_{t}/\xi_{t})$ from one realization. We compare the results
with the estimator in ZL11 (denoted by ZL11). The results are shown in Figure
1.

Figure 1: The ratios between $(\widehat{\xi}_{t})$ and $(\xi_{t})$ from one
simulation. We compare our proposed IPSN estimator with the estimator (ZL11)
in ZL11. The process $(\xi_{t})$ is generated from Pareto distribution
$P(\xi_{t}>x)=(x_{m}/x)^{\alpha}$ with the tail parameter $\alpha=2.2$.

We see from Figure 1 that the ratios between our estimated and the true
$(\xi_{t})$’s are remarkably close to one for all $t$. In contrast, due to the
factor structure in the data, the estimator in ZL11 is not consistent,
reflected in that the ratios between the estimated $(\xi_{t})$ process and the
true ones can be far away from one.

Next, we summarize the performance of the pilot estimators. We include our
IPSN pilot estimators and the following benchmark methods for comparison: the
robust estimators of fan2018large (FLW), as well as the sample covariance
matrix with its leading eigenvectors and eigenvalues (SAMPLE). For the
benchmark methods, we get the pilot estimators for $\mathbf{\Sigma}_{0}$ and
$\mathbf{\Lambda}_{0K}$ by normalizing $\widehat{\bm{\Sigma}}$ and
$\widehat{\mathbf{\Lambda}}_{K}$ with
$p/\operatorname{tr}(\widehat{\bm{\Sigma}})$. For the pilot estimators of
$\widehat{\bm{\Sigma}}_{0}$ and $\widehat{\bm{\Sigma}}$, we compute the error
in the maximum norm. For the estimators of $\mathbf{\Gamma}_{K}$, we evaluate
the maximum error scaled by $\sqrt{p}$,
$(\|\widehat{\mathbf{\Gamma}}_{K}-\mathbf{\Gamma}_{K}\|_{\max})\sqrt{p}.$ For
the estimators of $\mathbf{\Lambda}_{0,K}$ and $\mathbf{\Lambda}_{K}$, we
report the errors in ratios:
$\|\widehat{\mathbf{\Lambda}}_{0,K}\mathbf{\Lambda}_{0K}^{-1}-\mathbf{I}\|_{\max},\quad\text{and}\quad\|\widehat{\mathbf{\Lambda}}_{K}\mathbf{\Lambda}_{K}^{-1}-\mathbf{I}\|_{\max}.$
In Table 1, we report the performance of the pilot estimators from 100
replications.

Table 1: Performance of the pilot estimators. The heavy-tailedness index is set to be $\varepsilon=2$ or $0.2$. The dimension and the sample size are $p=500$ and $n=250$, respectively. Reported values are mean and standard deviations (in parentheses) from 100 replications for the error in maximum norm in estimating $\bm{\Sigma}_{0}$, $\bm{\Sigma}$, and ($\sqrt{p}\times$) $\bm{\Gamma}_{K}$, and the maximum norm error in ratios in estimating $\bm{\Lambda}_{0K}$ and $\bm{\Lambda}_{K}$. We include our IPSN and the following benchmark methods: the robust method proposed in fan2018large (FLW) and the sample covariance matrix (SAMPLE).  | $\bm{\Sigma}_{0}$ | $\bm{\Sigma}$ | $\bm{\Lambda}_{0K}$ | $\bm{\Lambda}_{K}$ | $\bm{\Gamma}_{K}$
---|---|---|---|---|---
$\varepsilon=2$ |  |  |  |  |
IPSN | 0.482 | 0.676 | 0.100 | 0.165 | 0.592
| (0.108) | (0.273) | (0.048) | (0.099) | (0.177)
FLW | 0.706 | 0.921 | 0.111 | 0.181 | 0.620
| (0.392) | (0.760) | (0.053) | (0.122) | (0.178)
SAMPLE | 0.728 | 0.923 | 0.135 | 0.204 | 0.872
| (0.394) | (0.768) | (0.076) | (0.147) | (0.488)
$\varepsilon=0.2$ |  |  |  |  |
IPSN | 0.485 | 2.105 | 0.100 | 0.623 | 0.598
| (0.109) | (2.706) | (0.048) | (0.886) | (0.182)
FLW | 1.820 | 3.371 | 0.141 | 0.656 | 0.641
| (1.289) | (9.982) | (0.062) | (0.891) | (0.185)
SAMPLE | 1.913 | 3.406 | 0.326 | 0.882 | 2.135
| (1.313) | (9.987) | (0.248) | (2.350) | (1.147)

We see from Table 1 that our IPSN estimators outperform the other estimators
in all measures. The advantage of IPSN over FLW and SAMPLE is more salient in
the estimation of the scatter matrices and with increasing heavy-tailedness in
the data. In particular, for the pilot estimators of the scatter matrix and
its leading eigenvectors and leading eigenvalues, when $\varepsilon$ changes
from 2 to 0.2, that is, when data get more heavy-tailed, the estimation errors
for IPSN remain almost the same. By contrast, the errors in the FLW and SAMPLE
estimators increase substantially. For pilot estimators of the covariance
matrix, the error of our IPSN estimator grows with the heavy-tailedness of the
model, but it still significantly outperforms the other methods.

Finally, we summarize the performance of the generic POET estimators in Table
2.

Table 2: Performance of generic POET estimators. The heavy-tailedness index is set to be $\varepsilon=2$ or $0.2$. The dimension and the sample size are $p=500$ and $n=250$, respectively. Reported values are mean and standardized deviation (in parentheses) from 100 replications for the error in relative Frobenius norm in estimating $\bm{\Sigma}_{0}$ and $\bm{\Sigma}$, and the error in $\ell_{2}$ norm in estimating $\bm{\Sigma}_{0u}$, $\bm{\Sigma}_{0u}$, $\bm{\Sigma}_{0}^{-1}$ and $\bm{\Sigma}^{-1}$. We include our GPOET-IPSN and the following benchmark methods: the robust method GPOET-FLW and the original POET based on sample covariance matrix (POET-S).  | $\bm{\Sigma}_{0}$ | $\bm{\Sigma}$ | $\bm{\Sigma}_{0u}$ | $\bm{\Sigma}_{u}$ | $\bm{\Sigma}_{0}^{-1}$ | $\bm{\Sigma}^{-1}$ | $\bm{\Sigma}_{0u}^{-1}$ | $\bm{\Sigma}_{u}^{-1}$
---|---|---|---|---|---|---|---|---
$\varepsilon=2$ |  |  |  |  |  |  |  |
GPOET-IPSN | 0.241 | 0.091 | 0.107 | 0.122 | 1.023 | 2.894 | 1.029 | 1.161
| (0.008) | (0.020) | (0.012) | (0.034) | (0.209) | (0.396) | (0.210) | (0.301)
GPOET-FLW | 0.328 | 0.127 | 0.545 | 0.632 | 1.533 | 2.995 | 1.542 | 1.621
| (0.103) | (0.085) | (0.768) | (1.360) | (0.354) | (0.290) | (0.357) | (0.318)
POET-S | 0.479 | 0.186 | 0.206 | 0.224 | 1.423 | 3.235 | 1.430 | 1.494
| (0.483) | (0.278) | (0.239) | (0.318) | (0.437) | (0.284) | (0.439) | (0.289)
$\varepsilon=0.2$ |  |  |  |  |  |  |  |
GPOET-IPSN | 0.242 | 0.219 | 0.107 | 0.266 | 1.023 | 7.155 | 1.028 | 5.448
| (0.009) | (0.289) | (0.013) | (0.371) | (0.210) | (2.567) | (0.211) | (2.417)
GPOET-FLW | 0.707 | 0.375 | 2.245 | 3.098 | 4.822 | 11.321 | 6.839 | 12.611
| (0.394) | (1.218) | (2.428) | (15.229) | (10.292) | (10.969) | (20.067) | (27.690)
POET-S | 3.331 | 1.381 | 0.549 | 0.412 | 4.925 | 11.616 | 4.972 | 9.888
| (4.122) | (5.096) | (0.567) | (0.434) | (11.121) | (2.812) | (11.203) | (2.996)

We see from Table 2 that our GPOET-IPSN estimators deliver the lowest
estimation error in all measures. Compared to GPOET-FLW or POET-S, the
advantage of GPOET-IPSN is especially evident for the estimation of
$\bm{\Sigma}_{0}$, $\bm{\Sigma}_{0u}$, their inverse matrices, and for the
more heavy-tailed setting when $\varepsilon=0.2$. When estimating the scatter
matrix, its idiosyncratic component and their inverses, the performance of
GPOET-IPSN is robust for different heavy-tailedness conditions. When
$\varepsilon$ decreases from 2 to 0.2, the data become more heavy-tailed, but
the estimation errors of GPOET-IPSN remain largely the same. By contrast, the
errors of the GPOET-FLW and POET-S estimators increase sharply in the more
heavy-tailed setting. In the estimation of $\bm{\Sigma}$ and
$\bm{\Sigma}^{-1}$, GPOET-IPSN also has lower estimation errors than the
benchmark methods.

In summary, the simulation results validate the theoretical properties of our
approach and demonstrate its substantial advantage in estimating the scatter
matrix and its inverse matrix. Such advantages are particularly valuable in
applications where only the scatter matrix matters while the scalar component
does not, as we will show in the empirical studies below.

##  4 Empirical Studies

###  4.1 Global Minimum Variance Portfolio Optimization

The estimation of the minimum variance portfolio (MVP) in the high-dimensional
setting has drawn considerable attention in recent years. It is widely used to
measure the empirical performance of covariance matrix estimators
(fan2012vast; ledoit2017nonlinear; DLZ19). In this subsection, we conduct
empirical analysis to evaluate the performance of our proposed estimator in
the MVP optimization.

The MVP is the portfolio that achieves the minimum variance among all
portfolios with the constraint that the summation of weights equals one.
Specifically, for $p$ asset returns with a covariance matrix of $\bm{\Sigma}$,
the MVP solves

$\mathop{\text{min}}\mathbf{w}^{T}\bm{\Sigma}\mathbf{w},\text{ subject to
}\mathbf{w}^{T}\mathbf{1}=1,$

where $\mathbf{w}$ is a length $p$ vector of portfolio weights, and
$\mathbf{1}$ is a length $p$ vector of ones. The solution of the MVP is

$\mathbf{w}^{*}=\frac{1}{\mathbf{1}^{T}\bm{\Sigma}^{-1}\mathbf{1}}\bm{\Sigma}^{-1}\mathbf{1}.$

It is straightforward to see that under elliptical models,

$\mathbf{w}^{*}=\frac{1}{\mathbf{1}^{T}\bm{\Sigma}_{0}^{-1}\mathbf{1}}\bm{\Sigma}_{0}^{-1}\mathbf{1}.$

In other words, using ${\bm{\Sigma}}^{-1}$ and ${\bm{\Sigma}}_{0}^{-1}$ are
equivalent in constructing the MVP. Using an estimator of the inverse of the
scatter matrix or the covariance matrix, $\widehat{\bm{\Sigma}}_{0}^{-1}$ or
$\widehat{\bm{\Sigma}}^{-1}$, we estimate the MVP by

$\widehat{\mathbf{w}}^{*}=\frac{1}{\mathbf{1}^{T}\widehat{\bm{\Sigma}}_{0}^{-1}\mathbf{1}}\widehat{\bm{\Sigma}}_{0}^{-1}\mathbf{1},\quad\text{or}\quad\widehat{\mathbf{w}}^{*}=\frac{1}{\mathbf{1}^{T}\widehat{\bm{\Sigma}}^{-1}\mathbf{1}}\widehat{\bm{\Sigma}}^{-1}\mathbf{1}.$
(4.1)

###  4.2 Data and Compared Methods

We use the daily returns of S&P 500 Index constituent stocks between January
1995 and December 2023 to construct MVPs. At the beginning of each month, we
use the historical returns of the stocks that stayed in the S&P 500 Index for
the past five years to estimate the portfolio weights. We evaluate the risk of
the portfolios based on their out-of-sample daily returns from January 2000 to
December 2023.

We estimate the MVP using equation (4.1) with the following estimators: our
GPOET estimator of the scatter matrix based on idiosyncratic-projected self-
normalization, GPOET-IPSN, the GPOET-FLW estimator, and the original POET
estimator based on the sample covariance matrix, POET-S. When performing the
generic POET procedure, the tuning parameter is chosen by cross-validation
with the criterion of minimizing the out-of-sample MVP risk. The number of
factors is estimated by the method in yu2019robust. We also include the equal-
weight portfolio $(1/p,...,1/p)^{T}$ as a benchmark portfolio, denoted by EW.

###  4.3 Out-of-Sample Performance

We use the standard deviation of out-of-sample returns to evaluate the out-of-
sample risk of portfolios. In addition, we perform the following test to
evaluate the statistical significance of the differences in risks between our
GPOET-IPSN method and benchmark portfolios:

$H_{0}:\sigma\geq\sigma_{0}\quad\text{ vs. }\quad H_{1}:\sigma<\sigma_{0},$
(4.2)

where $\sigma$ denotes the standard deviation of the GPOET-IPSN portfolio and
$\sigma_{0}$ denotes the standard deviation of a benchmark portfolio. We adopt
the test in LW11 and use the heteroskedasticity-autocorrelation-consistent
(HAC) standard deviation estimator of the test statistic therein.

Table 3 reports the risks of various portfolios.

Table 3: Out-of-sample risks of the portfolios. We report the annualized standard deviations. We construct the MVP portfolio using equation (4.1) with the covariance matrix estimated using generic POET based on our idiosyncratic-projected self-normalized method (GPOET-IPSN), the robust method GPOET-FLW, and the sample covariance matrix (POET-S). We also include the equal weight (EW) portfolio as a benchmark. The out-of-sample period is from January 2000 to December 2023. The risk is evaluated using out-of-sample daily returns of the portfolios. For all methods, the portfolios weights are estimated based on prior five-year daily excess returns of S&P 500 Index constituent stocks, reestimated monthly. The symbols *, ** and *** indicate statistical significance for testing (4.2) at 5%, 1% and 0.1% levels, respectively.  | Out-of-sample portfolio risk
---|---
Period | 2000–2023 | 2000 | 2001 | 2002 | 2003
GPOET-IPSN | 0.121 | 0.154 | 0.116 | 0.140 | 0.090
GPOET-FLW | 0.150 *** | 0.168 ** | 0.167 *** | 0.173 *** | 0.121 ***
POET-S | 0.150 *** | 0.152 | 0.121 | 0.140 | 0.095 **
EW | 0.217 *** | 0.192 *** | 0.208 *** | 0.269 *** | 0.179 ***
Period | 2004 | 2005 | 2006 | 2007 | 2008
GPOET-IPSN | 0.081 | 0.084 | 0.068 | 0.081 | 0.208
GPOET-FLW | 0.095 *** | 0.100 *** | 0.089 *** | 0.100 *** | 0.224
POET-S | 0.091 *** | 0.091 *** | 0.074** | 0.087 *** | 0.188
EW | 0.128 *** | 0.115 *** | 0.113 *** | 0.166 *** | 0.460 ***
Period | 2009 | 2010 | 2011 | 2012 | 2013
GPOET-IPSN | 0.133 | 0.083 | 0.103 | 0.082 | 0.084
GPOET-FLW * | 0.228 *** | 0.094 ** | 0.108 | 0.073 | 0.081
POET-S | 0.189 *** | 0.089 * | 0.114 ** | 0.084 | 0.089
EW | 0.357 *** | 0.210 *** | 0.271 *** | 0.148*** | 0.123 ***
Period | 2014 | 2015 | 2016 | 2017 | 2018
GPOET-IPSN | 0.076 | 0.096 | 0.099 | 0.072 | 0.105
GPOET-FLW | 0.081 | 0.114 *** | 0.128 *** | 0.076 | 0.116 *
POET-S | 0.089 | 0.120 *** | 0.116 ** | 0.073 | 0.119 ***
EW | 0.119 *** | 0.158 *** | 0.157 *** | 0.077 | 0.161 ***
Period | 2019 | 2020 | 2021 | 2022 | 2023
GPOET-IPSN | 0.097 | 0.269 | 0.106 | 0.141 | 0.125
GPOET-FLW | 0.131 *** | 0.317 ** | 0.184 *** | 0.179 ** | 0.176 ***
POET-S | 0.101 | 0.413 *** | 0.194 *** | 0.226 *** | 0.168 ***
EW | 0.137 *** | 0.404 *** | 0.144 *** | 0.239 *** | 0.152 ***

We see from Table 3 that our GPOET-IPSN portfolio achieves a substantially
lower risk than benchmark portfolios for the period between 2000 and 2023. The
statistical test results suggest that the differences in the risk between our
portfolio and benchmark portfolios are all statistically significantly
negative. In addition, when checking the risk performance in different years,
we see that the GPOET-IPSN portfolio performs robustly well. It outperforms
the benchmark portfolios by achieving the lowest risk in the majority of the
years (20 out of 24 years). The statistical test results show that the
difference between the risks of our GPOET-IPSN portfolio and other portfolios
is statistically significantly negative for most of the years.

##  5 Conclusion

We propose an idiosyncratic-projected self-normalization (IPSN) approach to
estimate the scatter matrix under high-dimensional elliptical factor models.
We show that our estimator achieves the sub-Gaussian convergence rate under
only $2+\varepsilon$th moment condition, a property that can not be achieved
by existing approaches. Moreover, we develop POET estimators of the scatter
matrix and the covariance matrix under a conditional sparsity assumption, and
show that they have higher convergence rates than the generic POET estimator
in fan2018large. Numerical studies demonstrate the clear advantages of our
proposed estimators over various benchmark methods.

##  6 Proofs

In this section we prove Theorems 1–3 and Propositions 1–2. We first give the
following lemmas. The proofs of Lemmas 4–8 are given in the Supplementary
Material (DZ24_supp).

###  6.1 Lemmas

###### Lemma 1

[Corollary 1 of avella2018robust] Suppose that $\delta\in(0,1)$,
$\varepsilon\in(0,1]$, $n>12\log(2\delta^{-1})$, and $Z_{1}$, …, $Z_{n}$ are
i.i.d. random variables with mean $g$ and bounded $1+\varepsilon$th moment,
i.e. $E(|Z_{1}-g|^{1+\varepsilon})=v<\infty$. Take
$H=\big{(}vn/\log(2\delta^{-1})\big{)}^{1/(1+\varepsilon)}$. Let
$\widehat{g}_{H}$ satisfies that
$\sum_{t=1}^{n}\Psi_{H}(Z_{t}-\widehat{g}_{H})=0$ with $\Psi_{H}$ being the
Huber function defined in equation (2.3). Then with probability at least
$1-2\delta$,

$|\widehat{g}_{H}-g|\leq
Cv^{1/(1+\varepsilon)}\Big{(}\frac{\log(2\delta^{-1})}{n}\Big{)}^{\varepsilon/(1+\varepsilon)}.$

###### Lemma 2

[Lemma 9 of ZL11] Suppose that $\mathbf{X}=(x_{1},...,x_{p})^{T}$, where
$x_{i}$’s are i.i.d. random variables such that $E(x_{1})=0$, $E(x_{1}^{2})=1$
and $E(x_{1}^{2k})<\infty$ for some $2\leq k<\infty$. Then there exists
$c_{k}>0$, depending only on $k$, $E(x_{1}^{4})$ and $E(x_{1}^{2k})$ such that
for any $p\times p$ nonrandom matrix $\mathbf{A}$,

$E\Big{(}|\mathbf{X}^{T}\mathbf{A}\mathbf{X}-\operatorname{tr}(\mathbf{A})|^{2k}\Big{)}\leq
c_{k}\Big{(}\operatorname{tr}(\mathbf{A}\mathbf{A}^{T})\Big{)}^{k}\leq
c_{k}p^{k}\|\mathbf{A}\|_{2}^{2k}.$

###### Lemma 3

[Lemma 3 of chatterjee2013assumptionless] Suppose that $x_{i}\sim
N(0,\sigma_{i}^{2})$, $i=1,...,m$. The $x_{i}$’s need not be independent. Let
$L:=\max_{1\leq i\leq m}\sigma_{i}$. Then

$E\Big{(}\max_{1\leq i\leq m}|x_{i}|\Big{)}\leq L\sqrt{2\log(2m)}.$

###### Lemma 4

Under Assumptions 1 and 2, if in addition $\log(p)=O(n)$, then

$\|\widehat{\bm{\mu}}-\bm{\mu}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)}.$

###### Lemma 5

Suppose that $\mathbf{z}_{t}=(z_{1},...,z_{p})^{T}\sim{N}(0,\mathbf{I})$. For
any $p\times p$ positive semi-definite matrix $\mathbf{\Sigma}$, we have

$E\Big{(}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}\mathbf{z}_{t}}{\mathbf{z}^{T}_{t}\mathbf{z}_{t}}\Big{)}=\frac{\operatorname{tr}(\bm{\Sigma})}{p},$
(6.1)

and

$E\Bigg{(}\frac{(\mathbf{z}_{t}^{T}\bm{\Sigma}\mathbf{z}_{t})^{2}}{(\mathbf{z}^{T}_{t}\mathbf{z}_{t})^{2}}\Bigg{)}\leq\frac{\operatorname{tr}(\bm{\Sigma}^{2})}{p}+\frac{\big{(}\operatorname{tr}(\bm{\Sigma})\big{)}^{2}}{p(p-1)}.$
(6.2)

Define

$\widetilde{\mathbf{X}}_{t}=\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t},\quad\text{
and
}\quad\widetilde{\mathbf{S}}_{0}=\frac{1}{n}\sum_{t=1}^{n}\widetilde{\mathbf{X}}_{t}\widetilde{\mathbf{X}}_{t}^{T}.$

Let
$\widetilde{\bm{\Lambda}}_{0K}=\operatorname{diag}(\widetilde{\lambda}_{0;1},...,\widetilde{\lambda}_{0;K})$
with $\widetilde{\lambda}_{0;1}\geq...\geq\widetilde{\lambda}_{0;K}$ be the
leading eigenvalues and
$\widetilde{\bm{\Gamma}}_{K}=(\widehat{\phi}_{1},...,\widehat{\phi}_{K})$ the
corresponding eigenvectors of $\widetilde{\mathbf{S}}_{0}$.

###### Lemma 6

Under Assumptions 1–3, if in addition, $\log(p)=o(n)$, then

$\|\widetilde{\mathbf{S}}_{0}-\bm{\Sigma}_{0}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)},$ (6.3)
$\|\widetilde{\bm{\Lambda}}_{0K}(\bm{\Lambda}_{0K})^{-1}-\mathbf{I}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)},\quad\text{and}$ (6.4)
$\|\widetilde{\bm{\Gamma}}_{K}-\bm{\Gamma}_{K}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{pn}}\Bigg{)}.$ (6.5)

Recall that $\lambda_{0;i}$ stands for the $i$th largest eigenvalue of
$\bm{\Sigma}_{0}$. Define

$\sigma_{a}=\frac{1}{p}\sum_{i=K+1}^{p}\lambda_{0;i}.$ (6.6)

Under Assumption 3, we have that $\sigma_{a}\asymp 1$.

###### Lemma 7

Under Assumptions 1–3, if in addition, $(\log(p))^{2+\gamma}=o(n)$ for some
$\gamma>0$ and $\log(n)=O(\log(p))$, then

$P\Bigg{(}\min_{1\leq t\leq
n}\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}>\frac{\sigma_{a}}{3}\Bigg{)}\to
1,$ (6.7)
$\frac{1}{n}\sum_{t=1}^{n}\Bigg{(}\frac{\xi_{t}^{2}p}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Bigg{)}^{2}=O_{p}\Bigg{(}\frac{\log
p}{p}+\frac{\log p}{n}\Bigg{)},$ (6.8)

and

$\displaystyle\frac{1}{p}\operatorname{tr}\Big{(}\frac{1}{n}\sum_{t=1}^{n}\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T}\Big{)}=\frac{1}{\sigma_{a}}+O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}+\sqrt{\frac{\log p}{p}}\Bigg{)}.$ (6.9)

###### Lemma 8

Under Assumptions 1–3, if in addition, $(\log(p))^{2+\gamma}=o(n)$ for some
$\gamma>0$ and $\log(n)=O(\log(p))$, then

$\|\widehat{\bm{\Sigma}}_{0}-\widetilde{\mathbf{S}}_{0}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p}}+\sqrt{\frac{\log p}{n}}\Bigg{)},$ (6.10)
$\|\widehat{\bm{\Lambda}}_{0K}(\widetilde{\bm{\Lambda}}_{0K})^{-1}-\mathbf{I}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p}}+\sqrt{\frac{\log p}{n}}\Bigg{)},\quad\text{and}$ (6.11)
$\|\widehat{\bm{\Gamma}}_{K}-\widetilde{\bm{\Gamma}}_{K}\|_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{p^{2}}}+\sqrt{\frac{\log p}{pn}}\Bigg{)}.$ (6.12)

###  6.2 Proof of Proposition 1

The desired results follow from Lemma 7. $\hfill\Box$

###  6.3 Proof of Theorem 1

The desired bounds follow from Lemmas 6 and 8. $\hfill\Box$

###  6.4 Proof of Proposition 2

Using the definition of $\widehat{E(\xi_{t}^{2})}$ and the inequality that for
any $x\geq y$,

$\Psi_{H}(x)-\Psi_{H}(y)\geq(x-y)\mathbf{1}_{\big{\\{}x<H\,\text{and}\,y>-H\big{\\}}},$
(6.13)

we have

$|\widehat{E(\xi_{t}^{2})}-E(\xi_{t}^{2})|\leq\frac{1}{\tilde{p}^{H}}\Bigg{|}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-E(\xi_{t}^{2})\Big{)}\Bigg{|},$

where $\tilde{p}^{H}=\sum_{t=1}^{n}\gamma_{t}/n$ for $(\gamma_{t})$ defined as
follows:

$\gamma_{t}=\begin{cases}&\mathbf{1}_{\\{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-\widehat{E(\xi_{t}^{2})}<H\quad\text{and}\quad\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-E(\xi_{t}^{2})>-H\\}}\quad\text{
if }\widehat{E(\xi_{t}^{2})}\leq E(\xi_{t}^{2})\\\
&\mathbf{1}_{\\{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-\widehat{E(\xi_{t}^{2})}>-H\quad\text{and}\quad\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-E(\xi_{t}^{2})<H\\}}\quad\text{
if }\widehat{E(\xi_{t}^{2})}\geq E(\xi_{t}^{2})\end{cases}.$

We will show that

$\Bigg{|}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}/p-E(\xi_{t}^{2})\Big{)}\Bigg{|}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}+\frac{1}{n^{\varepsilon/(2+\varepsilon)}}\Bigg{)},$ (6.14)

and

$\tilde{p}^{H}=1+o_{p}(1).$ (6.15)

The desired bound will follow from (6.14) and (6.15).

We start with (6.14). By (2.1), we have

$\displaystyle\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}=$
$\displaystyle\xi_{t}^{2}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{\mathbf{z}^{T}_{t}\mathbf{z}_{t}}+\frac{\|\widehat{\bm{\mu}}-\bm{\mu}\|_{2}^{2}}{p}+\frac{2}{\sqrt{p}}\xi_{t}\frac{(\bm{\mu}-\widehat{\bm{\mu}})^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}}{\|\mathbf{z}_{t}\|_{2}}$
(6.16) $\displaystyle=:$ $\displaystyle I_{t}+II+III_{t}.$

The definition of $\Psi_{H}(\cdot)$ implies that for any $a,b$,

$|\Psi_{H}(a)-\Psi_{H}(b)|\leq|a-b|.$

Hence, for any $(a_{i})_{1\leq i\leq n}$, $(b_{i})_{1\leq i\leq n}$,

$\Bigg{|}\sum_{t=1}^{n}\Psi_{H}(a_{i})-\sum_{t=1}^{n}\Psi_{H}(b_{i})\Bigg{|}\leq\sum_{t=1}^{n}|a_{i}-b_{i}|.$
(6.17)

By (6.16) and (6.17), we get

$\displaystyle\Bigg{|}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Bigg{)}\Bigg{|}$
(6.18) $\displaystyle\leq$
$\displaystyle\Bigg{|}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}I_{t}-E(\xi_{t}^{2})\Big{)}\Bigg{|}+II+\frac{1}{n}\sum_{t=1}^{n}|III_{t}|.$

By Lemma 5, and the assumptions that $E(\xi_{t}^{2+\varepsilon})<\infty$ and
$(\xi_{t})$ is independent with $(\mathbf{z}_{t})$, we have

$E(I_{t})=E(\xi_{t}^{2})\text{, and }E(I_{t}^{1+\varepsilon/2})<\infty.$
(6.19)

By (6.19) and equation (S31) in avella2018robust, for $H\asymp
n^{1/(1+\varepsilon/2)}$,

$\Bigg{|}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}I_{t}-E(\xi_{t}^{2})\Big{)}\Bigg{|}=O_{p}\Bigg{(}\frac{1}{n^{\varepsilon/(2+\varepsilon)}}\Bigg{)}.$
(6.20)

By Lemma 4,

$II=O_{p}\Bigg{(}\frac{\log p}{n}\Bigg{)}.$ (6.21)

By the Cauchy-Schwarz inequality,

$\displaystyle\frac{1}{n}\sum_{t=1}^{n}|III_{t}|\leq$
$\displaystyle\frac{1}{\sqrt{p}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\xi^{2}_{t}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\frac{\big{(}(\bm{\mu}-\widehat{\bm{\mu}})^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\big{)}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}}$
(6.22) $\displaystyle\leq$
$\displaystyle\frac{\|\bm{\mu}-\widehat{\bm{\mu}}\|_{2}}{\sqrt{p}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\xi^{2}_{t}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{\|\mathbf{z}_{t}\|^{2}_{2}}}.$

By the assumption that $E(\xi_{t}^{2+\varepsilon})<\infty$, (6.1) and Markov’s
inequality,

$\frac{1}{n}\sum_{t=1}^{n}\xi^{2}_{t}=O_{p}(1),\text{ and
}\,\frac{1}{n}\sum_{t=1}^{n}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{\|\mathbf{z}_{t}\|^{2}_{2}}=O_{p}(1).$
(6.23)

Combining (6.21), (6.22) and (6.23) yields that

$\frac{1}{n}\sum_{t=1}^{n}|III_{t}|=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)}.$ (6.24)

The desired bound (6.14) follows from (6.18), (6.20), (6.21) and (6.24).

Finally, about (6.15), if $\widehat{E(\xi_{t}^{2})}\leq E(\xi_{t}^{2})$, we
have that for all $n$ large enough,

$\displaystyle\gamma_{t}\geq$
$\displaystyle\mathbf{1}_{\big{\\{}\big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\big{|}\leq\frac{H}{2}\quad\text{and}\quad\widehat{E(\xi_{t}^{2})}\geq
0\big{\\}}}$ $\displaystyle\geq$ $\displaystyle
1-\mathbf{1}_{\big{\\{}\big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\big{|}>\frac{H}{2}\big{\\}}}-\mathbf{1}_{\big{\\{}\widehat{E(\xi_{t}^{2})}<0\big{\\}}},$

where in the first inequality we use $H\asymp
n^{1/(1+\varepsilon/2)}\to\infty$ so that $E(\xi_{t}^{2})<H/2$ for all large
$n$. Similarly, if $\widehat{E(\xi_{t}^{2})}\geq E(\xi_{t}^{2})$,

$\gamma_{t}\geq
1-\mathbf{1}_{\big{\\{}\big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\big{|}>\frac{H}{2}\big{\\}}}-\mathbf{1}_{\big{\\{}\widehat{E(\xi_{t}^{2})}>\frac{H}{2}\big{\\}}}.$

It follows that

$E(\gamma_{t})\geq
1-P\Big{(}\Big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Big{|}>\frac{H}{2}\Big{)}-P\Big{(}\widehat{E(\xi_{t}^{2})}>\frac{H}{2}\Big{)}-P\Big{(}\widehat{E(\xi_{t}^{2})}<0\Big{)}.$
(6.25)

We have

$\displaystyle
P\Bigg{(}\Big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Big{|}>\frac{H}{2}\Bigg{)}$
$\displaystyle\leq$ $\displaystyle
P\Bigg{(}|I_{t}-E(\xi_{t}^{2})|>\frac{H}{6}\Bigg{)}+P\Bigg{(}|II|>\frac{H}{6}\Bigg{)}+P\Bigg{(}|III_{t}|>\frac{H}{6}\Bigg{)},$

where terms $I_{t}$, $II$ and $III_{t}$ are defined in equation (6.16). By
(6.19), $H\asymp n^{1/(1+\varepsilon/2)}$, and Markov’ inequality, we have

$P\Bigg{(}|I_{t}-E(\xi_{t}^{2})|>\frac{H}{6}\Bigg{)}=O\Bigg{(}\frac{1}{H^{1+\varepsilon/2}}\Bigg{)}=O\Bigg{(}\frac{1}{n}\Bigg{)}\to
0.$

By Lemma 1, the assumptions that $\log p=O(n)$ and $H\asymp
n^{1/(1+\varepsilon/2)}\to\infty$, we have

$P\Bigg{(}II>\frac{H}{6}\Bigg{)}=o(1).$

By the Cauchy-Schwarz inequality, $|III_{t}|\leq\sqrt{II}\sqrt{I_{t}}$. We
then get that

$\displaystyle P\Bigg{(}|III_{t}|>\frac{H}{6}\Bigg{)}\leq$ $\displaystyle
P\Bigg{(}I_{t}>\frac{H}{6}\Bigg{)}+P\Bigg{(}II>\frac{H}{6}\Bigg{)}$
$\displaystyle=$ $\displaystyle o(1).$

Combining the results above yields that

$P\Bigg{(}\Big{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Big{|}>\frac{H}{2}\Bigg{)}=o(1).$
(6.26)

Using (6.13), we have

$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-\frac{H}{2}\Bigg{)}$
(6.27) $\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Bigg{)}$
$\displaystyle-\Big{(}\frac{H}{2}-E(\xi_{t}^{2})\Big{)}\frac{1}{n}\sum_{t=1}^{n}\mathbf{1}_{\Big{\\{}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})<H\,\text{and}\,\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}>-\frac{H}{2}\Big{\\}}}.$

By (6.26), we have

$P\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})<H\Bigg{)}\geq
1-P\Bigg{(}\Bigg{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Bigg{|}\geq\frac{H}{2}\Bigg{)}\to
1,$

and

$P\Bigg{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})>-\frac{H}{2}\Bigg{)}\geq
1-P\Bigg{(}\Bigg{|}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})\Bigg{|}\geq\frac{H}{2}\Bigg{)}\to
1.$

Therefore,

$P\Bigg{(}\Big{(}\frac{H}{2}-E(\xi_{t}^{2})\Big{)}\frac{1}{n}\sum_{t=1}^{n}\mathbf{1}_{\Big{\\{}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-E(\xi_{t}^{2})<H\,\text{and}\,\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}>-\frac{H}{2}\Big{\\}}}>\frac{H}{4}\Bigg{)}\to
1.$ (6.28)

Combining (6.14), (6.27) and (6.28) yields that

$P\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}-\frac{H}{2}\Big{)}<-\frac{H}{5}\Bigg{)}\to
1.$

Similarly, one can show that

$P\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Psi_{H}\Big{(}\frac{\|\mathbf{y}_{t}-\widehat{\bm{\mu}}\|_{2}^{2}}{p}\Big{)}>\frac{H}{5}\Bigg{)}\to
1.$

Using the definition of $\widehat{E(\xi_{t}^{2})}$ and the continuity of the
function $\Psi_{H}(\cdot)$, we then get that

$P\Bigg{(}0<\widehat{E(\xi_{t}^{2})}<\frac{H}{2}\Bigg{)}\to 1.$ (6.29)

Combining (6.25), (6.26) and (6.29) yields that

$E(\gamma_{t})=1-o(1).$

By Markov’s inequality and that $\tilde{p}_{H}\leq 1$, we get the desired
result (6.15). $\hfill\Box$

###  6.5 Proof of Theorems 2 and 3

The bounds in (2.15) follow from Theorem 1 and Theorem 2.1 of fan2018large.
The bounds in (2.16) follow from (2.15) and Proposition 2.$\hfill\Box$

## Supplementary Materials

Supplement to “Sub-Gaussian High-dimensional covariance matrix estimation
under elliptical factor model with $2+\varepsilon$th moment.”

This supplement contains the proofs of Propositions 1-2, and Theorems 1–3.

## References

  * Ahn and Horenstein, (2013) Ahn, S. C. and Horenstein, A. R. (2013).  Eigenvalue ratio test for the number of factors.  Econometrica, 81(3):1203–1227.
  * Avella-Medina et al., (2018) Avella-Medina, M., Battey, H. S., Fan, J., and Li, Q. (2018).  Robust estimation of high-dimensional covariance and precision matrices.  Biometrika, 105(2):271–284.
  * Bai and Ng, (2002) Bai, J. and Ng, S. (2002).  Determining the number of factors in approximate factor models.  Econometrica, 70(1):191–221.
  * Bollerslev et al., (2013) Bollerslev, T., Todorov, V., and Li, S. Z. (2013).  Jump tails, extreme dependencies, and the distribution of stock returns.  Journal of Econometrics, 172(2):307–324.
  * Cai and Liu, (2011) Cai, T. and Liu, W. (2011).  Adaptive thresholding for sparse covariance matrix estimation.  Journal of the American Statistical Association, 106(494):672–684.
  * Cai et al., (2016) Cai, T. T., Liu, W., and Zhou, H. H. (2016).  Estimating sparse precision matrix: Optimal rates of convergence and adaptive estimation.  The Annals of Statistics, 44(2):455–488.
  * Chatterjee, (2013) Chatterjee, S. (2013).  Assumptionless consistency of the lasso.  arXiv preprint arXiv:1303.5817.
  * Ding et al., (2023) Ding, Y., Li, Y., Liu, G., and Zheng, X. (2023).  Stock co-jump networks.  Journal of Econometrics.
  * Ding et al., (2021) Ding, Y., Li, Y., and Zheng, X. (2021).  High dimensional minimum variance portfolio estimation under statistical factor models.  Journal of Econometrics, 222(1):502–515.
  * (10) Ding, Y. and Zheng, X. (2024a).  Sub-gaussian high-dimensional covariance matrix estimation under elliptical factor models with 2+$\varepsilon$th moment.
  * (11) Ding, Y. and Zheng, X. (2024b).  Supplement to “sub-gaussian high-dimensional covariance matrix estimation under elliptical factor models with 2+$\varepsilon$th moment”.
  * Fan et al., (2008) Fan, J., Fan, Y., and Lv, J. (2008).  High dimensional covariance matrix estimation using a factor model.  Journal of Econometrics, 147(1):186–197.
  * Fan et al., (2012) Fan, J., Li, Y., and Yu, K. (2012).  Vast volatility matrix estimation using high-frequency data for portfolio selection.  Journal of the American Statistical Association, 107(497):412–428.
  * Fan et al., (2011) Fan, J., Liao, Y., and Mincheva, M. (2011).  High-dimensional covariance matrix estimation in approximate factor models.  The Annals of Statistics, 39(6):3320–3356.
  * Fan et al., (2013) Fan, J., Liao, Y., and Mincheva, M. (2013).  Large covariance estimation by thresholding principal orthogonal complements.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(4):603–680.
  * Fan et al., (2018) Fan, J., Liu, H., and Wang, W. (2018).  Large covariance estimation through elliptical factor models.  The Annals of Statistics, 46(4):1383.
  * Fan et al., (2021) Fan, J., Wang, W., and Zhu, Z. (2021).  A shrinkage principle for heavy-tailed data: High-dimensional robust low-rank matrix recovery.  The Annals of Statistics, 49(3):1239.
  * Han and Liu, (2014) Han, F. and Liu, H. (2014).  Scale-invariant sparse pca on high-dimensional meta-elliptical data.  Journal of the American Statistical Association, 109(505):275–287.
  * Han and Liu, (2018) Han, F. and Liu, H. (2018).  Eca: High-dimensional elliptical component analysis in non-gaussian distributions.  Journal of the American Statistical Association, 113(521):252–268.
  * Hubert et al., (2018) Hubert, M., Debruyne, M., and Rousseeuw, P. J. (2018).  Minimum covariance determinant and extensions.  Wiley Interdisciplinary Reviews: Computational Statistics, 10(3):e1421.
  * Ledoit and Wolf, (2011) Ledoit, O. and Wolf, M. (2011).  Robust performances hypothesis testing with the variance.  Wilmott, 2011(55):86–89.
  * Ledoit and Wolf, (2017) Ledoit, O. and Wolf, M. (2017).  Nonlinear shrinkage of the covariance matrix for portfolio selection: Markowitz meets goldilocks.  The Review of Financial Studies, 30(12):4349–4388.
  * Liu et al., (2012) Liu, H., Han, F., Yuan, M., Lafferty, J., and Wasserman, L. (2012).  High-dimensional semiparametric gaussian copula graphical models.  The Annals of Statistics, 40(4):2293–2326.
  * Maronna, (1976) Maronna, R. A. (1976).  Robust m-estimators of multivariate location and scatter.  The Annals of Statistics, pages 51–67.
  * Minsker, (2018) Minsker, S. (2018).  Sub-gaussian estimators of the mean of a random matrix with heavy-tailed entries.  The Annals of Statistics, 46(6A):2871–2903.
  * Minsker and Wei, (2020) Minsker, S. and Wei, X. (2020).  Robust modifications of u-statistics and applications to covariance estimation problems.  Bernoulli, 26(1):694–727.
  * Onatski, (2010) Onatski, A. (2010).  Determining the number of factors from empirical distribution of eigenvalues.  The Review of Economics and Statistics, 92(4):1004–1016.
  * Rousseeuw, (1984) Rousseeuw, P. J. (1984).  Least median of squares regression.  Journal of the American Statistical Association, 79(388):871–880.
  * Rousseeuw, (1985) Rousseeuw, P. J. (1985).  Multivariate estimation with high breakdown point.  Mathematical Statistics and Applications, 8(283-297):37.
  * Sun et al., (2020) Sun, Q., Zhou, W.-X., and Fan, J. (2020).  Adaptive huber regression.  Journal of the American Statistical Association, 115(529):254–265.
  * Tyler, (1987) Tyler, D. E. (1987).  A distribution-free m-estimator of multivariate scatter.  The Annals of Statistics, pages 234–251.
  * Van Aelst and Rousseeuw, (2009) Van Aelst, S. and Rousseeuw, P. (2009).  Minimum volume ellipsoid.  Wiley Interdisciplinary Reviews: Computational Statistics, 1(1):71–82.
  * Xue and Zou, (2012) Xue, L. and Zou, H. (2012).  Regularized rank-based estimation of high-dimensional nonparanormal graphical models.  The Annals of Statistics, 40(5):2541–2571.
  * Yu et al., (2019) Yu, L., He, Y., and Zhang, X. (2019).  Robust factor number specification for large-dimensional elliptical factor model.  Journal of Multivariate Analysis, 174:104543.
  * Zhang, (2016) Zhang, T. (2016).  Robust subspace recovery by tyler’s m-estimator.  Information and Inference: A Journal of the IMA, 5(1):1–21.
  * Zhang et al., (2016) Zhang, T., Cheng, X., and Singer, A. (2016).  Marčenko–pastur law for tyler’s m-estimator.  Journal of Multivariate Analysis, 149:114–123.
  * Zheng and Li, (2011) Zheng, X. and Li, Y. (2011).  On the estimation of integrated covariance matrices of high dimensional diffusion processes.  The Annals of Statistics, 39(6):3121–3151.

Supplement to “Sub-Gaussian High-Dimensional Covariance Matrix Estimation
under Elliptical Factor Model with $2+\varepsilon$th Moment”

Yi Ding and Xinghua Zheng

##  Appendix A Proofs

The theorem/proposition/lemma/equation numbers below refer to the main article
(DZ24). In the following proofs, $c$ and $C$ denote generic constants which do
not depend on $p$ and $n$ and can vary from place to place.

## Proof of Lemma 4

Under Assumptions 1 and 2, we have $\max_{1\leq i\leq
p}(\bm{\Sigma})_{ii}=O(1)$. The desired bound follows from Lemma 1 with
$\delta=O(1/\log p)$ and $\varepsilon=1$, Bonferroni’s inequality, and the
assumption that $\log(p)=O(n)$. $\hfill\Box$

## Proof of Lemma 5

Write the eigen decomposition of $\bm{\Sigma}$ as
$\bm{\Sigma}=\mathbf{\Gamma}\bm{\Lambda}\mathbf{\Gamma}^{T}$ with
$\bm{\Lambda}=\operatorname{diag}(\lambda_{1},...,\lambda_{p})$. Denote
$\mathbf{z}_{\Gamma;t}=\mathbf{\Gamma}^{T}\mathbf{z}_{t}=:(z_{\Gamma;1t},...,z_{\Gamma;pt})^{T}$.
Because $\mathbf{z}_{t}\sim N(0,\mathbf{I})$, and
$\mathbf{\Gamma}\mathbf{\Gamma}^{T}=\mathbf{\Gamma}^{T}\mathbf{\Gamma}=\mathbf{I}$,
we have that $\mathbf{z}_{\Gamma}\sim N(0,\mathbf{I})$, and

$\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}\mathbf{z}_{t}}{\mathbf{z}^{T}_{t}\mathbf{z}_{t}}=\sum_{i=1}^{N}\lambda_{i}\frac{z_{\Gamma;it}^{2}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{2}}.$
(A.1)

By symmetry,

$E\Bigg{(}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}\mathbf{z}_{t}}{\mathbf{z}^{T}_{t}\mathbf{z}_{t}}\Bigg{)}=\frac{1}{p}\operatorname{tr}(\bm{\Sigma}).$

Equation (6.1) follows.

Next, by (A.1), we have

$\displaystyle
E\Bigg{(}\frac{(\mathbf{z}_{t}^{T}\bm{\Sigma}\mathbf{z}_{t})^{2}}{(\mathbf{z}^{T}_{t}\mathbf{z}_{t})^{2}}\Bigg{)}=$
$\displaystyle\sum_{i,j=1}^{p}\lambda_{i}\lambda_{j}E\Bigg{(}\frac{z_{\Gamma;it}^{2}z_{\Gamma;jt}^{2}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}$
$\displaystyle\leq$
$\displaystyle\big{(}\operatorname{tr}(\bm{\Sigma}^{2})\big{)}E\Bigg{(}\frac{z_{\Gamma;1t}^{4}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}+\big{(}\operatorname{tr}(\bm{\Sigma})\big{)}^{2}E\Bigg{(}\frac{z_{\Gamma;1t}^{2}z_{\Gamma;2t}^{2}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}$
$\displaystyle\leq$
$\displaystyle\frac{\operatorname{tr}(\bm{\Sigma}^{2})}{p}+\frac{\big{(}\operatorname{tr}(\bm{\Sigma})\big{)}^{2}}{(p-1)p},$

where the last inequality holds because

$\sum_{1\leq i,j\leq
p}z_{\Gamma;it}^{2}z_{\Gamma;jt}^{2}=\|\mathbf{z}_{\Gamma;t}\|_{2}^{4},$

which implies that

$p(p-1)E\Bigg{(}\frac{z_{\Gamma;1t}^{2}z_{\Gamma;2t}^{2}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}=\sum_{1\leq
i\neq j\leq
p}E\Bigg{(}\frac{z_{\Gamma;it}^{2}z_{\Gamma;jt}^{2}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}\leq
1,$

and

$pE\Bigg{(}\frac{z_{\Gamma;1t}^{4}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}=\sum_{i=1}^{p}E\Bigg{(}\frac{z_{\Gamma;it}^{4}}{\|\mathbf{z}_{\Gamma;t}\|_{2}^{4}}\Bigg{)}\leq
1.$

Inequality (6.2) follows. $\hfill\Box$

## Proof of Lemma 6

The desired bounds are immediate results of Lemmas A.2, A.3 and A.5 of
DLZ19.$\hfill\Box$

For any vector $\mathbf{x}$, let $(\mathbf{x})_{i}$ be its $i$th entry of
$\mathbf{x}$, and for any matrix $\mathbf{A}$, let $(\mathbf{A})_{ij}$ be its
$(i,j)$th entry.

## Proof of Lemma 7

We start with (6.7). Using (2.1), we have

$\displaystyle\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}$
(A.2) $\displaystyle=$
$\displaystyle\frac{\|\widehat{\mathbf{P}}_{I}({\bm{\mu}}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}+\frac{\|\widehat{\mathbf{P}}_{I}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}+2\frac{({\bm{\mu}}-\widehat{\bm{\mu}})^{T}\widehat{\mathbf{P}}^{T}_{I}\widehat{\mathbf{P}}_{I}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}}{\sqrt{p}\xi_{t}\|\mathbf{z}_{t}\|_{2}}$
$\displaystyle=:$ $\displaystyle I_{t}+II_{t}+III_{t}.$

By Lemma 4 and the assumption that $\xi_{t}>c$,

$\max_{1\leq t\leq n}|I_{t}|=O_{p}\Bigg{(}\frac{\log p}{n}\Bigg{)}.$

Define $\mathbf{\Gamma}_{I}=(\phi_{K+1},...,\phi_{p})$. About term $II_{t}$,
we have

$\displaystyle II_{t}=$
$\displaystyle\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}+\frac{\|(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}+2\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{\Gamma}_{I}(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}}{\|\mathbf{z}_{t}\|^{2}_{2}}.$
(A.3)

We have
$E(\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}/p)=\sigma_{a}$,
and by Lemma 2,

$E\Bigg{(}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2s}\Bigg{)}\leq\frac{c_{s}}{p^{s}},\text{
for all }s\geq 1.$ (A.4)

Under the assumption that $\log n=O(\log p)$, there exists $\delta>0$ such
that $n\leq p^{\delta}$ for all $n$ large enough. We then set $s$ sufficient
large such that $s>\delta+1$ and get

$\displaystyle P\Bigg{(}\max_{1\leq t\leq
n}\Big{|}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{|}>\frac{\sigma_{a}}{5}\Bigg{)}$
(A.5) $\displaystyle\leq$ $\displaystyle
cnE\Bigg{(}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2s}\Bigg{)}$
$\displaystyle=$ $\displaystyle O\Bigg{(}\frac{n}{p^{s}}\Bigg{)}\to 0.$

By Lemma 2 again,

$\displaystyle P\Bigg{(}\max_{1\leq t\leq
n}\Big{|}\frac{\|\mathbf{z}_{t}\|_{2}^{2}}{p}-1\Big{|}>\frac{1}{5}\Bigg{)}$
(A.6) $\displaystyle\leq$ $\displaystyle
cnE\Bigg{(}\Big{(}\frac{\|\mathbf{z}_{t}\|_{2}^{2}}{p}-1\Big{)}^{2s}\Bigg{)}$
$\displaystyle=$ $\displaystyle O\Bigg{(}\frac{n}{p^{s}}\Bigg{)}\to 0.$

By Assumption 2 and Lemma 3, we have $E\Big{(}\max_{1\leq t\leq n,1\leq i\leq
p}|(\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t})_{i}|\Big{)}=O(\sqrt{\log(pn)})$.
Hence by the assumptions that $\log n=O(\log p)$ and
$(\log(p))^{2+\gamma}=o(n)$ for some $\gamma>0$, we have

$\displaystyle P\Big{(}\max_{1\leq t\leq n,1\leq i\leq
p}(\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t})_{i}\geq(\log(p))^{(1+\gamma)/2}\Big{)}$
$\displaystyle\leq$ $\displaystyle\frac{E\Big{(}\max_{1\leq t\leq n,1\leq
i\leq
p}|(\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t})_{i}|\Big{)}}{(\log(p))^{(1+\gamma)/2}}$
$\displaystyle=$ $\displaystyle O\Bigg{(}\frac{(\log
p+\log{n})^{1/2}}{(\log(p))^{(1+\gamma)/2}}\Bigg{)}=O\Bigg{(}\frac{1}{(\log(p))^{\gamma/2}}\Bigg{)}\to
0.$

In addition, by Theorem 4.1 of fan2018large111Note that the Spatial Kendall’s
tau is independent of the scalar process $(\xi_{t})$, see P.1402 in
fan2018large., we get that
$\|\widehat{\bm{\Gamma}_{K}}_{ED}-\bm{\Gamma}_{K}\|_{2}=O_{p}(\sqrt{(\log
p)/n})$. It follows that

$\|\mathbf{\Gamma}_{I}^{T}-\widehat{\mathbf{P}}_{I}\|_{2}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)}.$ (A.7)

Therefore,

$\displaystyle\max_{1\leq t\leq
n}\frac{\|(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}$
(A.8) $\displaystyle\leq$
$\displaystyle\|(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\|_{2}^{2}\max_{1\leq
t\leq n,1\leq i\leq p}|(\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t})_{i}^{2}|$
$\displaystyle=$ $\displaystyle O_{p}\Bigg{(}\frac{(\log
p)^{2+\gamma}}{n}\Bigg{)}=o_{p}(1),$

where the last equality holds by the assumption that
$(\log(p))^{2+\gamma}=o(n)$. By (A.5) and (A.8), we get

$\displaystyle\max_{1\leq t\leq
n}\frac{1}{p}\Big{|}\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{\Gamma}_{I}(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\Big{|}$
(A.9) $\displaystyle\leq$ $\displaystyle\frac{1}{p}\Bigg{(}\max_{1\leq t\leq
n}\|(\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T})\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}\Bigg{)}\Bigg{(}\max_{1\leq
t\leq
n}\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}\Bigg{)}$
$\displaystyle=$ $\displaystyle O_{p}\Bigg{(}\frac{(\log
p)^{1+\gamma/2}}{\sqrt{n}}\Bigg{)}=o_{p}(1).$

Combining (A.3), (A.5), (A.6), (A.8) and (A.9) yields

$P\Bigg{(}\max_{1\leq t\leq
n}|II_{t}-\sigma_{a}|>\frac{\sigma_{a}}{3}\Bigg{)}=o(1).$

About term $III_{t}$, we have

$\max_{1\leq t\leq n}|III_{t}|\leq\sqrt{\Big{(}\max_{1\leq t\leq
n}I_{t}\Big{)}\Big{(}\max_{1\leq t\leq n}II_{t}\Big{)}}=o_{p}(1).$

The desired result (6.7) follows.

Next, we show (6.8). We have

$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\Bigg{(}\frac{p\xi_{t}^{2}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Bigg{)}^{2}$
(A.10) $\displaystyle\leq$ $\displaystyle\frac{1}{\sigma_{a}^{2}\min_{1\leq
t\leq
n}\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{4}}{p^{2}\xi_{t}^{4}}}\sum_{t=1}^{n}\frac{1}{n}\Big{(}\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}-\sigma_{a}\Big{)}^{2}$
$\displaystyle\leq$ $\displaystyle\frac{2}{\sigma_{a}^{2}\min_{1\leq t\leq
n}\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{4}}{p^{2}\xi_{t}^{4}}}\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}-\sigma_{a}\Big{)}^{2}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}-\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}\Big{)}^{2}\Bigg{)}$
$\displaystyle=:$ $\displaystyle\frac{2}{\sigma_{a}^{2}\min_{1\leq t\leq
n}\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{4}}{p^{2}\xi_{t}^{4}}}(I+II).$

About term $I$, we have

$\displaystyle I\leq$
$\displaystyle\frac{2}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2}+\frac{2}{n}\sum_{t=1}^{n}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}\Big{(}\frac{\|\mathbf{z}_{t}\|_{2}^{2}}{p}-1\Big{)}^{2}$
$\displaystyle\leq$
$\displaystyle\frac{2}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2}+\frac{2C}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{z}_{t}\|_{2}^{2}}{p}-1\Big{)}^{2},$

where the last inequality holds because

$\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}\leq\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}\mathbf{\Gamma}_{I}\|_{2}=O(1).$

By Jensen’s inequality and (A.4), for all $k\geq 1$,

$\displaystyle
E\Bigg{(}\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2}\Bigg{)}^{k}\Bigg{)}$
$\displaystyle\leq$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}E\Bigg{(}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2k}\Bigg{)}$
$\displaystyle=$ $\displaystyle O\Bigg{(}\frac{1}{p^{k}}\Bigg{)}.$

Hence,

$\displaystyle
P\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{p}-\sigma_{a}\Big{)}^{2}\geq\frac{\log
p}{p}\Bigg{)}\leq\frac{c}{\log p}\to 0.$

Similarly,

$\displaystyle
P\Bigg{(}\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{\|\mathbf{z}_{t}\|_{2}^{2}}{p}-1\Big{)}^{2}\geq\frac{\log
p}{p}\Bigg{)}\leq\frac{c}{\log p}\to 0.$ (A.11)

It follows that

$I=O_{p}\Bigg{(}\frac{\log p}{p}\Bigg{)}.$ (A.12)

About term $II$, by (A.2) and (A.3), and the inequality that
$(\sum_{k=1}^{K}x_{k})^{2}\leq 2^{K-1}\sum_{k=1}^{K}x_{k}^{2}$, we have

$\displaystyle\Bigg{(}\frac{\|\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}^{1/2}\mathbf{z}_{t}\|_{2}^{2}}{\|\mathbf{z}_{t}\|_{2}^{2}}-\frac{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}{p\xi_{t}^{2}}\Bigg{)}^{2}$
(A.13) $\displaystyle\leq$ $\displaystyle
8\frac{\|\widehat{\mathbf{P}}_{I}({\bm{\mu}}-\widehat{\bm{\mu}})\|_{2}^{4}}{p^{2}\xi_{t}^{4}}$
$\displaystyle+8\frac{\|\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T}\|_{2}^{4}(\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t})^{2}}{\|\mathbf{z}_{t}\|_{2}^{4}}$
$\displaystyle+32\frac{\|{\bm{\mu}}-\widehat{\bm{\mu}}\|_{2}^{4}\cdot\|\widehat{\mathbf{P}}_{I}-\mathbf{\Gamma}_{I}^{T}\|_{2}^{4}(\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t})^{2}}{p^{2}\xi^{4}_{t}\|\mathbf{z}_{t}\|_{2}^{4}}$
$\displaystyle+32\frac{\|{\bm{\mu}}-\widehat{\bm{\mu}}\|_{2}^{4}(\mathbf{z}_{t}^{T}\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}\mathbf{\Gamma}_{I}\mathbf{z}_{t})^{2}}{p^{2}\xi^{4}_{t}\|\mathbf{z}_{t}\|_{2}^{4}}.$

By Lemma 2, we have that

$E\Bigg{(}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{p}\Bigg{)}=1,\quad\text{
and }\quad
E\Bigg{(}\Big{(}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{p}\Big{)}^{2}\Bigg{)}=O(1).$

Hence,

$\frac{1}{n}\sum_{t=1}^{n}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{p}=O_{p}(1),\quad\text{
and
}\quad\frac{1}{n}\sum_{t=1}^{n}\Bigg{(}\frac{\mathbf{z}_{t}^{T}\bm{\Sigma}_{0}\mathbf{z}_{t}}{p}\Bigg{)}^{2}=O_{p}(1).$
(A.14)

Because
$E(\mathbf{z}_{t}^{T}\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}\mathbf{\Gamma}_{I}\mathbf{z}_{t}/p)=\sigma_{a}=O(1)$,
we also have

$\frac{1}{n}\sum_{t=1}^{n}\frac{\mathbf{z}_{t}^{T}\mathbf{\Gamma}_{I}^{T}\bm{\Sigma}_{0}\mathbf{\Gamma}_{I}\mathbf{z}_{t}}{p}=O_{p}(1).$
(A.15)

Combining Lemma 4, (A.6), (A.7), (A.13), (A.14) and (A.15) yields that

$II=O_{p}\Bigg{(}\frac{\log p}{n}\Bigg{)}.$ (A.16)

The desired bound (6.8) follows from (6.7), (A.10), (A.12) and (A.16).

Finally, we show (6.9). Recall that

$\widehat{\mathbf{X}}_{t}=\frac{\sqrt{p}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}}(\mathbf{y}_{t}-\widehat{\bm{\mu}})=\frac{\sqrt{p}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}}\Bigg{(}\bm{\mu}-\widehat{\bm{\mu}}+\frac{\sqrt{p}{\xi_{t}}}{\|\mathbf{z}_{t}\|_{2}}\widetilde{\mathbf{X}}_{t}\Bigg{)}.$

We have

$\displaystyle\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T}=$
$\displaystyle\frac{p}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}(\bm{\mu}-\widehat{\bm{\mu}})(\bm{\mu}-\widehat{\bm{\mu}})^{T}$
$\displaystyle+\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}(\bm{\mu}-\widehat{\bm{\mu}})\widetilde{\mathbf{X}}_{t}^{T}$
$\displaystyle+\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}\widetilde{\mathbf{X}}_{t}(\bm{\mu}-\widehat{\bm{\mu}})^{T}$
$\displaystyle+\Big{(}\frac{p^{2}\xi^{2}_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Big{)}\widetilde{\mathbf{X}}_{t}\widetilde{\mathbf{X}}_{t}^{T}$
$\displaystyle+\frac{1}{\sigma_{a}}\widetilde{\mathbf{X}}_{t}\widetilde{\mathbf{X}}_{t}^{T}.$

Therefore, for any $1\leq i,j\leq p$,

$\displaystyle\Big{(}\frac{1}{n}\sum_{t=1}^{n}\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T}\Big{)}_{ij}=\frac{1}{\sigma_{a}}(\widetilde{\mathbf{S}}_{0})_{ij}$
$\displaystyle+\frac{1}{n}\sum_{t=1}^{n}\frac{p}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}(\bm{\mu}-\widehat{\bm{\mu}})_{i}(\bm{\mu}-\widehat{\bm{\mu}})_{j}$
$\displaystyle+\frac{1}{n}\sum_{t=1}^{n}\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}(\bm{\mu}-\widehat{\bm{\mu}})_{i}(\widetilde{\mathbf{X}}_{t})_{j}$
$\displaystyle+\frac{1}{n}\sum_{t=1}^{n}\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}(\widetilde{\mathbf{X}}_{t})_{i}(\bm{\mu}-\widehat{\bm{\mu}})_{j}$
$\displaystyle+\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{p^{2}\xi^{2}_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Big{)}(\widetilde{\mathbf{X}}_{t})_{i}(\widetilde{\mathbf{X}}_{t})_{j}$
$\displaystyle=:\frac{1}{\sigma_{a}}(\widetilde{\mathbf{S}}_{0})_{ij}$
$\displaystyle+I_{ij}+II_{ij}+III_{ij}+IV_{ij}.$

About term $I_{ij}$, by Lemma 4, (6.7) and the assumption that $\xi_{t}>c$, we
have

$\max_{1\leq i,j\leq p}|I_{ij}|=O_{p}\Bigg{(}\frac{\log p}{n}\Bigg{)}.$

About terms $II_{ij}$ and $III_{ij}$, we have

$\displaystyle\max_{1\leq i,j\leq p}\big{(}|II_{ij}|,|III_{ij}|\big{)}$
$\displaystyle\leq$ $\displaystyle\Bigg{(}\max_{1\leq t\leq
n}\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}\Bigg{)}\cdot\Bigg{(}\max_{1\leq
i,j\leq
p}\frac{1}{n}\sum_{t=1}^{n}|(\widetilde{\mathbf{X}}_{t})_{i}|\cdot|(\bm{\mu}-\widehat{\bm{\mu}})_{j}|\Bigg{)}$
$\displaystyle\leq$ $\displaystyle\Bigg{(}\max_{1\leq t\leq
n}\frac{p^{3/2}\xi_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}\|\mathbf{z}_{t}\|_{2}}\Bigg{)}\cdot\sqrt{\max_{1\leq
i\leq
p}\frac{1}{n}\sum_{t=1}^{n}\big{(}(\widetilde{\mathbf{X}}_{t})_{i}\big{)}^{2}}\cdot\|\bm{\mu}-\widehat{\bm{\mu}}\|_{\max}.$

It follows from (6.3), (6.7) and (A.6) that

$\max_{1\leq i,j\leq
p}\big{(}|II_{ij}|,|III_{ij}|\big{)}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)}.$

About term $IV_{ij}$, we have

$\displaystyle|IV_{ij}|\leq$
$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\frac{p}{\|\mathbf{z}_{t}\|_{2}^{2}}\Big{|}\frac{p\xi^{2}_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Big{|}\cdot|\widetilde{x}_{it}\widetilde{x}_{jt}|$
$\displaystyle+\frac{1}{n}\sum_{t=1}^{n}\frac{1}{\sigma_{a}}\Bigg{|}\frac{p}{\|\mathbf{z}_{t}\|_{2}^{2}}-1\Bigg{|}|\widetilde{x}_{it}\widetilde{x}_{jt}|,$

where $\widetilde{x}_{it}=(\widetilde{\mathbf{X}}_{t})_{i}$. By the Cauchy-
Schwarz inequality, we have

$\displaystyle\max_{1\leq i,j\leq p}|IV_{ij}|\leq$
$\displaystyle\Bigg{(}\max_{1\leq t\leq
n}\frac{p}{\|\mathbf{z}_{t}\|_{2}^{2}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\Big{|}\frac{p\xi^{2}_{t}}{\|\widehat{\mathbf{P}}_{I}(\mathbf{y}_{t}-\widehat{\bm{\mu}})\|_{2}^{2}}-\frac{1}{\sigma_{a}}\Big{|}^{2}}\Bigg{)}\cdot\sqrt{\max_{1\leq
i\leq p}\frac{1}{n}\sum_{t=1}^{n}\widetilde{x}_{it}^{4}}$
$\displaystyle+\frac{1}{\sigma_{a}}\sqrt{\frac{1}{n}\sum_{t=1}^{n}\Big{(}\frac{p}{\|\mathbf{z}_{t}\|_{2}^{2}}-1\Big{)}^{2}}\sqrt{\max_{1\leq
i\leq p}\frac{1}{n}\sum_{t=1}^{n}\widetilde{x}^{4}_{it}}.$

By Bernstein’s inequality, we have $\max_{1\leq i\leq
p}(\sum_{t=1}^{n}\widetilde{x}_{it}^{4}/n)=O_{p}(1)$. By (6.8), (A.6) and
(A.11), we then get

$\max_{1\leq i\leq p}|IV_{ij}|=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}+\sqrt{\frac{\log p}{p}}\Bigg{)}.$

Combining the results above yields

$\Big{\|}\Big{(}\frac{1}{n}\sum_{t=1}^{n}\widehat{\mathbf{X}}_{t}\widehat{\mathbf{X}}_{t}^{T}\Big{)}-\frac{1}{\sigma_{a}}\widetilde{\mathbf{S}}_{0}\Big{\|}_{\max}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}+\sqrt{\frac{\log p}{p}}\Bigg{)}.$ (A.17)

In addition, by (6.3),

$\frac{\operatorname{tr}(\widetilde{\mathbf{S}}_{0})}{p}=\frac{\operatorname{tr}(\bm{\Sigma}_{0})}{p}+O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}\Bigg{)}=1+O_{p}\Bigg{(}\sqrt{\frac{\log p}{n}}\Bigg{)}.$

The desired result (6.9) follows. $\hfill\Box$

##  Proof of Lemma 8

The desired bound (6.10) follows from (6.9) and (A.17). The bound (6.11) then
follows from (6.10), Assumption 3, Lemma 6 and Weyl’s Theorem.

Finally, about (6.12), by Assumption 3, (6.10), (6.11) and the sin$\theta$
Theorem,

$\|\widehat{\bm{\Gamma}}_{K}-\widetilde{\bm{\Gamma}}_{K}\|_{2}=O_{p}\Bigg{(}\sqrt{\frac{\log
p}{n}}+\sqrt{\frac{\log p}{p}}\Bigg{)}.$ (A.18)

By the definitions of $\widehat{\bm{\Gamma}}_{K}$ and
$\widetilde{\bm{\Gamma}}_{K}$, we have
$\widehat{\phi}_{k}=(1/\widehat{\lambda}_{0;k})\widehat{\bm{\Sigma}}_{0}\widehat{\phi}_{k}$
and
$\widetilde{\phi}_{k}=(1/\widetilde{\lambda}_{0;k})\widetilde{\mathbf{S}}_{0}\widetilde{\phi}_{k}$.
Hence, for all $1\leq k\leq K$,

$\displaystyle\|\widehat{\phi}_{k}-\widetilde{\phi}_{k}\|_{\max}\leq$
$\displaystyle\Bigg{|}1-\frac{\widehat{\lambda}_{0;k}}{\widetilde{{\lambda}}_{0;k}}\Bigg{|}\cdot\|\widehat{\phi}_{k}\|_{\max}$
(A.19)
$\displaystyle+\frac{1}{\widetilde{\lambda}_{0;k}}\|(\widehat{\bm{\Sigma}}_{0}-\widetilde{\mathbf{S}}_{0})\widehat{\phi}_{k}\|_{\max}$
$\displaystyle+\frac{1}{\widetilde{\lambda}_{0;k}}\|\widetilde{\mathbf{S}}_{0}(\widehat{\phi}_{k}-\widetilde{\phi}_{k})\|_{\max}.$

Note that

$\|(\widehat{\bm{\Sigma}}_{0}-\widetilde{\mathbf{S}}_{0})\widehat{\phi}_{k}\|_{\max}\leq\sqrt{p}\|\widehat{\bm{\Sigma}}_{0}-\widetilde{\mathbf{S}}_{0}\|_{\max}\|\widehat{\phi}_{k}\|_{2},$

and

$\|\widetilde{\mathbf{S}}_{0}(\widehat{\phi}_{k}-\widetilde{\phi}_{k})\|_{\max}\leq\sqrt{p}\|\widetilde{\mathbf{S}}_{0}\|_{\max}\|\widehat{\phi}_{k}-\widetilde{\phi}_{k}\|_{2}.$

By Assumption 3, Lemma 6, (6.10) and (6.11), we have
$\|\widehat{\bm{\Sigma}}_{0}\|_{\max}=O_{p}(1)$,
$\|\widehat{\bm{\Sigma}}_{0}\|_{\max}=O_{p}(1)$, and there exists $c>0$ such
that with probability approaching one, for all large $p$,

$\widehat{\lambda}_{0;K}>cp,\quad\text{
and}\quad\widetilde{\lambda}_{0;K}>cp.$ (A.20)

Because
$\widehat{\bm{\Gamma}}_{K}=\widehat{\bm{\Sigma}}_{0}\widehat{\bm{\Gamma}}_{K}\widehat{\bm{\Lambda}}_{0K}^{-1}$,
it follows that

$\|\widehat{\bm{\Gamma}}_{K}\|_{\max}\leq\sqrt{p}\|\widehat{\bm{\Sigma}}_{0}\|_{\max}\|\widehat{\bm{\Gamma}}_{K}\|_{2}\|\widehat{\bm{\Lambda}}_{0K}^{-1}\|_{2}=O_{p}\Bigg{(}\frac{1}{\sqrt{p}}\Bigg{)}.$
(A.21)

The desired bound (6.12) follows from (6.10), (6.11) and (A.18)–(A.21).
$\hfill\Box$