Higher Moment Estimation for Elliptically-distributed Data: Is it Necessary to Use a Sledgehammer to Crack an Egg?

TL;DR

This paper introduces a simple yet effective method for estimating moments of high-dimensional elliptical distributions without complex precision matrix estimation, demonstrating asymptotic efficiency and practical utility in finance.

Contribution

The paper proposes the marginal aggregation estimator (MAE) and block-wise aggregation estimator (BAE) for efficient high-dimensional moment estimation without precision matrix estimation.

Findings

MAE achieves asymptotic efficiency under mild sparsity conditions.

MAE and BAE perform well in simulations and financial data applications.

The methods simplify high-dimensional elliptical distribution analysis.

Abstract

Multivariate elliptically-contoured distributions are widely used for modeling economic and financial data. We study the problem of estimating moment parameters of a semi-parametric elliptical model in a high-dimensional setting. Such estimators are useful for financial data analysis and quadratic discriminant analysis. For low-dimensional elliptical models, efficient moment estimators can be obtained by plugging in an estimate of the precision matrix. Natural generalizations of the plug-in estimator to high-dimensional settings perform unsatisfactorily, due to estimating a large precision matrix. Do we really need a sledgehammer to crack an egg? Fortunately, we discover that moment parameters can be efficiently estimated without estimating the precision matrix in high-dimension. We propose a marginal aggregation estimator (MAE) for moment parameters. The MAE only requires estimating the diagonal of covariance matrix and is convenient to implement. With mild sparsity on the covariance structure, we prove that the asymptotic variance of MAE is the same as the ideal plug-in estimator which knows the true precision matrix, so MAE is asymptotically efficient. We also extend MAE to a block-wise aggregation estimator (BAE) when estimates of diagonal blocks of covariance matrix are available. The performance of our methods is validated by extensive simulations and an application to financial returns.

Figures

Tables

Table 1. Table 1 : Variance comparison of estimators (known ( 𝝁 , 𝚺 ) 𝝁 𝚺 (\boldsymbol{\mu},\boldsymbol{\Sigma}) ; ∗∗ means the constant is optimal).

	IE	ME	MAE	BAE
dominating term	$\frac{1}{n}{\left(\frac{r_{2m}}{r_m^2}-1\right)}^{\ast \ast }$	$\frac{1}{n}\left(\frac{r_{2m}{\eta }_{2m}}{r_m^2{\eta }_m^2}-1\right)$	$\frac{1}{n}{\left(\frac{r_{2m}}{r_m^2}-1\right)}^{\ast \ast }$	$\frac{1}{n}{\left(\frac{r_{2m}}{r_m^2}-1\right)}^{\ast \ast }$
2nd-order term	$\frac{1}{np}\frac{r_{2m}}{r_m^2}{h}_m{\left(p\right)}^{\ast \ast }$	—	$\frac{1}{np}\frac{r_{2m}}{r_m^2}{h}_m\left(1\right)$	$\frac{1}{np}\frac{r_{2m}}{r_m^2}{h}_m\left(k\right)$
correlation term	$0^{\ast \ast }$	$0^{\ast \ast }$	$\frac{C}{np}{\sum }_{1\le i\ne j\le p}{\left|{\Lambda }_{jj}\right|}^2$	$\frac{C}{np}{\sum }_{I\ne J\in 𝒜}{\Vert {𝚲}_{I,J}\Vert }_F^2$

Table 2. Table 2 : Empirical coverage probability that θ 2 subscript 𝜃 2 \theta_{2} lies in the 95 % percent 95 95\% confidence interval by Equation ( 11 ) for data following multivariate Gaussian or multivariate t 𝑡 t -distributions, across a variety of settings.

$n=100$	$p=$ 250	400	550	700	850	1000
Gaussian	92.0%	95.0%	93.5%	95.5%	95.5%	96.5%
Student’s $t$	96.5%	98.0%	94.5%	97.0%	96.0%	96.5%
$p=500$	$n=$ 50	100	150	200	250	300
Gaussian	95.5%	94.2%	93.5%	93.0%	95.5%	94.0%
Student’s $t$	98.0%	96.0%	95.5%	93.5%	94.5%	97.0%

Table 3. Table 3 : Pairwise correlations of ξ ^ t subscript ^ 𝜉 𝑡 \widehat{\xi}_{t} across GCIS sectors. Numbers ≥ .45 absent .45 \geq.45 are marked in circles.

	E	CD	CO	CS	F	HC	IN	IT	M	R	U
Energy (E)	–	.36	.43	.42		.37	.44	.35
Consumer Discretionary (CD)	.36	–	.30	.35	.43	.34	.37	.33	.33	.36	.35
Communication Services (CO)	.43	.30	–	.33	.43	.33	.40	.32	.39	.44	.41
Consumer Staples (CS)	.42	.35	.33	–	.42	.36	.37	.35	.38	.38	.43
Financials (F)		.43	.43	.42	–	.42		.40
Health Care (HC)	.37	.34	.33	.36	.42	–	.40	.39	.37	.36	.33
Industrials (IN)	.44	.37	.40	.37		.40	–	.39		.44
Information Technology (IT)	.35	.33	.32	.35	.40	.39	.39	–	.36	.33	.34
Materials (M)		.33	.39	.38		.37	.49	.36	–	.42	.43
Real Estate (R)		.36	.44	.38		.36	.44	.33	.42	–
Utilities (U)		.34	.41	.43		.33	.45	.34	.43	.46	–

Equations

$$\mathbf{Y}=\mathbf{\mu }+\xi {\mathbf{\Sigma }}^{1/2}\mathbf{U},$$

$$\mathbb{E}({\xi }^2)=p.$$

$${\theta }_m\equiv p^{-m}\mathbb{E}({\xi }^{2m}).$$

$${\hat{\theta }}_m^{\mathrm{I}}(\mathbf{\mu },\mathbf{\Omega })=\frac{1}{n{p}^m}\sum _{i=1}^n({\xi }_i^2{)}^m=\frac{1}{n{p}^m}\sum _{i=1}^n\{({\mathbf{Y}}_i-\mathbf{\mu }{)}^{\mathrm{T}}\mathbf{\Omega }({\mathbf{Y}}_i-\mathbf{\mu }){\}}^m.$$

$${\hat{\theta }}_m^{\mathrm{I}}(\hat{\mathbf{\mu }},\hat{\mathbf{\Omega }})=\frac{1}{n{p}^m}\sum _{i=1}^n\{({\mathbf{Y}}_i-\hat{\mathbf{\mu }}{)}^{\mathrm{T}}\hat{\mathbf{\Omega }}({\mathbf{Y}}_i-\hat{\mathbf{\mu }}){\}}^m,$$

$$Y_j={\mu }_j+\xi {\tilde{U}}_j,j=1,\dots ,p.$$

$$\mathbb{E}({\tilde{U}}_j^{2m})=p^{-m}{c}_m{\sigma }_{jj}^m,\text{mbox}where{c}_m=(2m-1)!!(p/2{)}^m\frac{\Gamma (p/2)}{\Gamma (p/2+m)}.$$

$${\hat{\theta }}_{m,j}^{\mathrm{M}}({\mu }_j,{\sigma }_{jj})=\frac{1}{n{p}^m}\sum _{i=1}^n\frac{(Y_{ij}-{\mu }_j{)}^{2m}}{\mathbb{E}({\tilde{U}}_j^{2m})}=\frac{1}{c_m{\sigma }_{jj}^m}\frac{1}{n}\sum _{i=1}^n(Y_{ij}-{\mu }_j{)}^{2m}.$$

$${\hat{\theta }}_m^{\mathrm{M}}(\mathbf{\mu },\mathrm{diag}(\mathbf{\Sigma }))=\frac{1}{p}\sum _{j=1}^p{\hat{\theta }}_{m,j}^{\mathrm{M}}({\mu }_j,{\sigma }_{jj})=\frac{1}{c_{m}np}\sum _{j=1}^p\left\{\frac{1}{{\sigma }_{jj}^m}\sum _{i=1}^n(Y_{ij}-{\mu }_j{)}^{2m}\right\}.$$

$${\hat{\theta }}_m^{\mathrm{M}}(\hat{\mathbf{\mu }},\mathrm{diag}(\hat{\mathbf{\Sigma }}))=\frac{1}{c_{m}np}\sum _{j=1}^p\left\{\frac{1}{{\hat{\sigma }}_{jj}^m}\sum _{i=1}^n(Y_{ij}-{\hat{\mu }}_j{)}^{2m}\right\},$$

$${\mathbf{Y}}_t={\mathbf{\mu }}_t+\mathbf{B}{\mathbf{f}}_t+{\xi }_t{\mathbf{\Sigma }}_t^{1/2}{\mathbf{U}}_t,t=1,\cdots ,T,$$

$$|{\hat{\theta }}_m^{\mathrm{M}}(\hat{\mathbf{\mu }},\mathrm{diag}(\hat{\mathbf{\Sigma }}))-{\theta }_m|=O_{\mathbb{P}}(n^{-1/2}).$$

$$\bigl{|}\widehat{\theta}^{\,\mathrm{M}}_{m}(\widehat{\boldsymbol{\mu}},\mathrm{diag}(\widehat{\boldsymbol{\Sigma}}))-\theta_{m}\bigr{|}\leq C_{\epsilon}\bigl{(}n^{-1/2}+\max\{\alpha_{n},\beta_{n}\}\bigr{)}.$$

$$\frac{\mathop{\rm var}\bigl{(}\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\mathrm{diag}(\boldsymbol{\Sigma}))\bigr{)}}{\theta_{m}^{2}}\leq\frac{1}{n}\frac{r_{2m}-r_{m}^{2}}{r_{m}^{2}}+\frac{1}{np}\frac{r_{2m}}{r^{2}_{m}}\frac{\eta_{2m}-\eta_{m}^{2}}{\eta_{m}^{2}}+\frac{C_{m}}{n}\frac{r_{2m}}{r^{2}_{m}\eta_{m}^{2}}\frac{\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}}{p^{2}}.$$

$$\frac{\mathop{\rm var}\bigl{(}\widehat{\theta}^{\,\mathrm{I}}_{m}(\boldsymbol{\mu},\boldsymbol{\Omega})\bigr{)}}{\theta_{m}^{2}}=\frac{1}{n}\frac{r_{2m}-r_{m}^{2}}{r_{m}^{2}}+\frac{1}{np}\frac{r_{2m}}{r_{m}^{2}}2m^{2}\bigl{[}1+O(p^{-1})\bigr{]}.$$

$$\frac{\mathop{\rm var}\bigl{(}\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\mathrm{diag}(\boldsymbol{\Sigma}))\bigr{)}}{\mathop{\rm var}\bigl{(}\widehat{\theta}^{\,\mathrm{I}}_{m}(\boldsymbol{\mu},\boldsymbol{\Omega})\bigr{)}}\to 1.$$

$$\frac{\sqrt{n}\Bigl{(}\widehat{\theta}^{\,\mathrm{M}}_{m,j}(\widehat{\mu}_{j},\widehat{\sigma}_{jj})-\theta_{m}\Bigr{)}}{\sqrt{\frac{c_{2m}}{c_{m}^{2}}\widehat{\theta}_{2m}-\widehat{\theta}_{m}^{2}}}\to_{d}N(0,1),$$

$${\hat{\theta }}_{m,j}^{\mathrm{M}}\pm \frac{q_{1-\alpha /2}}{\sqrt{n}}\sqrt{\frac{c_{2m}}{c_m^2}{\hat{\theta }}_{2m}-{\hat{\theta }}_m^2},$$

$$\frac{\mathop{\rm var}(\widehat{\theta}^{\,\mathrm{M}}_{m})}{\theta_{m}^{2}}\sim\frac{\mathop{\rm var}\bigl{(}\widehat{\theta}^{\,\mathrm{I}}_{m})}{\theta_{m}^{2}}\sim\frac{1}{n}\frac{r_{2m}-r_{m}^{2}}{r_{m}^{2}},\qquad\frac{\mathop{\rm var}(\widehat{\theta}^{\,\mathrm{M}}_{m,j})}{\theta_{m}^{2}}\sim\frac{1}{n}\frac{(\eta_{2m}/\eta_{m}^{2})r_{2m}-r_{m}^{2}}{r_{m}^{2}}.$$

$${\mathbf{Y}}_J\stackrel{(d)}{=}{\mathbf{\mu }}_J+B^{1/2}\xi \cdot {\mathbf{\Sigma }}_{JJ}^{1/2}{\mathbf{U}}_K,$$

$${\xi }^{2m}=\frac{\{({\mathbf{Y}}_J-{\mathbf{\mu }}_J{)}^{\mathrm{T}}{\mathbf{\Sigma }}_{JJ}^{-1}({\mathbf{Y}}_J-{\mathbf{\mu }}_J){\}}^m}{B^m}.$$

$$c_{1,K}^{\ast }=K,c_{m,K}^{\ast }=p\times \frac{K+2m-2}{p+2m-2}\times c_{m-1,K}^{\ast }\text{mbox}form\ge 2.$$

$${\hat{\theta }}_{m,J}^{\mathrm{B}}({\mathbf{\mu }}_J,{\mathbf{\Sigma }}_{JJ})$$

$$=\frac{1}{n{c}_{m,K}^{\ast }}\sum _{i=1}^n\{({\mathbf{Y}}_{i,J}-{\mathbf{\mu }}_J{)}^{\mathrm{T}}{\mathbf{\Sigma }}_{JJ}^{-1}({\mathbf{Y}}_{i,J}-{\mathbf{\mu }}_J){\}}^m.$$

$$\widehat{\theta}^{\,\mathrm{B}}_{m}\bigl{(}\boldsymbol{\mu},\;\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma})\bigr{)}=\frac{1}{|{\cal A}|}\sum_{J\in{\cal A}}\widehat{\theta}^{\,\mathrm{B}}_{m,J}(\boldsymbol{\mu}_{J},\boldsymbol{\Sigma}_{JJ}).$$

$$|{\hat{\theta }}_m^{\mathrm{B}}(\hat{\mathbf{\mu }},{\mathrm{diag}}_{\mathcal{A}}(\hat{\mathbf{\Sigma }}))-{\theta }_m|=O_{\mathbb{P}}(n^{-1/2}).$$

$$\bigl{|}\widehat{\theta}^{\,\mathrm{B}}_{m}(\widehat{\boldsymbol{\mu}},\mathrm{diag}_{\cal A}(\widehat{\boldsymbol{\Sigma}}))-\theta_{m}\bigr{|}\leq C_{\epsilon}\bigl{(}n^{-1/2}+\max\{\alpha_{n},\beta_{n}\}\bigr{)}.$$

$$\frac{\mathop{\rm var}\bigl{(}\widehat{\theta}_{m}^{\,\mathrm{B}}\bigl{(}\boldsymbol{\mu},\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma})\bigr{)}\bigr{)}}{\theta_{m}^{2}}\leq\frac{1}{n}\frac{r_{2m}-r_{m}^{2}}{r_{m}^{2}}+\frac{1}{np}\frac{r_{2m}}{r_{m}^{2}}\bar{h}_{m}({\cal A})+\frac{C_{m}}{n}\frac{1}{|{\cal A}|^{2}}\sum_{\begin{subarray}{c}I,J\in{\cal A}\\ I\neq J\end{subarray}}\|\boldsymbol{\Sigma}_{II}^{-1/2}\boldsymbol{\Sigma}_{IJ}\boldsymbol{\Sigma}_{JJ}^{-1/2}\|^{2}.$$

$${\bar{h}}_m(\mathcal{A})=h_m(k)=\frac{k\cdot \mathrm{var}({\chi }_k^{2m})}{(\mathbb{E}{\chi }_k^{2m}{)}^2}.$$

$$|\widetilde{\Sigma}(i,j)|/\sqrt{\widetilde{\Sigma}(i,i)\widetilde{\Sigma}(j,j)}>t,\qquad\mbox{for $1\leq i<j\leq p$}.$$

$$\mathcal{A}=\{\text{mbox}allcomponentsof{\mathcal{G}}_t\}.$$

Full text

For low-dimensional elliptical models, efficient moment estimators can be obtained by plugging in an estimate of the precision matrix. Natural generalizations of the plug-in estimator to high-dimensional settings perform unsatisfactorily, due to estimating a large precision matrix. Do we really need a sledgehammer to crack an egg? Fortunately, we discover that moment parameters can be efficiently estimated without estimating the precision matrix in high-dimension.

We propose a marginal aggregation estimator (MAE) for moment parameters. The MAE only requires estimating the diagonal of covariance matrix and is convenient to implement. With mild sparsity on the covariance structure, we prove that the asymptotic variance of MAE is the same as the ideal plug-in estimator which knows the true precision matrix, so MAE is asymptotically efficient. We also extend MAE to a block-wise aggregation estimator (BAE) when estimates of diagonal blocks of covariance matrix are available. The performance of our methods is validated by extensive simulations and an application to financial returns.

1 Introduction

The classical multivariate statistics is largely motivated by relaxing the Gaussian assumption, which is not satisfied in many applications. There is an extensive literature in finance on the tail-index estimates of stock returns; while being unimodal and symmetric, the empirical returns exhibit leptokurtosis, which means that they have heavier tails and flatter peaks than those of normal data (Fama, 1965; Bollerslev and Wooldridge, 1992; Eberlein and Keller, 1995; Frahm et al., 2003; Cizek et al., 2005). Empirical evidence of the violation of Gaussian assumption has also been observed in genomics (Liu et al., 2003; Posekany et al., 2011; Hardin and Wilson, 2009) and in bioimaging (Ruttimann et al., 1998). The family of multivariate elliptically contoured distributions (Kelker, 1970), which we shall call elliptical distributions in short, provides a natural generalization of multivariate Gaussian distributions. Recently, many statistical methods for elliptically distributed data have been proposed, including works on covariance matrix estimation (Fan et al., 2018), graphical modeling (Han and Liu, 2012), classification (Fan et al., 2015b), etc.

The elliptical distributions are typically used as a semi-parametric model. Given a mean vector $\boldsymbol{\mu}=(\mu_{1},\ldots,\mu_{p})^{\rm T}\in\mathbb{R}^{p}$, a covariance matrix $\boldsymbol{\Sigma}=(\sigma_{jk})_{1\leq j,k\leq p}\in\mathbb{R}^{p\times p}$ and a probability characteristic function $\phi:[0,\infty)\to\mathbb{R}$, we say a random vector $\bm{Y}=(Y_{1},\ldots,Y_{p})^{\rm T}$ has an elliptical distribution $\mathcal{E}(\boldsymbol{\mu},\boldsymbol{\Sigma},\phi)$ if

[TABLE]

where $\bm{U}$ is a random vector that is uniformly distributed on the unit sphere $\mathbb{S}^{p-1}$, and independent of $\bm{U}$, $\xi$ is a nonnegative random variable whose characteristic function is $\phi$. For model identifiability, we normalize $\xi$ such that

[TABLE]

Under (1)-(2), $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are the mean vector and covariance matrix of $\bm{Y}$, respectively. The variable $\xi$ determines which sub-family the distribution belongs to. When $\xi^{2}$ is a chi-square random variable, it belongs to the multivariate Gaussian sub-family, and when $\xi^{2}$ follows an $F$-distribution, it belongs to the multivariate $t$ sub-family or multivariate Cauchy sub-family. For most applications, the sub-family of the elliptical distribution is unknown, leaving the distribution of $\xi$ unspecified.

Although full knowledge of the distribution of $\xi$ is often not required, an estimate of its moment parameters is useful to statistical analysis and for understanding the tail of the distributions. One application is in quadratic classification. When data from two classes both follow elliptical distributions but have unequal covariance matrices, Fan et al. (2015b) showed that an estimate of $\mathbb{E}(\xi^{4})$ is desired for building a quadratic classifier. Another application is to capture the tail behavior of financial returns by estimating the leptokurtosis. Modeling the returns of a set of financial assets by an elliptical distribution, the leptokurtosis equals to $\{p(p+2)\}^{-1}\mathbb{E}(\xi^{4})-1$, so the problem reduces to estimating $\mathbb{E}(\xi^{4})$.

For any $m\geq 1$, define the $m$-th scaled even moment of $\xi$ by

[TABLE]

The first scaled even moment $\theta_{1}$ is $1$. In this paper, we are interested in estimating $\theta_{m}$ for any fixed $m\geq 2$, given independent and identically distributed (i.i.d.) samples $\bm{Y}_{1},\cdots,\bm{Y}_{n}$ from (1).

1.1 The plug-in estimators

We consider an ideal case where $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are known. Given $iid$ samples $\bm{Y}_{1},\cdots,\bm{Y}_{n}$ from an unknown elliptical distribution, each $\bm{Y}_{i}$ has a decomposition $\bm{Y}_{i}=\boldsymbol{\mu}+\xi_{i}\boldsymbol{\Sigma}^{1/2}\bm{U}_{i}$, and $\xi_{1},\ldots,\xi_{n}$ are $iid$ copies of $\xi$. Using the fact that $\bm{U}_{i}$ takes values on the unit sphere, we observe $\xi_{i}^{2}=(\bm{Y}_{i}-\boldsymbol{\mu})^{\rm T}\boldsymbol{\Omega}(\bm{Y}_{i}-\boldsymbol{\mu})$ for $i=1,\ldots,n$, where $\boldsymbol{\Omega}\equiv\boldsymbol{\Sigma}^{-1}$. Hence, in the ideal case, $\xi_{1},\ldots,\xi_{n}$ are directly observed. It motivates the following estimator of $\theta_{m}$:

[TABLE]

We call $\widehat{\theta}_{m}^{\,\mathrm{I}}(\boldsymbol{\mu},\boldsymbol{\Omega})$ the Ideal Estimator. The ideal estimator is not feasible in practice, and a natural modification is to plug in estimates of $(\boldsymbol{\mu},\boldsymbol{\Omega})$. This gives rise to the plug-in estimator:

[TABLE]

This estimator was proposed by Maruyama and Seo (2003) in the setting of a fixed dimension, where they used the sample mean to estimate $\boldsymbol{\mu}$ and the inverse of the sample covariance matrix to estimate $\boldsymbol{\Omega}$. In the modern high-dimensional settings where $p$ grows with $n$, one can on longer use the inverse of sample covariance matrix to estimate $\boldsymbol{\Omega}$; Fan et al. (2015b) proposed plugging in an estimator of $\boldsymbol{\Omega}$ from high-dimensional sparse precision matrix estimation methods, with stringent structural assumptions on $\boldsymbol{\Omega}$.

However, the plug-in estimators perform unsatisfactorily for high-dimensional settings due to the difficulty of estimating $\boldsymbol{\Omega}$. Existing methods of estimating $\boldsymbol{\Omega}$ only perform well under stringent conditions, such as the sub-Gaussian assumption on the distribution and/or structural assumptions on $\boldsymbol{\Omega}$ (e.g., sparsity). Especially, the structural assumption on $\boldsymbol{\Omega}$ is critical for the success of these methods. Figure 1 shows the performance of the plug-in estimator when the structural assumption required by $\widehat{\boldsymbol{\Omega}}$ is violated. We consider two estimators of $\boldsymbol{\Omega}$, the CLIME estimator (Cai et al., 2011) which requires sparsity of $\boldsymbol{\Omega}$, and the POET estimator (Fan et al., 2013) which assumes a factor structure with sparse covariance of the idiosyncratic component. On the left panel of Figure 1, we generate elliptical data with a sparse covariance matrix, $\boldsymbol{\Sigma}_{i,j}=a^{|i-j|}$, $1\leq i,j\leq p$, where $a$ controls the sparsity level and varies in $\{0.5,0.55,\ldots,0.85,0.9\}$. Here, the structural assumption of POET is not satisfied, and the associated plug-in estimator of $\theta_{2}$ performs unsatisfactorily. On the right panel, we generate data with a sparse precision matrix $\boldsymbol{\Omega}$, where each entry of the upper triangle of $\boldsymbol{\Omega}$ has a probability of $a$ to be nonzero,444We generate $\boldsymbol{\Omega}$ using fastclime.generator($\cdot$) in the R package clime, where the graph argument is set “random”. with $a$ chosen from $\{0.5,0.55,\ldots,0.85,0.9\}$. The assumption of CLIME is violated, so the associated plug-in estimator of $\theta_{2}$ has a unsatisfactory performance.

In fact, the philosophy of plug-in estimators is problematic. Estimating large precision matrices is a well-known difficult problem (even for Gaussian data), as one needs to estimate a large number of parameters. On the other hand, our problem only involves estimating one single parameter $\theta_{m}$. Intuitively, the latter should be much easier than the former. The plug-in estimators are realy using “a sledgehammer to crack an egg.”

1.2 The marginal aggregation estimator (MAE)

Is it possible to avoid using the “sledgehammer” of precision matrix estimation? We show that this is possible by a new marginal aggregation estimator. In model (1), letting $\widetilde{U}_{j}$ be the $j$-th coordinate of $\widetilde{\bm{U}}\equiv\boldsymbol{\Sigma}^{1/2}\bm{U}$, we have

[TABLE]

Our key observation is that each individual coordinate of $\bm{Y}$ contains information of $\xi$. It motivates us to construct an estimator of $\theta_{m}$ using only one coordinate of samples. Let $\sigma_{jj}$ be the $j$-th diagonal of $\boldsymbol{\Sigma}$. We notice that (6) implies $\xi^{2m}=(Y_{j}-\mu_{j})^{2m}/\widetilde{U}_{j}^{2m}$. The random variable $\widetilde{U}_{j}$ is unobserved, but its distribution is known once $\sigma_{jj}$ is given. It can be shown that (see Proposition 3.1)

[TABLE]

Inspired by (6)-(7), we introduce an estimator of $\theta_{m}$ using the marginal data $Y_{1j},\ldots,Y_{nj}$:

[TABLE]

We call $\widehat{\theta}^{\,\mathrm{M}}_{m,j}(\mu_{j},\sigma_{jj})$ the Marginal Estimator. It only requires knowledge of $(\mu_{j},\sigma_{jj})$ and successfully avoids precision matrix estimation. For each $1\leq j\leq p$, we can define a marginal estimator and we will show that all marginal estimator contains the same amount of information about $\theta_{m}$ (see Theorem 2.4). All these marginal estimators are unbiased, so taking their average gives rise to a new unbiased estimator:

[TABLE]

We call $\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\mathrm{diag}(\boldsymbol{\Sigma}))$ the Marginal Aggregation Estimator (MAE). The “aggregation” of marginal estimators helps reduce the asymptotic variance. Our proposed estimator is a natural plug-in version of (9) given by

[TABLE]

where $c_{m}$ is as in (7), $\widehat{\boldsymbol{\mu}}$ is an estimator of $\boldsymbol{\mu}$, and $\{\widehat{\sigma}_{jj}\}_{j=1}^{p}$ are the estimators of $\{\sigma_{jj}\}_{j=1}^{p}$.

Compared with the plug-in estimator (5), MAE is numerically more appealing, as it only needs to estimate the diagonal entries of $\boldsymbol{\Sigma}$. Back to the example in Figure 1, we implement MAE using sample mean as $\widehat{\boldsymbol{\mu}}$ and sample covariance matrix as $\widehat{\boldsymbol{\Sigma}}$. MAE significantly outperforms the plug-in estimators, even when the structural assumptions of the plug-in estimators are satisfied.

1.3 Organization of the paper

In Section 2, we study the theoretical properties of MAE. Under mild regularity conditions, we show that MAE is unbiased and root-$n$ consistent, regardless of the structure of $\boldsymbol{\Sigma}$. We also show that MAE is asymptotically efficient, with an asymptotic variance matching that of the ideal estimator when $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are given. We also discuss how to construct a confidence interval of $\theta_{m}$.

In Section 3, we generalize the idea of MAE to develop estimators of $\theta_{m}$ that use a small subset of the coordinates. We introduce the block-wise estimator and the blockwise aggregation estimator (BAE), analogous to the marginal estimator and MAE. These ideas help further reduce the estimation errors in the second order.

Section 4 validates the theoretical insight by extensive simulations. Section 5 gives an application of MAE to time series data. We consider an extension of model (1) to multivariate time series:

[TABLE]

where $\boldsymbol{\mu}_{t}$ is the time-varying mean, $\bm{f}_{t}\in\mathbb{R}^{K}$ is a vector of $K$ observed factors, and ${\bf B}$ is a $p\times K$ matrix of factor loadings. We extend MAE to a method for estimating the realized $\xi_{t}$. Its application to stock returns provides a new index that captures information of whole market. Section 6 contains conclusions and discussions. All the proofs are relegated to the appendix.

Notation: Throughout this paper, for any vector ${\bf v}$ and matrix ${\bf M}$, we let $\|{\bf v}\|$ denote the Euclidean norm of ${\bf v}$ and let $\|{\bf M}\|$, $\|{\bf M}\|_{F}$ and $\|{\bf M}\|_{\max}$ denote its spectral norm, Frobenius norm and entry-wise maximum norm, respectively. We use $\widehat{\theta}^{\,\mathrm{M}}_{m,j}(\mu_{j},\sigma_{jj})$, $\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\mathrm{diag}(\boldsymbol{\Sigma}))$, $\widehat{\theta}^{\,\mathrm{I}}_{m}(\boldsymbol{\mu},\boldsymbol{\Omega})$ and $\widehat{\theta}^{\,\mathrm{B}}_{m}(\boldsymbol{\mu},\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma}))$ to denote the Marginal Estimator, MAE, Ideal Estimator, and BAE (to be introduced), respectively, with given $(\boldsymbol{\mu},\boldsymbol{\Sigma})$; when $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are replaced by $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$, it means we plug in estimators of the mean vector and covariance matrix. We frequently use notations $(\theta_{m},c_{m},\eta_{m},r_{m})$, where $\theta_{m}$ is defined in (3), $c_{m}$ is defined in (7), $\eta_{m}$ and $r_{m}$ are defined in Definition 2.1. For all settings in this paper, $\eta_{m}$ is a constant, $(\theta_{m},c_{m},r_{m})$ depend on $p$ but are at the constant scale.

2 Theoretical properties of MAE

We study the asymptotic properties of MAE defined in (10), assuming both $(n,p)$ tend to infinity. First, we study the consistency of MAE. The following theorem shows that, when the distribution is marginally sub-Gaussian, if we plug in the sample mean and sample covariance matrix as $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$, then MAE is always root-$n$ consistent.

Theorem 2.1** (Root-$n$ consistency).**

Under model (1), suppose $\log^{2}(p)=o(n)$ and $\max_{1\leq j\leq p}\|Y_{j}-\mu_{j}\|_{\psi_{2}}\leq C$, where $\|\cdot\|_{\psi_{2}}$ denotes the sub-Gaussian norm.555For a random variable $X$, its sub-Gaussian norm is defined as $\|X\|_{\psi_{2}}=\sup_{k\geq 1}k^{-1}(\mathbb{E}|X|^{k})^{1/k}$. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the MAE in (10), where $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ are the sample mean vector and sample covariance matrix. Then,

[TABLE]

The root-$n$ consistency of MAE requires no conditions on either $\boldsymbol{\Sigma}$ or $\boldsymbol{\Omega}$. It confirms our previous insight that estimating moment parameters is an “easier” statistical problem than estimating large matrices. On the other hand, the plug-in estimators only perform well when the assumed structural assumptions (e.g., sparsity) on $\boldsymbol{\Sigma}$ or $\boldsymbol{\Omega}$ are satisfied.

Many distributions in the elliptical family are heavy-tailed and don’t satisfy the marginal sub-Gaussianity assumption. In these cases, we prefer to use robust estimators of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ (Fan et al., 2017; Sun et al., 2018+). They are M-estimators with robust loss functions or rank-based estimators. Compared to the sample mean and sample covariance estimators, these robust estimators lead to sharper bounds of $\|\widehat{\boldsymbol{\mu}}-\boldsymbol{\mu}\|_{\infty}$ and $\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\|_{\max}$ in the case of heavy-tailed data. The next theorem studies MAE with general mean/covariance estimators.

Theorem 2.2** (Consistency, with general mean/covariance estimators).**

Under model (1), suppose $\log^{2}(p)=o(n)$ and $\theta_{2m}\leq C$. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the MAE in (10). We assume the estimators $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ satisfy $\max_{1\leq j\leq p}|\widehat{\mu}_{j}-\mu_{j}|\leq\alpha_{n}$ and $\max_{1\leq j\leq p}|\widehat{\sigma}_{jj}-\sigma_{jj}|\leq\beta_{n}$ with probability $1-o(1)$, where $\alpha_{n}\to 0$ and $\beta_{n}\to 0$ as $n,p\to\infty$. Then, for any $\epsilon>0$, with probability $1-\epsilon$, there is a constant $C_{\epsilon}>0$ such that

[TABLE]

The typical error rate of robust estimators is $\alpha_{n}\asymp\sqrt{\log(p)/n}$ and $\beta_{n}\asymp\sqrt{\log(p)/n}$ (Fan et al., 2017; Sun et al., 2018+), so the associated MAE satisfies $|\widehat{\theta}^{\,\mathrm{M}}_{m}-\theta_{m}|=O_{\mathbb{P}}(\sqrt{\log(p)/n})$. Compared with the rate in Theorem 2.1, the extra $\sqrt{\log(p)}$ factor here is a price paid for heavy tails.

Next, we study the asymptotic variance of MAE. By Theorem 2.1, MAE is already rate-optimal. We would like to see whether it also achieves the optimal “constant”. We shall compare its asymptotic variance with that of the Ideal Estimator (4). Since the Ideal Estimator knows the true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$, for a fair comparison, we consider MAE with true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

Definition 2.1**.**

For any $k\geq 1$, let $\eta_{k}=\mathbb{E}[N(0,1)^{2k}]$ and $r_{k}=(\mathbb{E}\xi^{2k})/(\mathbb{E}\chi_{p}^{2k})$, where $\chi^{2}_{p}$ denotes the chi-square distribution with $p$ degrees of freedom.

The quantities $r_{k}$ capture the difference between moments of an elliptical distribution and moments of a multivariate Gaussian distribution with matching mean and covariance matrix. It depends on $p$ but is at the constant scale under our settings.

Theorem 2.3** (Variance).**

Under model (1), suppose $\log^{2}(p)=o(n)$ and $\theta_{2m}\leq C$. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the MAE in (9) where $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are given. Let $\boldsymbol{\Lambda}=[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}\boldsymbol{\Sigma}[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}$ be the correlation matrix. There is a constant $C_{m}>0$, independent of the distribution of $\xi$, such that

[TABLE]

When $m=2$, the equality holds with $C_{m}=72$.

The upper bound for the variance has three terms: The first term is $O(n^{-1})$; as we shall see, this term matches with the variance of the benchmark estimator. The second term is $O(n^{-1}p^{-1})$ and is negligible for diverging $p$. The third term is caused by correlations among different marginal estimators $\widehat{\theta}^{\,\mathrm{M}}_{m,j}$. This term is negligible as long as $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$; consider a special case where $\|\boldsymbol{\Sigma}\|$ is bounded, then $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=O(p)$; so the requirement of $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$ is mild. Indeed, if requires that the sparsity of correlation coefficients: $\sum_{i\not=j}\lambda_{ij}=o(p^{2})$, where $\boldsymbol{\Lambda}=(\lambda_{ij})$. The next proposition confirms that the asymptotic variance of MAE is the same as the asymptotic variance of the Ideal Estimator:

Proposition 2.1** (Comparison with benchmark).**

Let $\{\bm{Y}_{i}\}_{i=1}^{n}$ be iid samples of model (1). Suppose $\theta_{2m}\leq C$. For the Ideal Estimator in (4),

[TABLE]

As a result, if $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$, where $\boldsymbol{\Lambda}=[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}\boldsymbol{\Sigma}[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}$ is the correlation matrix, then

[TABLE]

Last, we construct confidence intervals of $\theta_{m}$. Since MAE is the average of $p$ strongly dependent marginal estimators, its asymptotic normality is hard to approach. We instead use the marginal estimator in (8) to construct confidence intervals.

Theorem 2.4** (Asymptotic normality).**

Under model (1), suppose $\log^{2}(p)=o(n)$ and $\max_{1\leq j\leq p}\|Y_{j}-\mu_{j}\|_{\psi_{2}}\leq C$, where $\|\cdot\|_{\psi_{2}}$ denotes the sub-Gaussian norm. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the Marginal Estimator in (8) for an arbitrary $1\leq j\leq p$, where $(\widehat{\mu}_{j},\widehat{\sigma}_{jj})$ are the sample mean and sample variance of $\{Y_{ij}\}_{i=1}^{n}$. Then,

[TABLE]

where $c_{k}=(2k-1)!!\,(p/2)^{k}\frac{\Gamma(p/2)}{\Gamma(p/2+k)}$ for $k\geq 1$, and $(\widehat{\theta}_{2m},\widehat{\theta}_{m})$ are consistent estimators of $(\theta_{2m},\theta_{m})$.

This theorem shows somewhat surprisingly that all marginal estimator contains the same amount of information about $\theta_{m}$. Given consistent estimators $(\widehat{\theta}_{2m},\widehat{\theta}_{m})$, the asymptotic level-$\alpha$ confidence interval of $\theta_{m}$ is

[TABLE]

where $q_{1-\alpha/2}$ is the $(1-\alpha/2)$-quantile of a standard normal. It doesn’t matter which of $1\leq j\leq p$ we use, as these marginal estimators have the same asymptotic variance. For the estimators $(\widehat{\theta}_{2m},\widehat{\theta}_{m})$, we suggest using MAE.

If we only need a point estimator but not a confidence interval, we prefer MAE to the Marginal Estimator, as MAE has a smaller variance in many scenarios. For example, when $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$, by plugging in the true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$,

[TABLE]

Since $\eta_{2m}>\eta_{m}^{2}$, the latter variance is strictly larger. In contrast, MAE is first-order efficient.

3 Extension to blockwise aggregation

In the construction of MAE, each marginal estimator only uses one coordinate of the samples. It is convenient to implement and gives rise to an estimator that is first-order efficient, provided that the third term in Theorem 2.3 is negligible. It turns out that, the second order term in the variance can be improved upon by using blockwise aggregation, and so is the third term, which is related to the correlation structure. Our simulation studies below show that the improvement is real. This motivates us to extend the marginal estimator to a blockwise estimator that uses a small number of coordinates of the samples and takes into account their correlation structures. We then generalize MAE to BAE — an aggregation of many blockwise estimators.

BAE can be applied to settings where the covariance matrix is approximately blockwise diagonal after row/column permutation. Figure 2 gives such an example, where the S&P 500 stocks divide into many small-size blocks according to sectors or industries of stocks and the stock returns within each block are correlated but admits block structure after taking out the market factor. BAE can take advantage of the within-block correlations and further improve MAE in the second order term.

3.1 A block-wise aggregation estimator (BAE)

We fix a block $J\subset\{1,2,\ldots,p\}$ and let $K=|J|$. For any vector ${\bf v}\in\mathbb{R}^{p}$ and matrix ${\bf M}\in\mathbb{R}^{p\times p}$, let ${\bf v}_{J}$ be the subvector of ${\bf v}$ containing the coordinates indexed by $J$ and let ${\bf M}_{JJ}$ be the submatrix of ${\bf M}$ containing the entries indexed by $J\times J$. By Fang and Zhang (1990), when $\bm{Y}$ follows an elliptical distribution (1), the subvector $\bm{Y}_{J}$ satisfies that

[TABLE]

where $B$ is a random variable that follows a beta distribution Beta$\,(\frac{K}{2},\frac{p-K}{2})$, the random vector ${\bm{U}}_{K}$ follows a uniform distribution on the unit sphere $\mathbb{S}^{K-1}$, and $(\xi,B,{\bm{U}}_{K})$ are mutually independent. Since $\|{\bm{U}}_{K}\|=1$,

[TABLE]

The random variable $B$ is not directly observable, but its expectation is known:

Proposition 3.1**.**

For each $m\geq 1$ and $1\leq K\leq p$, define $c^{*}_{m,K}=p^{m}\mathbb{E}(B^{m})$ with $B\sim{\rm Beta}\,(\frac{K}{2},\frac{p-K}{2})$. Then,

[TABLE]

Replacing $B^{m}$ by its expectation, we immediately have an estimator of $\theta_{m}$ based on $\{\bm{Y}_{i,J}\}_{i=1}^{n}$:

[TABLE]

We call $\widehat{\theta}^{\,\mathrm{B}}_{m,J}(\boldsymbol{\mu}_{J},\boldsymbol{\Sigma}_{JJ})$ the Blockwise Estimator. Now, given a collection of blocks ${\cal A}=\{J_{1},J_{2},\ldots,J_{N}\}$, we can define a blockwise estimator for each $J\in{\cal A}$ and then take their average:

[TABLE]

We call \widehat{\theta}^{\,\mathrm{B}}_{m}\bigl{(}\boldsymbol{\mu},\;\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma})\bigr{)} the Blockwise Aggregation Estimator (BAE). Here $\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma})$ denotes the collection of diagonal blocks $\boldsymbol{\Sigma}_{JJ}$ with $J\in{\cal A}$. Our final estimator is a plug-in version of BAE by plugging in an estimator $\widehat{\boldsymbol{\mu}}$ and estimators of those diagonal blocks of $\boldsymbol{\Sigma}$.

Since BAE only estimates the small-size diagonal blocks of $\boldsymbol{\Sigma}$ and does not need to estimate $\boldsymbol{\Omega}$, it inherits a nice property of MAE: root-$n$ consistency is guaranteed with no conditions on $\boldsymbol{\Sigma}$ or $\boldsymbol{\Omega}$.

Theorem 3.1** (Root-$n$ consistency).**

Fix $m\geq 2$ and $K\geq 1$. Under model (1), suppose $\log^{2m}(p)=o(n)$ and $\max_{1\leq j\leq p}\|Y_{j}-\mu_{j}\|_{\psi_{2}}\leq C$. We assume the minimum eigenvalue of any $K\times K$ diagonal block of $\boldsymbol{\Sigma}$ is lower bounded by $C$. Let ${\cal A}$ be a collection of nonrandom, non-overlapping blocks such that the size of each block is bounded by $K$. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the BAE in (15), where $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are estimated by the sample mean vector and sample covariance matrix. Then,

[TABLE]

Theorem 3.2** (Consistency, with general mean/covariance estimators).**

Fix $m\geq 2$ and $K\geq 1$. Under model (1), we assume $\log^{2m}(p)=o(n)$, $\theta_{2m}\leq C$, and the minimum eigenvalue of any $K\times K$ diagonal block of $\boldsymbol{\Sigma}$ is lower bounded by $C$. Let ${\cal A}$ be a collection of nonrandom, non-overlapping blocks where the size of blocks is bounded by $K$. Given iid samples $\{\bm{Y}_{i}\}_{i=1}^{n}$, consider the BAE in (15), where $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ satisfy $\|\widehat{\boldsymbol{\mu}}-\boldsymbol{\mu}\|_{\infty}\leq\alpha_{n}$ and $\max_{J\in{\cal A}}\|\widehat{\boldsymbol{\Sigma}}_{JJ}-\boldsymbol{\Sigma}_{JJ}\|\leq\beta_{n}$ with probability $1-o(1)$, with $\alpha_{n}\to 0$ and $\beta_{n}\to 0$ as $n,p\to\infty$. Then, for any $\epsilon>0$, with probability $1-\epsilon$, there is a constant $C_{\epsilon}>0$ such that

[TABLE]

We note that MAE is a special case of BAE, with all block size equal to $1$. The motivation of generalizing MAE to BAE is to better take advantage of correlation structures, and this is revealed by comparing the asymptotic variances of two methods; see Section 3.2 below. To implement BAE, we need to determine the collection of blocks, and in Section 3.3 we discuss how to select blocks.

3.2 Variance comparison

We compute the asymptotic variance of BAE and compare it with the asymptotic variances of MAE and Ideal Estimator. Same as before, in the variance calculation we assume $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are given.

Definition 3.1**.**

For each $k\geq 1$, let $h_{m}(k)=\frac{k\cdot\mathrm{var}(\chi_{k}^{2m})}{(\mathbb{E}\chi_{k}^{2m})^{2}}$, where $\chi^{2}_{k}$ denotes the chi-square distribution with $k$ degrees of freedom. Given a collection of blocks ${\cal A}$, let $\bar{h}_{m}({\cal A})=\frac{p}{|{\cal A}|^{2}}\sum_{J\in{\cal A}}\frac{h_{m}(|J|)}{|J|}$.

Theorem 3.3** (Variance of BAE).**

Let $\{\bm{Y}_{i}\}_{i=1}^{n}$ be $iid$ samples of model (1). Fix $m\geq 2$ and suppose $\theta_{2m}\leq C$. There exists a constant $\widetilde{C}_{m}>0$, independent of the distribution of $\xi$, such that for any collection ${\cal A}$ of non-overlapping blocks,

[TABLE]

The upper bound of the variance has three terms:

• •

The first term is $O(n^{-1})$, which also appears in the variance of MAE and Ideal Estimator. It is the dominating term of the variance.

• •

The second term is $O(p^{-1}n^{-1})$, where the constant in front of it is related to a quantity $\bar{h}_{m}({\cal A})$. We call $\bar{h}_{m}({\cal A})$ the block-division factor, as it is only a function of ${\cal A}$. To see how this factor changes with block size, let’s consider a special case where all blocks have an equal size $k$ and $p$ is a multiple of $k$. Then,

[TABLE]

It is a monotone decreasing function of $k$ (see Figure 3). Hence, increasing the block size leads to a reduction of this term, which indicates that the second order efficiency of MAE can be improved with $m>1$.

• •

The last term comes from the correlations among estimators associated with different blocks. It doesn’t exist for the Ideal Estimator, but both MAE and BAE have this extra term. For MAE, all off-diagonal entries of $\boldsymbol{\Sigma}$ contribute to this term. However, for BAE, only off-diagonal blocks contribute. Especially, when $\boldsymbol{\Sigma}$ is blockwise diagonal with respect to ${\cal A}$, this extra term becomes zero. Again, increasing the block size leads to a reduction of this term.

From MAE to BAE, we can see that the dominating term in the variance bound remains the same, but the other two terms are reduced and the performance still improves. However, we cannot use too large blocks, because BAE needs to invert an estimate of $\boldsymbol{\Sigma}_{J,J}$ and the error of estimating $\widehat{\boldsymbol{\Sigma}}_{J,J}$ increases as the block size increases.

We now give a more thorough comparison of four estimators, the Ideal Estimator (IE) $\widehat{\theta}^{\,\mathrm{I}}_{m}$, the Marginal Estimator (ME) $\widehat{\theta}^{\,\mathrm{M}}_{m,j}$, the MAE $\widehat{\theta}^{\,\mathrm{M}}_{m}$, and the BAE $\widehat{\theta}^{\,\mathrm{B}}_{m}$; see Table 1. We conclude that

• •

IE has the optimal variance, but it works unsatisfactorialy in the real case of unknown $(\boldsymbol{\mu},\boldsymbol{\Sigma})$, as it requires estimating $\boldsymbol{\Omega}$.

• •

ME avoids estimating $\boldsymbol{\Omega}$ and works in the real case, but its asymptotic variance is non-optimal.

• •

MAE aggregates a number of ME’s and achieves the optimal variance when $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$.

• •

Compared with MAE, BAE relaxes the condition of $\|\boldsymbol{\Lambda}-{\bf I}\|_{F}^{2}=o(p^{2})$ and reduces the second-order term of the variance.

From ME to BAE, we have used two methodological ideas: to aggregate “local” estimators and to use a block of coordinates in each “local” estimator. Both help reduce the variance of the estimator, with the first idea playing a more significant role.

Remark 1. IE and MAE are special cases of BAE with equal-size blocks of $k=1$ and $k=p$, respectively. We note that $h_{m}(1)=\frac{\eta_{2m}-\eta_{m}^{2}}{\eta_{m}^{2}}$ and $h_{m}(p)=2m^{2}[1+O(p^{-1})]$, so Theorem 3.3 matches with the variance bounds of MAE (Theorem 2.3) and the IE (Proposition 2.1).

Remark 2 (multivariate Gaussian). Let’s consider a special case where the data are multivariate Gaussian but the user doesn’t know and still applies the estimators in this paper. For Gaussian distributions, the first term in the variance bound disappears, so the estimators considered here all have a faster rate of convergence as $O(p^{-1}n^{-1})$. This is the only case where a large $p$ helps, i.e., “dimensionality is a blessing.” Moreover, the difference between MAE and BAE is more prominent, as the second term in the variance bound is now dominating. Figure 4 displays the error bound according to Theorem 3.3 for the case of $\boldsymbol{\Sigma}={\bf I}$ and $\boldsymbol{\Sigma}$ being a blockwise diagonal matrix with $2\times 2$ blocks whose off-diagonal element is $\rho$. The results favor BAE, especially for the blockwise $\boldsymbol{\Sigma}$ with large within-block off-diagonals.

3.3 Construction of blocks

We provide two approaches of selecting the blocks. The first approach works well when the true $\boldsymbol{\Sigma}$ is approximately block-wise diagonal, such as example on the returns of the S&P 500 components(see Figure 2). The second approach is a random scheme and works for general settings.

BAE1: Constructing blocks from a raw estimate of $\boldsymbol{\Sigma}$. Let $\widetilde{\boldsymbol{\Sigma}}$ be a raw estimate of $\boldsymbol{\Sigma}$; for example, it can be the sample covariance matrix or the robust estimator of $\boldsymbol{\Sigma}$ in Section 4. Fixing a threshold $t\in(0,1)$, we define a graph ${\cal G}_{t}$ with nodes $\{1,2\cdots,p\}$, where there is an undirected edge between nodes $i$ and $j$ if and only if the estimated absolute correlation exceeds $t$, namely,

[TABLE]

The nodes of this graph uniquely partitions into components (a component of a graph is a maximal connected subgraph). We propose using

[TABLE]

See Figure 5 for an illustration of this procedure.

This approach guarantees that all blocks are non-overlapping. Numerical evidence suggests that it performs well with an appropriate choice of $t$, especially when the true $\boldsymbol{\Sigma}$ is blockwise diagonal. However, the threshold $t$ is a tuning parameter, and it can be inconvenient to select $t$ in a data-driven fashion. Below, we introduce a tuning-free approach.

BAE2: Randomly selecting pairs as blocks. In this approach, we let

[TABLE]

This approach is designed for block size equal to $2$, and the obtained blocks may overlap. Although it sounds ad-hoc, this approach has an appealing numerical performance. When the number of pairs are sampled sufficiently large, by the law of large numbers, it approaches the all pairwise aggregation estimator and this explains why the approach has an appealing numerical performance. This approach can easily be extended to blocks of any size that is smaller than $n$ so long as the estimated covariance matrix for each block can be easily inverted and estimated well.

4 Simulations

We investigate the performance of estimators on extensive simulations. To have realistic simulation settings, we use a $\boldsymbol{\Sigma}$ calibrated from stock returns. The calibration procedure is the same as that in Fan et al. (2015c) and Fan et al. (2013). Fix $p$. We take the daily returns of $p$ companies in S&P 500 index with the largest market capitalization from July 1st, 2013 to June 29th, 2018 (data were downloaded from the COMPUSTAT website). We fit the Fama-French three-factor model to the excess returns $\{\bm{Y}_{t}\}_{t=1}^{T}$:

[TABLE]

where ${\bf B}\in\mathbb{R}^{p\times 3}$ is the factor loading matrix, $\bm{f}_{t}\in\mathbb{R}^{3}$ denotes the Fama-French factors with covariance matrix $\mathop{\rm cov}(\bm{f}_{t})\in\mathbb{R}^{3\times 3}$ and ${\bf u}_{t}$ is the idiosyncratic component. This factor model induces a covariance structure for ${\bf Y}_{t}$:

[TABLE]

where $\boldsymbol{\Sigma}_{u}$ is the covariance matrix of idiosyncratic noise ${\bf u}_{t}$. We downloaded the factors $\{\bm{f}_{t}\}_{t=1}^{T}$ from the Kenneth French data library and used the method in Fan et al. (2013) with the recommended threshold (for estimating sparse $\boldsymbol{\Sigma}_{u}$) to get and estimate $\widehat{\boldsymbol{\Sigma}}_{Y}$. We then use $\widehat{\boldsymbol{\Sigma}}_{Y}$ as the true $\boldsymbol{\Sigma}$ to generate data from model (1).

When implementing the estimators, we plug in two different estimators of $(\boldsymbol{\mu},\boldsymbol{\Sigma})$. The first choice is to use sample mean and sample covariance matrix. The second choice is to use robust M-estimators, called adaptive Huber estimator (Fan et al., 2017; Sun et al., 2018+), which are designed for heavy-tailed data. These estimators lead to better large-deviation bounds. In detail, for a tuning parameter $\tau>0$ chosen by cross-validation, we estimate $\boldsymbol{\mu}$ by $(\widehat{\mu}_{1},\ldots,\widehat{\mu}_{p})^{{\rm T}}$, where

[TABLE]

the Huber loss. We estimate $\boldsymbol{\Sigma}$ by $(\widehat{\sigma}_{jk})_{1\leq j,k\leq p}$, where

[TABLE]

Here, each tuning parameter $\tau_{jk}$ is selected via cross-validation using the data $\{(Y_{ij},Y_{ik})\}_{i=1}^{n}$.

Experiment 1: Performance of MAE.

Fix $m=2$. We consider four sub-experiments:

• •

Experiments 1.1 and 1.3: We fix $p=500$ and let $n$ vary in $\{50,100,150,200,250,300\}$. The data follow multivariate Gaussian distributions (Experiment 1.1) or multivariate $t$-distributions with degrees of freedom equal to $4.5$ (Experiment 1.3).

• •

Experiments 1.2 and 1.4: We fix $n=100$ and let $p$ vary in $\{250,400,550,700,850,1000\}$. The data follow multivariate Gaussian distributions (Experiment 1.2) or multivariate $t$-distributions with degrees of freedom equal to $4.5$ (Experiment 1.4).

In all settings, $p>n$, so we focus on the challenging case of high-dimensionality. For each setting, we compare four estimators:

• •

$\widehat{\theta}^{\,\mathrm{I}}(\boldsymbol{\mu},\boldsymbol{\Sigma})$: Ideal Estimator, which knowns $(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

• •

$\widehat{\theta}^{\,\mathrm{M}}(\boldsymbol{\mu},\boldsymbol{\Sigma})$: MAE with given $(\boldsymbol{\mu},\boldsymbol{\Sigma})$.

• •

$\widehat{\theta}^{\,\mathrm{M}}(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$: MAE, where $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are estimated using the sample mean/covariance matrix in Experiment 1.1&1.2 and using the aforementioned robust-M estimators for Experiment 1.3&1.4.

• •

$\widehat{\theta}^{\,\mathrm{I}}(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}}_{P})$: Plug-in Ideal Estimator, with plugged-in estimators of $(\boldsymbol{\mu},\boldsymbol{\Sigma})$. We use the sample mean to estimate $\widehat{\boldsymbol{\mu}}$ and use POET (Fan et al., 2013) (with a default threshold) to estimate $\boldsymbol{\Sigma}$.

The results are presented in Figure 6, where the $y$-axis is $\log\{(\widehat{\theta}_{2}/\theta_{2}-1)^{2}\}$, based on the average over $200$ repetitions. As we have expected, the Ideal Estimator always gives the lowest error, however, such an estimator is not practically feasible. Instead, we plug estimates of $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ into the Ideal Estimator to make it practically feasible, then it has an unsatisfactory performance; this confirms our previous insight about the drawback of the plug-in estimator. Our proposed MAE works well, always significantly better than the plug-in estimator. The performance of MAE becomes better as the sample size $n$ grows, and its performance stays relatively stable as the dimension $p$ grows. This is desirable: our proposed estimator doesn’t face any curse of dimensionality. The results are similar for the multivariate Gaussian data and the multivariate $t$-data, except that for Gaussian data, MAE with $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ even outperforms MAE with true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$. One possible reason is the self-normalization phenomenon: An estimator, when divided by its sample variance, gives better performance than that divided by the true variance.

Experiment 2: Confidence Interval.

For each of the experiments above: Experiments 1.1, 1.2, 1.3 and 1.4, we calculate the probability that the true value of $\theta_{2}$ lies in the confidence interval derived in Theorem 2.4 and presented in Equation (11). In Table 2, we see that for a $95\%$ confidence interval, the empirical coverage probabilities are close to the confidence level.

Experiment 3: Performance of BAE.

We study whether BAE, which uses a block of coordinates at a time and takes advantage of the correlation structure, can further improve the performance of MAE. The four sub-experiments, Experiments 3.1-3.4, have the same settings as those of Experiments 1.1-1.4. When implementing BAE, we use the second approach in Section 3.3 to choose the blocks; note that the blocks all have a size $2$ and may overlap. We use the sample mean/covariance to estimate $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ for multivariate Gaussian data and the robust M-estimators for multivariate $t$ data. Since we focus on the comparison between MAE and BAE, we do not report the errors of the Ideal Estimator and plug-in estimator in this experiment.

The results are presented in Figure 7. First, we can see that BAE improves the performance of MAE, especially when $p$ is large. Second, the self-normalization phenomenon is also observed: BAE with $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ even outperforms BAE with true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$, especially for Gaussian data.

5 Application: Estimating realized $\xi_{t}$ in a time series

Given the returns of a panel of stocks, we are interested in extending the idea of MAE to provide a daily risk index for the whole panel of stocks. We cast it as the problem of estimating the realized $\xi_{t}$ in a multivariate time series with elliptically-distributed noise. Let $\bm{Y}_{1},\ldots,\bm{Y}_{T}\in\mathbb{R}^{p}$ be the returns of $p$ stocks during a time period of $T$ days. We extend model (1) to an elliptical model for multivariate time series

[TABLE]

where $\boldsymbol{\mu}_{t}$ is the time-varying mean, $\bm{f}_{t}\in\mathbb{R}^{K}$ is a vector of $K$ factors, and ${\bf B}$ is a $p\times K$ matrix of factor loadings. We are interested in estimating the daily realized $\xi_{t}$.

Our method has four steps:

1. Estimate $\boldsymbol{\mu}_{t}$. For daily or higher frequency data, we set $\widehat{\boldsymbol{\mu}}_{t}\equiv{\bf 0}$, since it is commonly believed that the short-time returns are not predictable. For weekly or monthly data, we estimate $\boldsymbol{\mu}_{t}$ by the weekly or monthly average.

2. Obtain the factor-adjusted returns $\widehat{\bm{Z}}_{t}$. Let $\widehat{\bm{f}}_{t}\in\mathbb{R}^{K}$ contain either observed factors or data-drive factors from PCA (Fan et al., 2013). We then follow the approach in Fan et al. (2013) to get $\widehat{{\bf B}}$, the estimated factor loading matrix. Let

[TABLE]

3. Estimate $\boldsymbol{\Sigma}_{t}$. We assume $\boldsymbol{\Sigma}_{t}$ is a diagonal matrix and estimate its diagonal elements by fitting an ARCH model on each coordinate of $\bm{Z}_{t}$. In detail, for each $1\leq j\leq p$, let $Z_{t}(j)$ be the $j$-th coordinate of $\bm{Z}_{t}$. We assume there is idiosyncratic noise $\{\epsilon_{t}(j)\}_{t=1}^{T}$ such that

[TABLE]

where $k$ is the order of ARCH model and $(a_{0},\ldots,a_{k})$ are parameters. We estimate $(a_{0},\ldots,a_{k})$ using the conditional maximum likelihood estimator and then construct $\{\widehat{\lambda}_{t}(j)\}_{t=1}^{T}$. Let

[TABLE]

4. Estimate $\xi_{t}$. We adapt the idea of MAE to the current setting. Let $\bm{Z}_{t}=\bm{Y}_{t}-{\bf B}\bm{f}_{t}$. Our model becomes $\bm{Z}_{t}=\xi_{t}\boldsymbol{\Sigma}_{t}^{1/2}\bm{U}_{t}$, i.e., the $j$-th component of $\bm{Z}_{t}$ is $Z_{t}(j)=\xi_{t}(\boldsymbol{\Sigma}^{1/2}\bm{U}_{t})_{j}$. It follows that

[TABLE]

where $\Sigma_{t}(j,j)$ is the $j$-th diagonal of $\boldsymbol{\Sigma}_{t}$. Here, the last equality is due to $c_{1}=1$ in Equation (7). We approximate $(\bm{Z}_{t},\boldsymbol{\Sigma}_{t})$ by $(\widehat{\bm{Z}}_{t},\widehat{\boldsymbol{\Sigma}}_{t})$ and get a marginal estimator of $\xi_{t}^{2}$: $\widehat{\xi}^{2}_{t,j}=p\widehat{Z}_{t}(j)/\widehat{\Sigma}_{t}(j,j)$. We then aggregate them:

[TABLE]

In Section D of the appendix, we investigate the performance of our estimator in simulations. Under a variety of settings, our estimated curve of $\widehat{\xi}_{t}$ fits the true curve of $\xi_{t}$ very well. See details therein.

We applied our estimator to the S&P 500 stock returns. We took the daily returns of $300$ stocks from the S&P 500 index with the largest market capitalization, from July 1, 2008 to June 29, 2012. Each stock is assigned a Global Industry Classification Standard (GCIS) code. The GCIS code divides $300$ stocks into eleven sectors: Energy, Consumer Discretionary, Communication Services, Consumer Staples, Financials, Health Care, Industrials, Information Technology, Materials, Real Estate, and Utilities. We applied our estimator to stocks in each sector. When implementing our method, we set $\widehat{\boldsymbol{\mu}}_{t}\equiv{\bf 0}$ in Step 1, used three observed Fama-French factors as $\widehat{\bm{f}}_{t}$ in Step 2, and set the order of ARCH model to $k=2$ in Step 3.

The curves of estimated $\widehat{\xi}_{t}$ for $11$ sectors are displayed in Figure 8 (the curves are smoothed by taking a moving average on a $65$-day window). The estimated $\widehat{\xi}_{t}$ for all sectors largely synchronize, reaching their peaks during the 2008 financial crisis. In the crisis, the estimated $\widehat{\xi}_{t}$ for the Financials sector is significantly larger than that of other sectors. The large value of $\widehat{\xi}_{t}$ for the Financials sector remains in the post-crisis period until May, 2009. We also computed the pairwise correlations among $\widehat{\xi}_{t}$ of $11$ sectors, as shown in Table 3. It suggests that the $\widehat{\xi}_{t}$ for the Energy sector and the Financials sector are highly correlated with each other. These two sectors are also highly correlated with sectors of Materials, Real Estate, and Utilities. In comparison, for the Consumer Discretionary sector and Information Technology sector, their $\widehat{\xi}_{t}$ are less correlated with those of other sectors.

6 Discussion

In this paper, we consider the problem of estimating the even moments of $\xi$ in an elliptical distribution $\bm{Y}=\boldsymbol{\mu}+\xi\boldsymbol{\Sigma}^{1/2}\bm{U}$. A natural idea is the plug-in estimator (Maruyama and Seo, 2003; Fan et al., 2015b), which requires an estimator $\widehat{\boldsymbol{\Omega}}$ of the precision matrix and whose performance crucially relies on structural assumptions on $\boldsymbol{\Omega}$ or $\boldsymbol{\Sigma}$. Instead, we propose a marginal aggregation estimator (MAE) that only needs to estimate the diagonal of $\boldsymbol{\Sigma}$. Our approach validates the insight that estimating a large precision matrix is statistically more challenging than estimating a moment parameter—it is unnecessary to use the sledge hammer to crack an egg. We prove that MAE is root-$n$ consistent, under no conditions on $\boldsymbol{\Sigma}$ or $\boldsymbol{\Omega}$. We also show that MAE achieves the first-order efficiency, with an asymptotic variance matching with the variance of an ideal estimator when $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are given. We further generalize MAE to a block-wise aggregation estimator (BAE) that needs to estimate small-size diagonal blocks of $\boldsymbol{\Sigma}$. BAE takes advantage of correlations among coordinates and improves MAE on the second-order efficiency. Our proposed estimators are conceptually simple and easy to implement.

Estimating the moment parameters of an elliptical distribution is useful in quadratic discriminant analysis (Fan et al., 2015b) and estimating tail behavior of financial returns (Fama, 1965; Bollerslev and Wooldridge, 1992; Eberlein and Keller, 1995; Frahm et al., 2003; Cizek et al., 2005). In an application on the stock returns, we propose a multivariate time series model with factor structures and elliptically distributed idiosyncratic noise. We extend MAE to an estimator for estimating the day-to-day value of $\xi_{t}$. We apply the method to stocks of each industry sector. It produces an “tail index” for each industry sector. These tail indices reveal interesting difference among industry sectors, especially during the financial crisis.

The study leaves a few open questions for future work. The first is how to improve the estimators for heavy tailed data. Our current approach plugs into MAE the robust estimators of mean and covariance matrix. Instead, we may construct a robust M-estimator for simultaneously estimating $(\theta_{m},\mu_{j},\sigma_{jj})$ with marginal data and then aggregate these marginal estimators of $\theta_{m}$ in a similar way. We hope such an approach helps remove the $\sqrt{\log(p)}$-factor in the error rate of Theorem 2.2. The second is the optimal strategy of constructing blocks in BAE. There is a trade-off in choosing the blocks: With larger blocks, it reduces the variance of the estimator when true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are plugged in, but at the same time, the errors of estimating diagonal blocks of $\boldsymbol{\Sigma}$ increase. How to construct the blocks in a data-driven way is an interesting question. Third, the current theory for BAE assumes non-overlapping blocks. The results can be extended to overlapping blocks, with nontrivial efforts. We leave it for future work. The last problem is to extend our estimators to time dependent data, where the distribution of $\xi$ have change-points. For financial data, such change-points may relate to financial boom or crisis. We propose a kernel-smoothed version of MAE: Given data $\{\bm{Y}_{t}\}_{t=1}^{n}$, for a kernel function $K_{h}(\cdot)$ with bandwidth $h$, let

[TABLE]

We can similarly define the one-sided versions of the kernel estimator. We can combine these estimators with change-point detection methods, which we leave for future work.

Appendix A Proof of main results

A.1 Proof of Theorem 2.1

Write for short $\widehat{\theta}^{\,\mathrm{M}}_{m}=\widehat{\theta}^{\,\mathrm{M}}_{m}(\widehat{\boldsymbol{\mu}},\mathrm{diag}(\widehat{\boldsymbol{\Sigma}}))$ and $\widetilde{\theta}^{\,\mathrm{M}}_{m}=\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\mathrm{diag}(\boldsymbol{\Sigma}))$. By Theorem 2.3, $\widetilde{\theta}_{m}^{\,\mathrm{M}}$ is unbiased and satisfies

[TABLE]

We note that $(\eta_{m},C_{m})$ are constants, $(\theta_{m},r_{m},r_{2m})$ are bounded above/below by constants, and all entries of the correlation matrix $\boldsymbol{\Lambda}$ are bounded by $1$. Hence, the right hand side is $O(n^{-1})$, and it implies

[TABLE]

To show the claim, it suffices to show that

[TABLE]

Below, we show (18). Write for short $X_{ij}=(Y_{ij}-\mu_{j})/\sqrt{\sigma_{jj}}$, for $1\leq i\leq n,1\leq j\leq p$. For any $k\geq 0$, let $S_{kj}=\frac{1}{n}\sum_{i=1}^{n}X^{k}_{ij}$. Using these notations,

[TABLE]

At the same time, noticing that $(\widehat{\mu}_{j}-\mu_{j})/\sqrt{\sigma_{jj}}=S_{1j}$ and $(Y_{ij}-\widehat{\mu}_{j})/\sqrt{\sigma_{jj}}=X_{ij}-S_{1j}$, we have

[TABLE]

Combining the above gives

[TABLE]

To bound the right hand side of (19), we define an event. By (6), $Y_{ij}=\mu_{j}+\xi_{i}(\boldsymbol{\Sigma}^{1/2}\bm{U}_{i})_{j}$. Let $\boldsymbol{\Lambda}=[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}\boldsymbol{\Sigma}[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}$. Then, $X_{ij}=\frac{Y_{ij}-\mu_{j}}{\sqrt{\sigma_{jj}}}=\xi_{i}(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})_{j}$. It follows that

[TABLE]

Note that $\mathbb{E}X_{ij}=(\mathbb{E}\xi_{i}^{k})\mathbb{E}[(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})^{k}_{j}]$. At the same time, since $X_{ij}\sim N(0,1)$ when $\xi_{i}^{2}\sim\chi_{p}^{2}$, it holds that $\mathbb{E}[N^{k}(0,1)]=(\mathbb{E}\chi_{p}^{k})\mathbb{E}[(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})^{k}_{j}]$. Together, we have $\mathbb{E}(X^{k}_{ij})=\mathbb{E}[N^{k}(0,1)]\cdot[(\mathbb{E}\xi_{i}^{k})/(\mathbb{E}\chi_{p}^{k})]$. Our assumption of $\theta_{2m}\leq C$ guarantees $(\mathbb{E}\xi_{i}^{k})/(\mathbb{E}\chi_{p}^{k})\leq C$ for $1\leq k\leq 4m$. It follows that $\mathbb{E}(X_{ij}^{k})\leq C$ and $\mathrm{var}(X_{ij}^{k})\leq C$ for $1\leq k\leq 2m$. As a result,

[TABLE]

Using the marginal sub-Gaussianity, for any $\epsilon>0$, there exists a constant $C>0$ such that, with probability $\geq 1-\epsilon$,

[TABLE]

Let $B$ be the event that (24) holds. To show (18), it suffices to show that

[TABLE]

We now show (25). Consider $(I_{2})$ and $(I_{3})$. By (22) and using that $(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})_{j}$ has a symmetric distribution, we have $\mathbb{E}S_{kj}=0$ for any odd $k$. As a result, over the event $B$, $|S_{1j}|\leq C\sqrt{(\log p)/n}$, $|S_{(2m-1)j}|\leq C\sqrt{(\log p)/n}$ and $|S_{(2m-k)j}|\leq C$, for all $1\leq j\leq p$ and $1\leq k\leq 2m$. Additionally, since \widehat{\sigma}_{jj}/\sigma_{jj}=\frac{1}{n}\sum_{i=1}^{n}\bigl{(}\frac{Y_{ij}-\widehat{\mu}_{j}}{\sqrt{\sigma_{jj}}}\bigr{)}^{2}=\frac{1}{n}\sum_{i=1}^{n}(X_{ij}-S_{1j})^{2}=S_{2j}-S_{1j}^{2}, where $\mathbb{E}S_{2j}=1$, it holds that $\sigma_{jj}/\widehat{\sigma}_{jj}\leq C$ over the event $B$. It follows that

[TABLE]

Consider $(I_{1})$. Since $\widehat{\sigma}_{jj}/\sigma_{jj}=S_{2j}-S_{1j}^{2}$, we write

[TABLE]

Over the event $B$, $\max_{j}|S_{1j}|\leq C\sqrt{(\log p)/n}$ and $\max_{1\leq j\leq p}|\widehat{\sigma}_{jj}/\sigma_{jj}-1|\leq C\sqrt{(\log p)/n}$. This means, for all $1\leq j\leq p$, $\widehat{\sigma}_{jj}/\sigma_{jj}$ is contained in a diminishing neighborhood of $1$. We use Taylor expansion of the function $(1+x)^{-m}-1$. It gives

[TABLE]

where the third line is due to $\max_{1\leq j\leq p}|S_{1j}|\leq O(\sqrt{(\log p)/n})$ over the event $B$ and the fourth line is due to (22). By (22) and (28),

[TABLE]

Write $R_{ij}=(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})_{j}$ for short. Then,

[TABLE]

We introduce positive random variables $\{\omega_{i}\}_{i=1}^{n}$ such that $\omega_{i}^{2}\overset{iid}{\sim}\chi_{p}^{2}$ and that $\{\omega_{i}\}_{i=1}^{n}$ are independent of $\{(\xi_{i},\bm{U}_{i}):1\leq i\leq n\}$. Then, $Z_{i}\equiv\omega_{i}(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})\sim N(0,I_{p})$. For even integers $s,t$ and $1\leq j,j^{\prime}\leq p$,

[TABLE]

For all $s,t$ such that $s+t\leq 4m$, the left hand side is uniformly bounded by a constant. Additionally, by elementary probability, $\mathbb{E}(\omega_{i}^{s+t})\asymp p^{(s+t)/2}$. It follows that

[TABLE]

In particular, by taking $s=2\ell$ and $t=0$ in the above, we have $\mathbb{E}R_{ij}^{2\ell}\leq Cp^{-\ell}$ for all $0\leq\ell\leq 2m$. Additionally, $\theta_{s}=p^{-s}\mathbb{E}\xi^{2s}$ by definition, so the assumption $\theta_{2m}\leq C$ guarantees

[TABLE]

Using (36)-(37), we first bound $|\sum_{i=1}^{n}Q_{ii}|$. It is seen that

[TABLE]

As a result,

[TABLE]

We then bound $|\sum_{i\neq k}Q_{ik}|$. Consider $(i,k,i^{\prime},k^{\prime})$ such that $i\neq k$ and $i^{\prime}\neq k^{\prime}$. By (35), $\mathbb{E}Q_{ik}=0$ for $i\neq k$. Therefore, if $(i,k,i^{\prime},k^{\prime})$ are mutually distinct, $\mathbb{E}(Q_{ik}Q_{i^{\prime}k^{\prime}})=0$. It follows that

[TABLE]

By (35) and (36)-(37),

[TABLE]

Moreover, the total number of such distinct $(i,k,k^{\prime})$ is $O(n^{3})$. It follows that

[TABLE]

Pluging (38) and (39) into (32) gives

[TABLE]

We further plug (26) and (40) into (19). It gives (25). The proof is now complete. ∎

A.2 Proof of Theorem 2.2

Similar to the proof of Theorem 2.1, let $\widetilde{\theta}^{\,\mathrm{M}}_{m}$ and $\widehat{\theta}^{\,\mathrm{M}}_{m}$ denote the MAE with true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and estimates $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$; here, $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ are not necessarily the sample mean and sample covariance matrix. It follows from Theorem 2.3 that $\mathbb{E}[(\widetilde{\theta}^{\,\mathrm{M}}_{m}-\theta_{m})^{2}]\leq Cn^{-1}$. By Markov’s inequality, for any constant $C_{1}>0$,

[TABLE]

Hence, given $\epsilon>0$, we can choose an appropriate $C_{1}>0$ such that the above probability is bounded by $\epsilon/3$.

Below, we bound $|\widehat{\theta}^{\,\mathrm{M}}_{m}-\widetilde{\theta}^{\,\mathrm{M}}_{m}|$. Letting $X_{ij}=(Y_{ij}-\mu_{j})/\sqrt{\sigma_{jj}}$ and $\widehat{X}_{ij}=(Y_{ij}-\widehat{\mu}_{j})/\sqrt{\widehat{\sigma}_{jj}}$, we have

[TABLE]

where

[TABLE]

It follows that

[TABLE]

First, we consider $(J_{1})$. By direct calculations,

[TABLE]

where $S_{kj}=\frac{1}{n}\sum_{i=1}^{n}X^{k}_{ij}$ for $k\geq 0$. Under our assumption, $\max_{j}|\sigma_{jj}/\widehat{\sigma}_{jj}|\lesssim 1$, and $|\sqrt{\widehat{\sigma}_{jj}/\sigma_{jj}}-1|\leq C|\widehat{\sigma}_{jj}-\sigma_{jj}|$. Moreover, by similar technique in the proof of Theorem 2.3, we can prove that, $\frac{1}{p}\sum_{j=1}^{p}\mathbb{E}|S_{kj}|\leq C$, for $1\leq k\leq 2m$. As a result, for any $\epsilon>0$, there exists $C_{2}>0$ such that, $\frac{1}{p}\sum_{j=1}^{p}|S_{(2m)j}|\leq C_{2}$ simultaneously for $1\leq k\leq 2m$, with probability $1-\epsilon/3$. On this event,

[TABLE]

Next, we consider $(J_{2})$. By our assumption, $|\Delta_{ij}|\leq\alpha_{n}+\beta_{n}|X_{ij}|$. It follows that

[TABLE]

Plugging it into the definition of $(J_{2})$, we have

[TABLE]

Again, we can easily prove that $\frac{1}{np}\sum_{j=1}^{p}\sum_{i=1}^{n}\mathbb{E}|X_{ij}|^{k}\leq C$ for all $1\leq k\leq 2m$. It follows that, for any $\epsilon>0$, there exists $C_{3}>0$, such that $\frac{1}{np}\sum_{j=1}^{p}\sum_{i=1}^{n}|X_{ij}|^{k}\leq C_{3}$ simultaneously for all $1\leq k\leq 2m$. On this event,

[TABLE]

Combining (42)-(43) gives $|\widehat{\theta}^{\,\mathrm{M}}_{m}-\widetilde{\theta}^{\,\mathrm{M}}_{m}|\leq C\max\{\alpha_{n},\beta_{n}\}$. We further combine it with (41). It gives the claim.

A.3 Proof of Theorem 2.3

Write for short $\widehat{\theta}^{\,\mathrm{M}}_{m}=\widehat{\theta}^{\,\mathrm{M}}_{m}(\boldsymbol{\mu},\boldsymbol{\Omega})$ and $\widehat{\theta}^{\,\mathrm{M}}_{m,j}=\widehat{\theta}^{\,\mathrm{M}}_{m,j}(\mu_{j},\sigma_{jj})$. First, we show that $\widehat{\theta}^{\,\mathrm{M}}_{m}$ is unbiased. Recall that $\widehat{\theta}^{\,\mathrm{M}}_{m}=p^{-1}\sum_{j=1}^{p}\widehat{\theta}^{\,\mathrm{M}}_{m,j}$. It suffices to show $\widehat{\theta}^{\,\mathrm{M}}_{m,j}$ is unbiased for each $1\leq j\leq p$. Recall that

[TABLE]

By the form of elliptical distribution, $\bm{Y}_{i}-\boldsymbol{\mu}=\xi_{i}\widetilde{\bm{U}}_{i}$, where $\xi_{i}$ and $\widetilde{\bm{U}}_{i}$ are independent of each other. We have seen in Section 1.2 that $\mathbb{E}\widetilde{U}_{ij}^{2m}=p^{-m}c_{m}\sigma_{jj}^{m}$. It follows that

[TABLE]

Plugging it into (44) gives

[TABLE]

This proves that each $\widehat{\theta}^{\,\mathrm{M}}_{m,j}$ is unbiased. It follows that $\widehat{\theta}^{\,\mathrm{M}}_{m}$ is also unbiased.

Next, we calculate the variance of $\widehat{\theta}^{\,\mathrm{M}}_{m}$. For each $1\leq i\leq n$, let $\bm{W}^{(i)}=(W_{1}^{(i)},\ldots,W_{p}^{(i)})^{T}$, where $W_{j}^{(i)}=(Y_{ij}-\mu_{j})^{2m}/\sigma_{jj}^{m}$, $1\leq j\leq p$. Noting that $\{{\bf W}^{(i)}\}_{i=1}^{n}$ are $iid$ random vectors, we have

[TABLE]

It suffices to calculate the variance in the case of $n=1$. From now on, we fix $n=1$. Let $\bm{Y}=\boldsymbol{\mu}+\xi\bm{U}$ be the observed realization of the elliptical distribution. Write

[TABLE]

We now calculate $\mathrm{var}(W_{j})$ and $\mathrm{cov}(W_{j},W_{k})$. Recalling that $\widetilde{\bm{U}}=\boldsymbol{\Sigma}^{1/2}\bm{U}$, we define random vectors

[TABLE]

where $\chi_{p}^{2}$ is a chi-square random variable independent of $\widetilde{\bm{U}}$. Since the multivariate normal distribution is a special elliptical distribution with $\xi\sim\chi_{p}^{2}$, we immediately have $\bm{Z}\sim N({\bf 0},\boldsymbol{\Sigma})$. It follows that $\mathbb{E}\widetilde{Z}_{j}^{2m}=\sigma_{jj}^{-m}(\mathbb{E}Z_{j}^{2m})=\sigma_{jj}^{-m}(\mathbb{E}\chi_{p}^{2m})(\mathbb{E}\widetilde{U}_{j}^{2m})$. At a result, for all $m\geq 1$,

[TABLE]

It follows that

[TABLE]

Similarly, since $Y_{j}-\mu_{j}=\xi\widetilde{U}_{j}$ and $Z_{j}=\chi_{p}^{2}\widetilde{U}_{j}$, we have

[TABLE]

Therefore,

[TABLE]

Combining (48) and (51) and noting that $\tilde{Z}_{j}\sim N(0,1)$ for all $1\leq j\leq p$, we rewrite

[TABLE]

As a result,

[TABLE]

Moreover, since $\mathbb{E}\widetilde{U}_{ij}^{2m}=p^{-m}c_{m}\sigma_{jj}^{m}$ and $\mathbb{E}\widetilde{Z}_{j}^{2m}=\sigma_{jj}^{-m}(\mathbb{E}\chi_{p}^{2m})(\mathbb{E}\widetilde{U}_{j}^{2m})$, we have

[TABLE]

Plugging it into (54) gives

[TABLE]

This is for the case of $n=1$. For a general $n$, we combine it with (46) to get

[TABLE]

What remains is to calculate the variance of $\sum_{j=1}^{p}\widetilde{Z}_{j}^{2m}$. By definition,

[TABLE]

Here $\boldsymbol{\Lambda}$ coincides with the correlation matrix of the elliptical distribution. It is seen that

[TABLE]

where $\beta_{m}(\Lambda_{jk})$ denotes the covariance between $X_{1}^{2m}$ and $X_{2}^{2m}$ when $(X_{1},X_{2})^{{\rm T}}$ follows a bivariate normal distribution with covariances $\mathop{\rm var}(X_{1})=\mathop{\rm var}(X_{2})=1$ and $\mathop{\rm cov}(X_{1},X_{2})=\Lambda_{jk}$. The following lemma is proved in Section B.1:

Lemma A.1**.**

Let $\bm{X}=(X_{1},X_{2})^{T}$ be a bivariate normal random vector satisfying $\mathbb{E}(X_{1}^{2})=\mathbb{E}(X_{2}^{2})=1$ and $\mathop{\rm cov}(X_{1},X_{2})=\rho$. Let $\eta_{m}=\mathbb{E}[N(0,1)^{2m}]$ and $\beta_{m}(\rho)=\mathop{\rm cov}(X_{1}^{2m},X_{2}^{2m})$ for $m\geq 2$. Define

[TABLE]

Then, for all $m\geq 2$,

[TABLE]

As a result, $\beta_{m}(\rho)=72\rho^{2}$ for $m=2$, and $\beta_{m}(\rho)\leq C_{m}\rho^{2}$ for $m\geq 3$, where $C_{m}>0$ is a constant that only depends on $m$.

By Lemma A.1,

[TABLE]

Plugging it into (57) gives

[TABLE]

Moreover, for $m=2$, the equality holds for $C_{m}=72$. Since $\eta_{m}=3$ and $\eta_{2m}=105$, we have

[TABLE]

A.4 Proof of Proposition 2.1

Write $\widehat{\theta}^{\,{\rm I}}_{m}=\widehat{\theta}^{\,{\rm I}}_{m}(\boldsymbol{\mu},\boldsymbol{\Omega})$ for short. By definition, $\widehat{\theta}^{\,{\rm I}}_{m}=\frac{1}{np^{m}}\sum_{i=1}^{n}\xi^{2m}_{i}$, and $\theta_{m}=p^{-m}\mathbb{E}(\xi^{2m})$. Therefore

[TABLE]

We divide both sides by $\theta^{2}_{m}$ and note that $\theta_{m}=p^{-m}(\mathbb{E}\xi^{2m})=p^{-m}r_{m}(\mathbb{E}\chi_{p}^{2m})$. It follows that

[TABLE]

By elementary statistics, $\mathbb{E}\chi_{p}^{2m}=\prod_{j=0}^{m-1}(p+2j)$. As a result,

[TABLE]

Plugging (62) into (59) gives the claim.

A.5 Proof of Theorem 2.4

Fix $1\leq j\leq p$. Using the Slutsky’s lemma, we only need to prove

[TABLE]

Write for short $\widehat{\theta}^{\,\mathrm{M}}_{m,j}=\widehat{\theta}^{\,\mathrm{M}}_{m,j}(\widehat{\mu}_{j},\widehat{\sigma}_{jj})$. Let $X_{ij}=(Y_{ij}-\mu_{j})/\sqrt{\sigma_{jj}}$ and $S_{kj}=\frac{1}{n}\sum_{i=1}^{n}X^{k}_{ij}$, for $1\leq i\leq n$ and $k\geq 0$. Then, $\widehat{\mu}_{j}=S_{1j}$, $\widehat{\sigma}_{jj}=S_{2j}-S_{1j}^{2}$, and

[TABLE]

It follows that

[TABLE]

Let $\bm{S}=(S_{1j},S_{2j},\ldots,S_{(2m)j})^{{\rm T}}$. Below, we first derive the asymptotic normality of $\bm{S}$, then we use the delta method to prove (67).

First, we study the random vector $\bm{S}$. It is not hard to see that $\mathbb{E}S_{kj}=\mathbb{E}X_{ij}^{k}$. By (6), $X_{ij}=\xi_{i}(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})_{j}$, where $\{(\xi_{i},\bm{U}_{i})\}_{i=1}^{n}$ are mutually indepependent and $\boldsymbol{\Lambda}=[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}\boldsymbol{\Sigma}[\mathrm{diag}(\boldsymbol{\Sigma})]^{-1/2}$ is the correlation matrix. Since $X_{ij}\sim N(0,1)$ when $\xi_{i}\sim\chi_{p}^{2}$, the symmetry of $N(0,1)$ implies that $(\boldsymbol{\Lambda}^{1/2}\bm{U}_{i})_{j}$ has a symmetric distribution. Hence, $\mathbb{E}X_{ij}^{k}=0$ for an odd $k$. For an even $k=2s$, by definition of $c_{m}$ in (7), $\mathbb{E}[(\boldsymbol{\Lambda}^{1/2}\bm{U})^{2s}_{j}]=p^{-s}c_{s}$; also, $\mathbb{E}(\xi_{i}^{2s})=p^{s}\theta_{s}$; combining them gives $\mathbb{E}X_{ij}^{2s}=\mathbb{E}(\xi_{i}^{2s})\mathbb{E}[(\boldsymbol{\Lambda}^{1/2}\bm{U})^{2s}_{j}]=c_{s}\theta_{s}$. It follows that

[TABLE]

Moreover, $\mathop{\rm cov}(S_{kj},S_{\ell j})=\frac{1}{n}\mathop{\rm cov}(X_{ij}^{k},X_{ij}^{\ell})=\frac{1}{n}[\mathbb{E}X_{ij}^{k+\ell}-(\mathbb{E}X_{ij}^{k})(\mathbb{E}X_{ij}^{\ell})]$. It follows that

[TABLE]

By classical central limit theorem,

[TABLE]

Next, we prove (67). Define a function $h:\mathbb{R}^{2m}\to\mathbb{R}$ by $h({\bf x})=\sum_{k=0}^{2m}\gamma_{k}x_{1}^{k}x_{2m-k}$. By (71),

[TABLE]

Note that $\frac{\partial}{\partial x_{1}}h({\bf x})=\sum_{k=1}^{2m}k\gamma_{k}x_{1}^{k-1}x_{2m-k}$, and $\frac{\partial}{\partial x_{k}}h({\bf x})=\gamma_{2m-k}x_{1}^{2m-k}$ for $k\neq 1$. Combining them with (71) and (72) gives

[TABLE]

We then apply the delta method and obtain

[TABLE]

By (68), $\widehat{\theta}^{\,\mathrm{M}}_{m,j}=\frac{\sigma_{jj}^{m}}{\widehat{\sigma}_{jj}^{m}}\cdot c_{m}^{-1}h(\bm{S})$. Since $\frac{\sigma_{jj}^{m}}{\widehat{\sigma}_{jj}^{m}}\to 1$ in probability, using the Slutsky’s lemma, we have

[TABLE]

This proves (67).

A.6 Proof of Theorem 3.1

Write for short \widehat{\theta}^{\,\mathrm{B}}_{m}=\widehat{\theta}^{\,\mathrm{B}}_{m}\bigl{(}\widehat{\boldsymbol{\mu}},\;\mathrm{diag}_{\cal A}(\widehat{\boldsymbol{\Sigma}})\bigr{)} and \widetilde{\theta}^{\,\mathrm{B}}_{m}=\widehat{\theta}^{\,\mathrm{B}}_{m}\bigl{(}\boldsymbol{\mu},\;\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma})\bigr{)}. It follows from Theorem 3.3 that $\mathbb{E}[(\widetilde{\theta}^{\,\mathrm{B}}_{m}-\theta_{m})^{2}]=O(n^{-1/2})$. This implies $|\widetilde{\theta}^{\,\mathrm{B}}_{m}-\theta_{m}|=O_{\mathbb{P}}(n^{-1/2})$. Hence, it suffices to show

[TABLE]

First, we derive an expression of $\widehat{\theta}^{\,\mathrm{B}}_{m}-\widetilde{\theta}^{\,\mathrm{B}}_{m}$. Let $\bm{X}_{i,J}=\boldsymbol{\Sigma}_{J,J}^{-1/2}(\bm{Y}_{i,J}-\boldsymbol{\mu}_{J})$ and $\widehat{\bm{X}}_{i,J}=\widehat{\boldsymbol{\Sigma}}_{J,J}^{-1/2}(\bm{Y}_{i,J}-\widehat{\boldsymbol{\mu}}_{J})$ for all $1\leq i\leq n$ and $J\in{\cal A}$. Then,

[TABLE]

Let $\bm{S}_{1,J}=\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i,J}$ and $\bm{S}_{2,J}=\frac{1}{n}\sum_{i=1}^{n}\bm{X}_{i,J}\bm{X}_{i,J}^{{\rm T}}$. By direct calculations,

[TABLE]

Define an event $B$ such that

[TABLE]

It is not hard to see that the event $B$ holds with probability $1-o(1)$ (see the proof of Theorem 2.1 for similar arguments). On the event $B$, noting that $\mathbb{E}\bm{S}_{2,J}={\bf I}_{|J|}$, we have

[TABLE]

It follows that

[TABLE]

Over the event $B$, $|\Delta_{i,J}|\leq Cn^{-1/2}\log(n\vee p)$. As a result,

[TABLE]

Plugging it into (76), we obtain

[TABLE]

Next, we bound $(K_{1})$ and $(K_{2})$. Note that $\bm{S}_{2,J}-\mathbb{E}\bm{S}_{2,J}=\frac{1}{n}\sum_{k=1}^{n}[\bm{X}_{k,J}\bm{X}_{k,J}-\mathbb{E}(\bm{X}_{k,J}\bm{X}_{k,J})]$. This allows us to re-write

[TABLE]

It is not hard to see that $\mathbb{E}|Q_{ii}|\leq C|{\cal A}|$ and that $\mathbb{E}|Q_{ik}Q_{i^{\prime}k^{\prime}}|\leq C|{\cal A}|^{2}$ when $\{i,k,i^{\prime},k^{\prime}\}$ has at least two distinct values (see the proof of Theorem 2.1 for similar arguments). As a result,

[TABLE]

Moreover, noting that $\mathbb{E}Q_{ik}=0$ for $i\neq k$, we have $\mathbb{E}(Q_{ik}Q_{i^{\prime}k^{\prime}})=0$ for $\{i,k,i^{\prime},k^{\prime}\}$ that are mutually distinct. It follows that

[TABLE]

Combining the above gives

[TABLE]

Similarly, since $\bm{S}_{1,J}=\frac{1}{n}\sum_{k=1}^{n}\bm{X}_{i,J}$, we re-write

[TABLE]

Then, $\mathbb{E}R_{ik}=0$ for $i\neq k$, $\mathbb{E}|R_{ii}|\leq C|{\cal A}|$, and $\mathbb{E}(R_{ik}R_{i^{\prime}k^{\prime}})\leq C|{\cal A}|^{2}$ when $\{i,k,i^{\prime},k^{\prime}\}$ has at least two distinct values. As a result,

[TABLE]

We immediately have

[TABLE]

Plugging (86)-(87) into (82) gives (75). The claim then follows.

A.7 Proof of Theorem 3.2

Similar to the proof of Theorem 2.1, let $\widetilde{\theta}^{\,\mathrm{B}}_{m}$ and $\widehat{\theta}^{\,\mathrm{B}}_{m}$ denote the BAE with true $(\boldsymbol{\mu},\boldsymbol{\Sigma})$ and estimates $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$; here, $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$ may not be the sample mean and sample covariance matrix. By Theorem 3.3, $\mathbb{E}[(\widetilde{\theta}^{\,\mathrm{M}}_{m}-\theta_{m})^{2}]\leq Cn^{-1}$. It follows from the Markov’s inequality that, for any $\epsilon>0$, there is a constant $C_{\epsilon}>0$ such that, with probability $1-\epsilon/2$,

[TABLE]

To show the claim, it suffices to show that, there is a constant $C^{\prime}_{\epsilon}>0$ such that with probability $1-\epsilon/2$,

[TABLE]

We now show (88). Let $\bm{X}_{i,J}=\boldsymbol{\Sigma}_{J,J}^{-1/2}(\bm{Y}_{i,J}-\boldsymbol{\mu}_{J})$ and $\widehat{\bm{X}}_{i,J}=\widehat{\boldsymbol{\Sigma}}_{J,J}^{-1/2}(\bm{Y}_{i,J}-\widehat{\boldsymbol{\mu}}_{J})$. Then,

[TABLE]

By direct calculations,

[TABLE]

As a result,

[TABLE]

Introduce

[TABLE]

Then, (94) can be rewritten as

[TABLE]

First, we study the main terms in (99). Note that $\widetilde{\bm{S}}^{(m)}_{2,J}$ is the sample covariance matrix of $\{\|\bm{X}_{i,J}\|^{m-1}\bm{X}_{i,J}:1\leq i\leq n\}$, and $\widetilde{\bm{S}}_{1,J}^{(m)}$ is the sample mean of $\{\|\bm{X}_{i,J}\|^{2m-2}\bm{X}_{i,J}:1\leq i\leq n\}$. Using similar calculations as in the proof of Theorem 3.3, we can prove that

[TABLE]

Combining it with the Markov inequality, for any $\epsilon>0$, there is $C>0$ such that, with probability $1-\epsilon/4$, \bigl{\|}\frac{1}{|{\cal A}|}\sum_{J\in{\cal A}}\widetilde{\bm{S}}^{(m)}_{2,J}\bigr{\|}\leq C and \bigl{\|}\frac{1}{|{\cal A}|}\sum_{J\in{\cal A}}\widetilde{\bm{S}}^{(m)}_{1,J}\bigr{\|}\leq C. On this event, the sum of the first two terms in (99) is bounded in absolute value by

[TABLE]

Next, we study the remainder terms in (99). By (89) and our assumption on $(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}})$, we have

[TABLE]

It follows that $\|\Delta_{i,J}\|^{k}\leq C\beta^{k}_{n}\|\bm{X}_{i,J}\|^{2k}+C\alpha^{k}_{n}\|\bm{X}_{i,J}\|^{k}$. Then,

[TABLE]

Using similar calculations as in the proof of Theorem 3.3, we can prove that $\frac{1}{n|{\cal A}|}\sum_{J\in{\cal A}}\sum_{i=1}^{n}\mathbb{E}\|\bm{X}_{i,J}\|^{k}\leq C$, for all $1\leq k\leq 4m$. It follows from the Markov inequality that, for a constant $C>0$, with probability $1-\epsilon/4$, $\frac{1}{n|{\cal A}|}\sum_{J\in{\cal A}}\sum_{i=1}^{n}\|\bm{X}_{i,J}\|^{k}\leq C$, for all $1\leq k\leq 2m$. On this event,

[TABLE]

Combining (101) and (103) gives $|\widehat{\theta}^{\,\mathrm{B}}_{m}-\widetilde{\theta}^{\,\mathrm{B}}_{m}|\leq C\max\{\alpha_{n},\beta_{n}\}$. This proves (88), and the claim follows immediately.

A.8 Proof of Theorem 3.3

Fix a collection ${\cal A}$ of blocks. Write for short $\widehat{\theta}_{m}^{B}=\widehat{\theta}_{m}^{B}(\boldsymbol{\mu},\mathrm{diag}_{\cal A}(\boldsymbol{\Sigma}))$. For preparation, first, we verify that $\widehat{\theta}_{m}^{B}$ is an unbiased estimator. For any $J\in{\cal A}$, by (12) and the fact that $\|\bm{U}_{|J|}\|=1$, we have

[TABLE]

As a result,

[TABLE]

In particular, it implies that

[TABLE]

Therefore, $\widehat{\theta}_{m}^{B}$ is unbiased. Additionally, we have

[TABLE]

Second, we introduce an alternative expression of $\theta_{m}\cdot c^{*}_{m,|J|}$. Consider the special case $\xi\sim\chi_{p}^{2}$. Since $\bm{Y}\sim N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ in this case, we then have $\bm{Y}_{J}\sim N(\boldsymbol{\mu}_{J},\boldsymbol{\Sigma}_{J,J})$ and $(\bm{Y}_{J}-\boldsymbol{\mu}_{J})^{{\rm T}}\boldsymbol{\Sigma}_{J,J}^{-1}(\bm{Y}_{J}-\boldsymbol{\mu}_{J})\sim\chi^{2}_{|J|}$. Hence, in (104), the left hand side equals to $\mathbb{E}\chi^{2m}_{|J|}$. At the same time, the right hand side is equal to $\theta_{m}\cdot c_{m,|J|}^{*}=p^{-m}\mathbb{E}\xi^{2m}\cdot c_{m,|J|}^{*}=p^{-m}\mathbb{E}\chi_{p}^{2m}\cdot c_{m,|J|}^{*}$. Equating the left/right hand sides gives

[TABLE]

We combine it with the definition of $\theta_{m}=p^{-m}\mathbb{E}\xi^{2m}$ and $r_{m}=\mathbb{E}\xi^{2m}/\mathbb{E}\chi_{p}^{2m}$. It implies that

[TABLE]

We now show the claim. For $J\in{\cal A}$, let $W_{J}=[(\bm{Y}_{J}-\boldsymbol{\mu}_{J})^{{\rm T}}\boldsymbol{\Sigma}_{J,J}^{-1}(\bm{Y}_{J}-\boldsymbol{\mu}_{J})]^{m}$. By (105)-(106),

[TABLE]

Consider $(I)$. Combining (104) and (106), we have

[TABLE]

Hence,

[TABLE]

where the last two lines are from Definition 3.1.

Consider $(II)$. Fix $I$ and $J$. Note that

[TABLE]

We have had an expression of $\mathbb{E}W_{I}$ as in (110). We still need to get an expression of $\mathbb{E}(W_{I}W_{J})$. For the set $I\cup J$, we apply (12) and find that

[TABLE]

where $B$ is a Beta distribution with parameters $\frac{|I|+|J|}{2}$ and $\frac{p-(|I|+|J|)}{2}$. Let $\widetilde{\bm{U}}_{I}$ and $\widetilde{\bm{U}}_{J}$ be the vectors formed by the first $|I|$ coordinates and the last $|J|$ coordinates of $\boldsymbol{\Sigma}^{1/2}_{I\cup J,I\cup J}\bm{U}_{|I|+|J|}$, respectively. We then have $W_{I}=\xi^{2m}B^{m}\|\boldsymbol{\Sigma}_{II}^{-1/2}\widetilde{\bm{U}}_{I}\|^{2m}$ and $W_{J}=\xi^{2m}B^{m}\|\boldsymbol{\Sigma}_{JJ}^{-1/2}\widetilde{\bm{U}}_{J}\|^{2m}$. As a result,

[TABLE]

We then use the cross-moments of multivariate normal distributions to get the last term above. Let $\xi_{0}^{2}\sim\chi_{p}^{2}$ be a random variable independent of $B$ and $\bm{U}_{|I|+|J|}$. The random vector

[TABLE]

It follows that

[TABLE]

Write $\widetilde{\bm{Z}}_{1}=\boldsymbol{\Sigma}_{II}^{-1/2}\bm{Z}_{I}$ and $\widetilde{\bm{Z}}_{2}=\boldsymbol{\Sigma}_{JJ}^{-1/2}\bm{Z}_{J}$. Note that

[TABLE]

Combining (115) and (116) gives

[TABLE]

We now combine (110) and (118) and note that $\|\widetilde{\bm{Z}}_{1}\|^{2}\sim\chi_{|I|}^{2}$ and $\|\widetilde{\bm{Z}}_{2}\|^{2}\sim\chi_{|J|}^{2}$. It yields

[TABLE]

As a result,

[TABLE]

We now plug (111) and (119) into (107). It gives

[TABLE]

What remains is to bound the last term. Since the random vectors $\widetilde{\bm{Z}}_{1}$ and $\widetilde{\bm{Z}}_{2}$ jointly follow a multivariate normal distribution as dictated in (117), we can apply the following lemma:

Lemma A.2**.**

Let $\bm{Z}_{1}$ and $\bm{Z}_{2}$ be two random vectors such that

[TABLE]

Then, for a constant $\widetilde{C}_{m}>0$ that only depends on $m$ but is independent of $(k_{1},k_{2})$,

[TABLE]

We combine Lemma A.2 with (117) and them plug it into (121). It follows that

[TABLE]

This proves the claim.

Appendix B Supplementary proofs

B.1 Proof of Lemma A.1

Let \theta=\arcsin\bigl{(}\mathop{\rm sign}(\rho)\cdot\sqrt{|\rho|}\bigr{)}\in[-\frac{\pi}{2},\frac{\pi}{2}]. We then have $\sin\theta=\mathop{\rm sign}(\rho)\cdot\sqrt{|\rho|}$ and $\cos\theta=\sqrt{1-|\rho|}$. Let $U_{1},U_{2},V$ be iid $N(0,1)$ random variables. It is easy to see that

[TABLE]

For notation simplicity, we omit the superscript $(d)$ in all equations. It follows that

[TABLE]

Then,

[TABLE]

Note that for random variables $(X,Y,W_{1},W_{2})$, when $X$, $Y$ and $(W_{1},W_{2})$ are mutually independent, $\mathop{\rm cov}(XW_{1},YW_{2})=\mathbb{E}X\cdot\mathbb{E}Y\cdot\mathop{\rm cov}(W_{1},W_{2})$. Plugging it into the above expression, we obtain

[TABLE]

Using our previous notations, $\eta_{m}$ is the $2m$-th moment of $N(0,1)$. By elementary statistics, $\eta_{m}=(2m-1)!!=\prod_{j=0}^{m-1}(1+2j)$. Using this formula, we can prove $\mathbb{E}V^{2s}-\mathbb{E}V^{2k_{1}}\mathbb{E}V^{2k_{2}}\geq 0$. Hence,

[TABLE]

At the same time, we note that $\cos^{2}\theta=1-|\rho|$ and $\sin^{2}\theta=|\rho|$. It follows that

[TABLE]

The claim then follows.

B.2 Proof of Lemma A.2

Suppose the rank of $\boldsymbol{\Gamma}$ is $k\leq\min\{k_{1},k_{2}\}$. Let $\boldsymbol{\Gamma}={\bf H}_{1}\boldsymbol{\Lambda}{\bf H}_{2}^{{\rm T}}$ be the singular value decomposition of $\boldsymbol{\Gamma}$. We note that all singular values have an absolute value no larger than $1$. For $\ell=1,2$, let $\widetilde{{\bf H}}_{\ell}\in\mathbb{R}^{k_{\ell},k_{\ell}-k}$ be such that $[{\bf H}_{\ell},\widetilde{{\bf H}}_{\ell}]$ form an orthogonal basis of $\mathbb{R}^{k_{\ell}}$. Define

[TABLE]

It is easy to see that ${\bf A}_{\ell}{\bf A}^{\prime}_{\ell}={\bf I}-{\bf H}_{\ell}\boldsymbol{\Lambda}{\bf H}_{\ell}^{\prime}$. Let $\bm{X}_{1}\sim N({\bf 0},{\bf I}_{k_{1}})$, $\bm{X}_{2}\sim N({\bf 0},{\bf I}_{k_{2}})$, and $\bm{W}\sim N({\bf 0},{\bf I}_{k})$ be mutually independent random variables. We claim that

[TABLE]

This can be verified by computing the covariance matrix of the right hand side. We shall omit the superscript $(d)$ in all equations for notation simplicity. Write $\bm{X}_{\ell}=(\bm{X}_{\ell 1}^{{\rm T}},\bm{X}_{\ell 2}^{{\rm T}})^{{\rm T}}$, corresponding to the first $k_{\ell}$ and the last $(k_{\ell}-k)$ coordinates, respectively, $\ell=1,2$. It follows that

[TABLE]

where the third line is from the zero mean and mutual independence of $(\bm{X}_{\ell 1},\bm{X}_{\ell 2},\bm{W})$ and the last line is due to that ${\bf H}_{\ell}^{\prime}{\bf H}_{\ell}={\bf I}_{k}$ and $\widetilde{{\bf H}}_{\ell}\widetilde{{\bf H}}_{\ell}^{\prime}={\bf I}_{k_{\ell}-k}$. Since $(U_{1},U_{2},V)$ are mutually independent, it follows that

[TABLE]

It is not hard to see that $\mathop{\rm cov}(V^{j_{1}},V^{j_{2}})\geq 0$. Hence, \mathop{\rm cov}\bigl{(}\|\bm{Z}_{1}\|^{2m},\|\bm{Z}_{2}\|^{2m}\bigr{)}\geq 0. Furthermore, since all entries of the diagonal matrix $\boldsymbol{\Lambda}$ are between [math] and $1$, we have

[TABLE]

where $X_{\ell}(j)$’s and $W(j)$’s are all $iid$ standard normal variables. In particular,

[TABLE]

Plugging these results into (127) gives

[TABLE]

We note that $m$ is bounded, but $(k_{1},k_{2},k)$ can grow with $(n,p)$. Note that $\mathbb{E}\chi_{k}^{2m}=\prod_{j=0}^{m-1}(k+2j)$ for all $k,m\geq 1$. As a result,

[TABLE]

Therefore,

[TABLE]

Noticing that $\boldsymbol{\Lambda}$ is a diagonal matrix containing the singular values of $\boldsymbol{\Gamma}$, we have proved the claim.

Appendix C The case of multivariate Gaussian distributions

We present a corollary about the errors of MAE and BAE for the special case of multivariate Gaussian distributions. Here $R_{n}(\widehat{\theta}_{2})=\mathbb{E}[(\widehat{\theta}_{2}-\theta_{2})^{2}/\theta_{2}^{2}]$. The proof is elementary and omitted.

Corollary C.1**.**

Let $\bm{Y}_{1},\cdots,\bm{Y}_{n}$ be i.i.d. samples of $N(\boldsymbol{\mu},\boldsymbol{\Sigma})$. For a constant integer $k\geq 2$, we assume the blocks in BAE are $J_{i}=\{(i-1)k+1,(i-1)k+2,\cdots,\min\{ik,p\}\}$, $1\leq i\leq\lceil p/k\rceil$.

• •

Suppose $\boldsymbol{\Sigma}={\bf I}_{p}$. Then, $R_{n}(\widehat{\theta}_{2}^{\,\mathrm{I}})\sim\frac{8}{np}$, $R_{n}(\widehat{\theta}_{2}^{\,\mathrm{M}})\sim\frac{32}{3np}$, and $R_{n}(\widehat{\theta}^{\,\mathrm{B}}_{2})\sim\frac{8(k+3)}{(k+2)np}$.

• •

Suppose $\boldsymbol{\Sigma}$ is a block-wise diagonal matrix with $2\times 2$ blocks, where each block has diagonals $1$ and off-diagonals $\rho\in(-1,1)$. Let $k=2$ in BAE. Then, $R_{n}(\widehat{\theta}_{2}^{\,\mathrm{I}})\sim\frac{8}{np}$, $R_{n}(\widehat{\theta}_{2}^{\,\mathrm{M}})\sim\frac{8(4+3\rho^{2}+\rho^{4})}{3np}$, and $R_{n}(\widehat{\theta}^{\,\mathrm{B}}_{2})\sim\frac{10}{np}$.

Appendix D Simulations for the estimator in Section 5

We conducted simulations to investigate the performance of the estimator of realized $\xi_{t}$ in Section 5.

In the first experiment, we generate $\{\bm{Y}_{t}\}_{t=1}^{T}$ iid from model (1) with a constant covariance matrix $\boldsymbol{\Sigma}$. The covariance is set to be $\Sigma_{ij}=0.3^{|i-j|}$, which is approximately banded. We fix $T=100$ and let $p$ varies. The results are displayed in Figure 9, where we study both cases of multivariate Gaussian data and multivariate $t_{4.5}$ data. We see that the estimated values are very close to the true values in all the cases.

In the second experiment, we generate data using the calibrated covariance matrix from S&P500 stock returns as in Section 4. In this case, the covariance matrix is heavily non-sparse, however, our estimator still works very well, no matter for Gaussian data or heavy-tailed data with multivariate $t$-distributions.