# Mean Test with Fewer Observation than Dimension and Ratio Unbiased Estimator
for Correlation Matrix

Tiefeng Jiang
School of Statistics
University of Minnesota
224 Church St SE Minneapolis, MN 55455
<EMAIL_ADDRESS>Ping Li
Cognitive Computing Lab
Baidu Research
10900 NE 8th St. Bellevue, WA 98004
<EMAIL_ADDRESS>
###### Abstract

Hotelling’s T-squared test is a classical tool to test if the normal mean of a
multivariate normal distribution is a specified one or the means of two
multivariate normal means are equal. When the population dimension is higher
than the sample size, the test is no longer applicable. Under this situation,
in this paper we revisit the tests proposed by Srivastava and Du (2008), who
revise the Hotelling’s statistics by replacing Wishart matrices with their
diagonal matrices. They show the revised statistics are asymptotically normal.
We use the random matrix theory to examine their statistics again and find
that their discovery is just part of the big picture. In fact, we prove that
their statistics, decided by the Euclidean norm of the population correlation
matrix, can go to normal, mixing chi-squared distributions and a convolution
of both. Examples are provided to show the phase transition phenomenon between
the normal and mixing chi-squared distributions. The second contribution of
ours is a rigorous derivation of an asymptotic ratio-unbiased-estimator of the
squared Euclidean norm of the correlation matrix.

##  1 Introduction

Among many statistical hypothesis testing problems, the Hotelling’s $T^{2}$
tests (Hotelling, 1931) are classic ones to study if the mean of a
multivariate normal distribution is equal to the given one, or if the means of
two normal distributions are equal. Let us quickly review the two problems.
First, assume $\\{{\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}\\}$ is a
random sample from a $p$-dimensional normal distribution
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$, where ${\boldsymbol{\mu}}$
is the mean vector and ${\boldsymbol{\Sigma}}$ is the covariance matrix.
Consider the test
$H_{0}:{\boldsymbol{\mu}}=\boldsymbol{0}~{}~{}\textrm{versus}~{}~{}H_{1}:{\boldsymbol{\mu}}\not=\boldsymbol{0}.$
The Hotelling $T^{2}$ test statistic is defined by
$T^{2}=(\bar{{\boldsymbol{X}}}-\mu)^{T}\hat{{\bf
S}}^{-1}(\bar{{\boldsymbol{X}}}-\mu)$, where $\bar{{\boldsymbol{X}}}$ and
$\hat{{\bf S}}$ are the sample mean and the sample covariance matrix. For the
two population test, assume that
$\\{{\boldsymbol{X}}_{i1},\cdots,{\boldsymbol{X}}_{in_{i}}\\}$ for $i=1,2$ are
two independent random samples from
$N_{p}({\boldsymbol{\mu}}_{1},{\boldsymbol{\Sigma}})$ and
$N_{p}({\boldsymbol{\mu}}_{2},{\boldsymbol{\Sigma}})$, respectively. We aim to
test $H_{0}:{\boldsymbol{\mu}}_{1}={\boldsymbol{\mu}}_{2}\ \ \mbox{versus}\ \
H_{1}:{\boldsymbol{\mu}}_{1}\neq{\boldsymbol{\mu}}_{2}.$ The corresponding
Hotelling’s $T^{2}$ test statistic is given by
$\frac{n_{1}n_{2}}{n_{1}+n_{2}}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})^{T}\hat{{\bf
S}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})$ where
$\bar{{\boldsymbol{X}}}_{i}$ is the sample mean vector of the $i$-th sample
and $\hat{{\bf S}}$ is the pooled sample covariance matrix.

As a consequence of the likelihood ratio test, the above two Hotelling $T^{2}$
tests have very nice properties. For example, for the one-sample test, $T^{2}$
is invariant under linear transformations. Also, $T^{2}$ is a uniformly most
powerful test under the general linear group; see, for example, p. 211 from
Muirhead (1982) or p. 190 from Anderson (2003).

Despite its important role in classical statistics, Hotelling’s $T^{2}$ have
some limitations. First, in order to guarantee that the sample covariance
matrix is invertible, the sample size have to be larger than the population
dimension. These are not true for some modern data in which the population
dimension is larger or even much larger than the sample size. For example, for
DNA microarray data, thousands of gene expression levels are often measured on
a small number of subjects. From http://genomics-pubs.princeton.edu/oncology/,
which is a popular colon dataset, one can see $p=2000$, $n_{1}=22$ and
$n_{2}=40$. Second, as Bai and Saranadasa (1996) have showed, the Hotelling
two-sample test is inconsistent as the population dimension and the sample
size are comparable.

To accommodate data with the feature of large $p$ small $n$ for the above one-
sample and two-sample testing problems, Dempster (1958, 1960) study the so-
called “non-exact” test. Bai and Saranadasa (1996) construct a new test
statistic by removing the inverse matrix from the definition of $T^{2}$ given
earlier. Starting from this century, investigators begin to conceive new test
statistics. Srivastava and Du (2008), Srivastava (2009), Chen and Qin (2010),
and Dong et al. (2016) replace the inverse of the sample covariance matrix by
mild quantities. Lopes et al. (2011) and Srivastava et al. (2016) use a sort
of dimension reduction method to lower the population dimension and then use
Hotelling’s $T^{2}$ tests. For data with a certain of sparse nature, Cai et
al. (2014) establish tests by measuring the maximum componentwise mean
difference of appropriately transformed observations. Other research related
to sparsity can be seen from Zhong et al. (2013), Chen et al. (2014), Wang et
al. (2015), Gregory et al. (2015) and Guo and Chen (2016). The researchers
Biswas and Ghosh (2014), Chang et al. (2017), Chakraborty and Chaudhuri
(2017), and Xue and Yao (2020) consider the above test under non-normal
assumptions. Other contributions include Park and Ayyala (2013) and Wu et al.
(2006) for the consideration of scale-invariant tests and Gretton et al.
(2012) for a kernel-based discrepancy measure. Feng et al. (2016) work on a
two-sample location problem via a multivariate-sign-based high-dimensional
test. Xu et al. (2016) and Zhang et al. (2020a) construct test statistics for
the two-sample test problem by $l_{q}$-norm. Recently, built upon Chen et al.
(2011), the authors (Li et al., 2020) investigate the Hotelling’s $T^{2}$ by a
ridge-regularized method. Zhang et al. (2020b) extend the tests by Wu et al.
(2006) and Zhang et al. (2020a) through a modification of the test by
Srivastava and Du (2008).

In this paper, we will revisit the procedure by Srivastava and Du (2008) for
the one-sample and two-sample problems of testing means. We will briefly
review them next and state our motivation.

Let $\\{{\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}\\}$ be a random
sample from a $p$-dimensional normal distribution
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with correlation matrix
$\mathbb{R}$. The sample mean and the sample covariance matrix defined by

$\displaystyle\bar{{\boldsymbol{X}}}=\frac{1}{n}\sum_{i=1}^{n}{\boldsymbol{X}}_{i}\
\ \mbox{and}\ \ \hat{{\bf
S}}=\frac{1}{n}\sum_{i=1}^{n}({\boldsymbol{X}}_{i}-\bar{{\boldsymbol{X}}})({\boldsymbol{X}}_{i}-\bar{{\boldsymbol{X}}})^{T}.$
(1)

Let $\hat{{\bf D}}$ be the diagonal matrix of $\hat{{\bf S}}$ and $\hat{{\bf
R}}$ be the sample correlation matrix defined by $\hat{{\bf R}}=\hat{{\bf
D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2}.$ For the testing problem
$H_{0}:{\boldsymbol{\mu}}=\boldsymbol{0}$ vs
$H_{1}:{\boldsymbol{\mu}}\not=\boldsymbol{0}$, Srivastava and Du (2008)
propose the following test statistic.

$\displaystyle T_{SD}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-p(n-1)(n-3)^{-1}}{\sqrt{2[\mathrm{tr}(\hat{{\bf
R}}^{2})-p^{2}(n-1)^{-1}]}}.$ (2)

Under certain conditions, they show $T_{SD}$ converges to a normal
distribution.

For two population case, the test is
$H_{0}:{\boldsymbol{\mu}}_{1}={\boldsymbol{\mu}}_{2}$ vs
$H_{1}:{\boldsymbol{\mu}}_{1}\neq{\boldsymbol{\mu}}_{2}$. Assume that
$\\{{\boldsymbol{X}}_{i1},\cdots,{\boldsymbol{X}}_{in_{i}}\\}$ for $i=1,2$ are
two independent random samples from
$N_{p}({\boldsymbol{\mu}}_{1},{\boldsymbol{\Sigma}})$ and
$N_{p}({\boldsymbol{\mu}}_{2},{\boldsymbol{\Sigma}})$, respectively. Let
$\bar{{\boldsymbol{X}}}_{i}$ be the sample mean for the $i$-th sample and
$\hat{{\bf S}}$ be the pooled sample covariance matrix defined by

$\displaystyle\hat{{\bf
S}}=\frac{1}{n_{1}+n_{2}}\Big{[}\sum_{j=1}^{n_{1}}(\mathbb{X}_{1j}-\bar{\mathbb{X}}_{1})(\mathbb{X}_{1j}-\bar{\mathbb{X}}_{1})^{T}+\sum_{j=1}^{n_{2}}(\mathbb{X}_{2j}-\bar{\mathbb{X}}_{2})(\mathbb{X}_{2j}-\bar{\mathbb{X}}_{2})^{T}\Big{]}.$
(3)

Assume $\hat{{\bf D}}$ is the diagonal matrix of $\hat{{\bf S}}$ and
$\hat{{\bf R}}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2}$ is the
pooled sample correlation matrix. Srivastava and Du (2008) consider the
following statistic defined by

$\displaystyle
T^{\prime}_{SD}=\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-\frac{(n_{1}+n_{2}-2)p}{n_{1}+n_{2}-4}}{\sqrt{2\big{[}\mathrm{tr}(\hat{{\bf
R}}^{2})-\frac{p^{2}}{n_{1}+n_{2}-2}\big{]}}}.$ (4)

Under certain conditions, they show $T^{\prime}_{SD}$ converges to a normal
distribution. In fact, to improve the convergence speeds, they actually add a
term $c_{p,n}=1+\frac{\mathrm{tr}(\hat{{\bf R}}^{2})}{p^{3/2}}$ under the
squared roots in the denominators of $T_{SD}$ and $T_{SD}^{\prime}$,
respectively. Based on their assumptions, $c_{p,n}$ goes to one. However, we
will study a more general case in which $c_{p,n}$ may not go to one, and
sometimes it even goes to infinity. This is the reason we dump the term
$c_{p,n}$ from both $T_{SD}$ and $T_{SD}^{\prime}$, respectively.

Evidently, $T_{SD}$ and $T^{\prime}_{SD}$ have a nice property of scale-
invariance, i.e., they are not changed if data are multiplied by a constant.
Also, they can be directly computed. Our motivation in this paper to reexamine
the tests by Srivastava and Du (2008) has three folds. (a) The conditions to
guarantee the central limit theorems of $T_{SD}$ and $T^{\prime}_{SD}$ are
somehow stringent; see (8). Also the sample size and population dimension have
to satisfy that $p=o(n^{2})$. To make the method more applicable, we hope to
relax the condition imposed on ${\boldsymbol{\Sigma}}$ as well as that $p$ on
$n$. (b) We would like to understand an interesting observation by Zhang et
al. (2020b) from their simulation: the distribution of $T^{\prime}_{SD}$
sometimes looks like a normal curve, and other times it looks like a chi-
square curve. (c) We plan to give a rigorous proof of the major ingredient of
this theory, that is, an asymptotic ratio-unbiased-estimator of
$\mbox{tr}(\mathbb{R}^{2})$, where $\mathbb{R}$ is the population correlation
matrix of $N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ aforementioned.

Now we state our findings. For (a), we have obtained the asymptotic
distributions of $T_{SD}$ and $T^{\prime}_{SD}$ for arbitrary
${\boldsymbol{\Sigma}}$ in Theorems 1 and 2, respectively. In particular,
Theorem 1 holds for two extreme cases: the independent case with
$\mathbb{R}=\mathbb{I}_{p}$ and the most dependent case, i.e., all entries of
the population vector are identical. Our theory says the restriction on $p$
and $n$ will be $p=o(n^{2})$ if the entries of the population vector are not
far from independent. However, for dependent or very dependent case, our
conclusion holds as long as $p$ is not more than a polynomial order of $n$.
For (b), we successfully understand the observation by Zhang et al. (2020b).
In fact, there are indeed transition phenomena of the limiting distributions
of $T_{SD}$ and $T^{\prime}_{SD}$. They are sometimes normal, mixing chi-
squared or the sum of two independent random variables, one has normal
distribution and the other has a mixing chi-squared distribution. This can be
quickly seen from (11) and (17). As for (c), for arbitrary
${\boldsymbol{\Sigma}}$ we have proved rigorously that an asymptotic ratio-
unbiased-estimator of $\mbox{tr}(\mathbb{R}^{2})$ is
$\mbox{tr}(\hat{\mathbb{R}}^{2})-\frac{p(p-1)}{n-1}$ in Theorem 3. Especially
the theorem is true for two extreme cases: ${\boldsymbol{\Sigma}}$ is diagonal
or proportional to a matrix whose entries are all equal to $1$. In addition,
some of our partial calculations and heuristics indicate that the above phase
transition from the Gaussian to mixing chi-squared distributions is possibly a
universal phenomenon. We will present three testing procedures to justify this
claim.

The above solution is conducted through the understanding of the sample
correlation matrix $\hat{{\bf R}}$, a special random matrix, defined below
(1). Unlike the Gaussian orhtogonal/unitary/symplectic ensemble or Wishart
matrices, $\hat{{\bf R}}$ lacks the orthogonal-invariant property. As a
consequence, investigating the sample correlation matrix always cost extra
energy than working on other popular matrices. We employ the machinery for
$\hat{{\bf R}}$ developed by Jiang (2004), Jiang (2009), Cai et al. (2013),
Jiang and Yang (2013), and Fan and Jiang (2019). In particular, we extend the
method conceived by Jiang (2004) and Fan and Jiang (2019) to prove a weak law
of large numbers for $\mbox{tr}(\hat{\mathbb{R}}^{2})$ en route to obtain an
asymptotic ratio-unbiased-estimator. More elaboration are provided in Section
4.

The remaining of the paper is organized as follows. The one-sample mean and
two-sample mean problems are studied in Sections 2 and 3, respectively. In
Section 4, we give an asymptotic ratio-unbiased-estimator of
$\mbox{tr}(\mathbb{R}^{2})$, which is more applicable than the one by
Srivastava and Du (2008) and Srivastava (2009) (their result lacks a rigorous
proof although it has been used in literature). The concluding remarks and
discussions are presented in Section 5, in which we particularly point out our
findings on phase transitions between normal and mixing chi-squared
distributions may also exist for some other testing procedures. Finally, the
proofs are given in Section 6.

##  2 One Sample Mean Test for Large $p$ and Small $n$

Let $\\{{\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}\\}$ be a random
sample from $N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$. Consider the
test that

$\displaystyle
H_{0}:{\boldsymbol{\mu}}=\boldsymbol{0}~{}~{}\textrm{versus}~{}~{}H_{1}:{\boldsymbol{\mu}}\not=\boldsymbol{0}.$
(5)

Let $\mathbb{D}$ be the diagonal matrix of $\mathbb{\Sigma}.$ Then the
$p\times p$ population correlation matrix is

$\displaystyle\mathbb{R}=\mathbb{D}^{-1/2}\mathbb{\Sigma}\mathbb{D}^{-1/2}.$
(6)

Similarly, let $\hat{{\bf D}}$ be the diagonal matrix of $\hat{{\bf S}}$ from
(1). Then the sample correlation matrix $\hat{{\bf R}}$ is defined by

$\displaystyle\hat{{\bf R}}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf
D}}^{-1/2}.$ (7)

Srivastava and Du (2008) and Srivastava (2009) obtain a result on $T_{SD}$
from (2) as follows.

###### Result 2.1

Assume $n=O(p^{\zeta})$, $\frac{1}{2}<\zeta\leq 1$ and

$\displaystyle
0<\lim_{p\to\infty}\frac{\mbox{tr}(\mathbb{R}^{i})}{p}<\infty,~{}~{}~{}~{}i=1,2,3,4.$
(8)

If ${\boldsymbol{\mu}}=\boldsymbol{0}$ then $T_{SD}\to N(0,1)$ in distribution
as $p\to\infty.$

A quick comment is that (8) holds automatically for $i=1$ since all of the
diagonal entries of $\mathbb{R}$ are equal to $1$, and hence
$\mbox{tr}(\mathbb{R})=p$. In order to make approximation better, we now
revise the statistic $T_{SD}$ slightly. Set

$\displaystyle T_{p,1}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-pn(n-3)^{-1}}{\sqrt{2\big{|}\mathrm{tr}(\hat{{\bf
R}}^{2})-p(p-1)(n-1)^{-1}\big{|}}}.$ (9)

We are doing so because, by Lemma 6.3, the major contribution of the mean of
$\mbox{tr}(\hat{\mathbb{R}}^{2})$ is
$\mbox{tr}(\mathbb{R}^{2})+p(p-1)(n-1)^{-1}$. Also, by Theorem 3,
$\mathrm{tr}(\hat{{\bf R}}^{2})-p(p-1)(n-1)^{-1}>0$ as $n$ and $p$ are large
enough. On the other hand, with probability one, $\mathrm{tr}(\hat{{\bf
R}}^{2})-p(p-1)(n-1)^{-1}\neq 0$ because $\mathrm{tr}(\hat{{\bf R}}^{2})$ is a
continuous function of Gaussian random variables. A discussion between
$T_{p,1}$ and $T_{SD}$ will be elaborated shortly. Review the Frobenius norm,
sometimes also called the Euclidean norm,
$\|\mathbb{A}\|_{F}:=[\mbox{tr}(\mathbb{A}^{T}\mathbb{A})]^{1/2}$ for any
matrix $\mathbb{A}$. For mathematical rigor, we assume the sample sizes $n$
depends on $p$. According to Hu and Bai (2016), the behavior of the statistic
$T_{SD}$ is not known if ${\boldsymbol{\Sigma}}$ has spikes. The following
gives a complete characterization of the properties of $T_{SD}$ and $T_{p,1}$
in terms of the spikes of the correlation matrix $\mathbb{R}.$

######  THEOREM 1

Suppose ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ is a random sample
from $N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$. Let $\mathbb{R}$ be
the correlation matrix as in (6) with eigenvalues
$\lambda_{1}\geq\cdots\geq\lambda_{p}\geq 0$. Assume

(a) $\lim_{p\to\infty}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}=\rho_{i}\geq 0$
for all $i\geq 1$;

(b) $\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ and
$\lim_{p\to\infty}\frac{p}{n^{a}}=0$ for some constant $a>0$.
If ${\boldsymbol{\mu}}=\boldsymbol{0}$, then $T_{SD}\to
b\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$ in
distribution, where $\xi_{0},\xi_{1},\xi_{2},\cdots$ are i.i.d. $N(0,1)$ and
$b=(1-\sum_{i=1}^{\infty}\rho_{i}^{2})^{1/2}$. The same conclusion also holds
for $T_{p,1}$.

Condition (a) considers the possibility that $\mathbb{R}$ may have spikes. If
$\rho_{1}=0$, then $\rho_{i}=0$ for every $i\geq 1$ due to the monotonicity of
$\rho_{i}$, and we say there are no spikes in $\mathbb{R}$. In this situation
we have $T_{p,1}\to N(0,1)$ by Theorem 1. Using essentially the Fatou lemma,
we have checked the given conditions actually imply
$\sum_{i=1}^{\infty}\rho_{i}^{2}\leq 1$; see the proof of Lemma 4.2. If
$\sum_{i=1}^{\infty}\rho_{i}^{2}=1$, then $T_{p,1}$ converges to a mixing chi-
square distribution. If $\sum_{i=1}^{\infty}\rho_{i}^{2}\in(0,1)$, then the
asymptotic distribution is a blend of a Gaussian distribution and a mixing
chi-squared distribution. Condition (b) characterizes the restriction between
the sample size $n$ and population dimension $p$.

Observe that (8) with $i=3$ implies that
$\lambda_{1}^{3}\leq\mbox{tr}(\mathbb{R}^{3})=O(p)$. Thus,
$\lambda_{1}=O(p^{1/3})$. Since all of the diagonal entries of $\mathbb{R}$
are identical to $1$ then $\mbox{tr}(\mathbb{R}^{2})\geq p$. This shows that
$\rho_{1}=0$ where $\rho_{1}$ is from Theorem 1. Also, (8) with $i=2$ implies
that $\|\mathbb{R}\|$ and $\sqrt{p}$ have the same order, thus
$\frac{p}{n\|\mathbb{R}\|_{F}}\to 0$ as $p\to\infty$. This together with the
condition “$n=O(p^{\zeta})$, $\frac{1}{2}<\zeta\leq 1$” implies (b) from
Theorem 1. So our theorem is more general than Result 2.1 by Srivastava and Du
(2008).

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
an $AR(1)$ model with $\mathbb{R}=(\gamma^{|i-j|})$ and the absolute values of
$\gamma=\gamma_{p}$ staying away from $1$. By using the Gersgorin disc theorem
[see, e.g., p. 344 from Horn and Johnson (2012)], the largest eigenvalue or
$\mathbb{R}$ is of order $O(1)$. Hence condition (a) of Theorem 1 holds with
$\rho_{1}=0$. If condition (b) also holds, then both $T_{SD}$ and $T_{p,1}$ go
to the standard normal distribution. The same conclusion is also valid for a
banded correlation matrix $\mathbb{R}=(r_{ij})$ with $r_{ij}=0$ for $|j-i|\geq
t$ where $t=t_{p}=o(\sqrt{p})$. In this case, the largest eigenvalue or
$\mathbb{R}$ is of order $o(\sqrt{p})$. Similar results can be obtained for
other patterned matrices including Toeplitz matrices, Hankel matrices and
symmetric circulant matrices; see, e.g., Brockwell et al. (2016).

Now let us look at some special features of Theorem 1. First, we do not need
the population matrix ${\boldsymbol{\Sigma}}$ to be invertible as required in
Hotelling’s $T^{2}$ test (Hotelling, 1931). The largest discrepancy between
$n$ and $p$ from Srivastava and Du (2008) is that $p=o(n^{2})$. Our range is
that $p$ can be at any polynomial order of $n$ provided
$p=o(n\|\mathbb{R}\|_{F})$. If the entries of the population vector are not
very far from independence, in the sense that $\|\mathbb{R}\|_{F}$ is in the
order of $\sqrt{p}$, then the restriction $p=o(n\|\mathbb{R}\|_{F})$ is
reduced to $p=o(n^{2})$. In the case that the entries of the population vector
are very dependent such that $\|\mathbb{R}\|_{F}$ is in the order of $p$, then
$p$ is allowed to take any polynomial order of $n$. To convince our readers
for the dependent case and to make further discussions, we next study two
extreme cases: independence and most dependence. The derivation of the results
below does not use any techniques and steps from the proof of Theorem 1.

######  PROPOSITION 1

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with
${\boldsymbol{\mu}}=\mathbb{0}$. Assume $n=n_{p}\to\infty$ as $p\to\infty$.
The following hold.

(i) Let all of the $p^{2}$ entries of ${\boldsymbol{\Sigma}}$ be identical.
Then both $T_{SD}$ and $T_{p,1}$ converge to
$\frac{1}{\sqrt{2}}[\chi^{2}(1)-1]$ in distribution as $p\to\infty$ regardless
of the relative speeds of $n$ and $p$.

(ii) Let ${\boldsymbol{\Sigma}}$ be a diagonal matrix whose diagonal entries
are all positive, equivalently, $\mathbb{R}=\mathbb{I}_{p}$. Then

$\displaystyle T_{SD}\to\begin{cases}\eta,&\text{if $p/n^{2}\to 0$};\\\
\eta+\sqrt{h/2},&\text{if $p/n^{2}\to h$};\\\ \infty,&\text{if
$p/n^{2}\to\infty$}\end{cases}$ (10)

in distribution, where $\eta\sim N(0,1)$. However, $T_{p,1}\to N(0,1)$ as
$p\to\infty$ regardless of the speeds of $n$ and $p$ going to infinity.

In case (i) we see Theorem 1 holds without the assumption that $p=o(n^{a})$
for some constant $a>0$. The conclusion for case (ii) says that the condition
$\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ in Theorem 1 is sharp, and
$T_{p,1}$ is better than $T_{SD}$. In particular, we do not need the
assumption $\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ to assure
$T_{p,1}\to N(0,1)$. This is simply a coincidence because of the special
structure of ${\boldsymbol{\Sigma}}$ or $\mathbb{R}$. The condition
$\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ is essentially required at
handling the denominator of $T_{p,1}$ for arbitrary $\mathbb{R}$; see (9) and
Theorem 3. Figure 1 clearly shows the existence of shifts between the density
curves of $T_{SD}$ and $N(0,1)$. However, with the ratio $p/n^{2}$ becoming
smaller, the shift diminishes gradually.

Figure 1: These pictures are based on Proposition 1, from $10^{5}$
simulations. It shows the centers of $T_{p,1}$ and $N(0,1)$ are always close,
and the existence of shifts between density curves of $T_{SD}$ and $N(0,1)$.
The shifts are smaller as $p/n^{2}$ become smaller.

Practically, given the population correlation $\mathbb{R}$, we need to
determine if it has spikes and how many spikes. Interested readers are
referred to Fan et al. (2020) and Morales-Jimenez et al. (2021) for procedures
to do so. We next present some examples. Their verification are presented in
Section 6.4.

###### Example 2.1

Given $r\in[0,1]$, define
$\mathbb{A}_{m}:=(1-r)\mathbb{I}_{m}+r\mathbb{J}_{m}$ for any $m\geq 1$, where
$\mathbb{I}_{p}$ and $\mathbb{J}_{p}$ denote the $p\times p$ identity matrix
and the $p\times p$ matrix of ones. Set $\mathbb{R}=\mathbb{A}_{p}$. Then
$\lambda_{1}=1+(p-1)r$ and $\lambda_{2}=\cdots=\lambda_{p}=1-r$. Assume
$r=r_{p}$ and $\lim_{p\to\infty}\sqrt{p}\cdot r=c\geq 0$. Then

$\displaystyle T_{p}\to\begin{cases}\xi_{0},&\text{if $c=0$};\\\
\frac{1}{\sqrt{c^{2}+1}}\xi_{0}+\frac{c}{\sqrt{2(c^{2}+1)}}(\xi_{1}^{2}-1),&\text{if
$c\in(0,\infty)$};\\\ \frac{1}{\sqrt{2}}\cdot(\xi_{1}^{2}-1),&\text{if
$c=\infty$}\end{cases}$ (11)

under condition $p=o(n^{2})$ for the case $c\in[0,\infty)$ and under the
conditions $nr\to\infty$ and $p=o(n^{a})$ for some constant $a>0$ for the case
$c=\infty$. Here $\xi_{0}$ and $\xi_{1}$ are i.i.d. $N(0,1)$ and $T_{p}$
stands for $T_{SD}$ or $T_{p,1}$. The phase transition occurs between the
normal and mixing chi-squared distributions as the entries of the correlation
matrix change their values.

We make Figures 2 and 3 based on Example 2.1 to compare the performances of
our $T_{p,1}$ and $T_{SD}$ by Srivastava and Du (2008). For different values
of $n,p,r$, the statistics $T_{p,1}$ and $T_{SD}$ are simulated for $10^{5}$
times. The elaboration is given below.

Figure 2 is designed for a small value of $r$ with $r=0.1$. Look at the first
two pictures on the top row, that is, the ones with $(n,p)=(10,50)$ and
$(n,p)=(20,100)$, respectively. Neither $T_{p,1}$ nor $T_{SD}$ perform well.
But this is expected because the sample sizes are very small. However, by
fixing the rate of $p/n$ and let $n$ and $p$ grow gradually, the
approximations become better and better. In particular our theoretical curves
are always close to the empirical ones, whereas the normal ones stay farther
or much farther from the empirical ones.

Figure 3 is designed for a big value of $r$ with $r=0.5$. Except the first two
pictures on the top row, which correspond to $(n,p)=(10,50)$ and
$(n,p)=(20,100)$, respectively, our theoretical curves match the empirical
ones well. As the explanation aforementioned, that the two simulations do not
behave well is understandable, simply because the sample sizes are very small.
Our theoretical curves are constantly close to the empirical ones, the normal
ones are nowhere close to the empirical ones.

In summary, the simulation indicates our approximation stated in Theorem 1
that
$(1-\sum_{i=1}^{\infty}\rho_{i}^{2})^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$
to $T_{SD}$ and $T_{p,1}$, respectively, outperforms the normal approximation
from Result 2.1.

Figure 2: Example 2.1. For small $r=0.1$, the empirical curves in blue are
close to our theoretical curves in red. As sample size $n$ becomes larger, our
approximation becomes better; the normal approximation (black curve) is no
longer valid.

Figure 3: Example 2.1. Compared to Figure 2 with $r=0.1$, we take $r=0.5$
here. The empirical curves in blue are close to our theoretical curves in red.
The normal approximation (black curve) from Result 2.1 is no longer valid.

###### Example 2.2

For $r\in(0,1)$, set $m=[p^{r}]$. Recall $\mathbb{A}_{m}$ from Example 2.1.
Define $\mathbb{R}$ by

$\displaystyle\mathbb{R}=\begin{pmatrix}\mathbb{A}_{m}&\mathbb{0}\\\
\mathbb{0}&\mathbb{I}_{p-m}\end{pmatrix}_{p\times p},$

where the two “$r$” from “$m=[p^{r}]$” and “$\mathbb{A}_{m}$” are the same
one. Assuming ${\boldsymbol{\mu}}=\boldsymbol{0}$, we have

$\displaystyle T_{p}\to\begin{cases}\xi_{0}&\text{if $0<r<\frac{1}{2}$};\\\
\frac{2}{\sqrt{5}}\xi_{0}+\frac{1}{\sqrt{10}}\cdot(\xi_{1}^{2}-1),&\text{if
$r=\frac{1}{2}$};\\\ \frac{1}{\sqrt{2}}\cdot(\xi_{1}^{2}-1),&\text{if
$\frac{1}{2}<r<1$}\end{cases}$ (12)

under condition $p=o(n^{2})$ for $0<r\leq\frac{1}{2}$ and $p=o(n^{1/(1-r)})$
for $\frac{1}{2}<r<1$, where $\xi_{0}$ and $\xi_{1}$ are i.i.d. $N(0,1)$ as
$r=\frac{1}{2}$. Here $T_{p}$ stands for $T_{SD}$ or $T_{p,1}$. Obviously,
there is a phase transition at $r=\frac{1}{2}$ as $r$ runs between $0$ and
$1$. The phase transition appears at both relative sizes of $n$ and $p$
together with the change of the values of matrix entries.

###### Example 2.3

Set $m=[\log(p+2)].$ Given $\tau\geq 0$, define integer
$p^{\prime}=p-[\tau\sqrt{p}]$, $\lambda_{i}=1+\tau
2^{-i}\big{(}1-2^{-m}\big{)}^{-1}\sqrt{p}$ for $1\leq i\leq m$ and
$\lambda_{i}=1$ for $i=m+1,\cdots,p^{\prime}-1$,
$\lambda_{p^{\prime}}=1+[\tau\sqrt{p}]-\tau\sqrt{p}\in[0,1]$, and the rest of
$\lambda_{i}$ are identical to $0$. By Lemma 4.3, there exists a correlation
matrix $\mathbb{R}$ such that $\mathbb{R}$ has eigenvalues
$\lambda_{i},\,1\leq i\leq p.$ Assuming ${\boldsymbol{\mu}}=\boldsymbol{0}$,
we have

$\displaystyle T_{p}\to\begin{cases}\xi_{0},&\text{if $\tau=0$};\\\
\sqrt{\frac{3}{\tau^{2}+3}}\xi_{0}+\sqrt{\frac{3\tau^{2}}{2(\tau^{2}+3)}}\sum_{i=1}^{\infty}\frac{1}{2^{i}}(\xi_{i}^{2}-1),&\text{if
$\tau\in(0,\infty)$};\\\
\sqrt{\frac{3}{2}}\sum_{i=1}^{\infty}\frac{1}{2^{i}}(\xi_{i}^{2}-1),&\text{if
$\tau\to\infty$}\end{cases}$ (13)

in distribution as $p\to\infty$ under condition $p=o(n^{2})$. Here
$\xi_{0},\xi_{1},\cdots,\xi_{d}$ are i.i.d. $N(0,1).$ Here $T_{p}$ stands for
$T_{SD}$ or $T_{p,1}$. The statistic $T_{p}$ behaves like a rainbow, which has
a Gaussian distribution at one end, and a mixing chi-squared distribution at
the horizon, and a hybrid in between.

An interesting remark is that, unlike the eigenvalues of a covariance matrix,
the eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{p}\geq 0$ of
the $p\times p$ correlation matrix $\mathbb{R}$ cannot be arbitrary. By Lemma
4.3, they have to satisfy the so-called “majorization” property:
$\lambda_{1}+\cdots+\lambda_{k}\geq k$ for every $1\leq k\leq p$ and
$\lambda_{1}+\cdots+\lambda_{p}=p$. Conversely, for any
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{p}\geq 0$ with the
“majorization” property, there is always a correlation matrix with eigenvalues
$\lambda_{i},\,1\leq i\leq p$.

##  3 Two Sample Mean Test for Large $p$ and Small $n$

We will study the two-sample mean testing problem in the high-dimensional
setting in this section. As before, assume that
$\\{{\boldsymbol{X}}_{i1},\cdots,{\boldsymbol{X}}_{in_{i}}\\}$ for $i=1,2$ are
two independent random samples with sizes $n_{1}$ and $n_{2}$, and from
$p$-variate normal distributions
$N_{p}({\boldsymbol{\mu}}_{1},{\boldsymbol{\Sigma}})$ and
$N_{p}({\boldsymbol{\mu}}_{2},{\boldsymbol{\Sigma}})$, respectively. We wish
to test

$\displaystyle H_{0}:{\boldsymbol{\mu}}_{1}={\boldsymbol{\mu}}_{2}\ \
\mbox{versus}\ \ H_{1}:{\boldsymbol{\mu}}_{1}\neq{\boldsymbol{\mu}}_{2}.$ (14)

Let $\bar{{\boldsymbol{X}}}_{i}$ be the sample mean of the $i$-th sample with
$i=1,2$ and $\hat{{\bf S}}$ be the pooled sample covariance matrix defined in
(3). Set $n=n_{1}+n_{2}$. Srivastava and Du (2008) consider the statistic
$T^{\prime}_{SD}$ from (4). Recall $\hat{{\bf D}}$ is the diagonal matrix of
$\hat{{\bf S}}$ in (3) and

$\displaystyle\hat{{\bf R}}=\hat{{\bf R}}_{p}=\hat{{\bf D}}^{-1/2}\hat{{\bf
S}}\hat{{\bf D}}^{-1/2}$ (15)

is the pooled sample correlation matrix. Similar to the discussion before (9)
on the modification of $T_{SD}$, we make a little change for $T_{SD}^{\prime}$
as follows. Define

$\displaystyle
T_{p,2}=\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-\frac{(n_{1}+n_{2}-1)p}{n_{1}+n_{2}-4}}{\sqrt{2\big{|}\mathrm{tr}(\hat{{\bf
R}}^{2})-\frac{p(p-1)}{n_{1}+n_{2}-2}\big{|}}}.$ (16)

Srivastava and Du (2008) derive the limiting distribution of $T_{SD}^{\prime}$
under conditions similar to those from Result 2.1 and an extra assumption
“$\frac{n_{1}}{n_{1}+n_{2}}\to c\in(0,1)$”. On the other hand, Hu and Bai
(2016) ask for the properties of $T_{SD}^{\prime}$ as the variances of the
entries of the population vector are quite different. We will provide a
solution for $T^{\prime}_{SD}$ and $T_{p,2}$ next. For mathematical rigor, we
assume that both sample sizes $n_{1}$ and $n_{2}$ depend on $p$.

######  THEOREM 2

Let $\\{{\boldsymbol{X}}_{i1},\cdots,{\boldsymbol{X}}_{in_{i}}\\}$ for $i=1,2$
be two independent random samples from
$N_{p}({\boldsymbol{\mu}}_{1},{\boldsymbol{\Sigma}})$ and
$N_{p}({\boldsymbol{\mu}}_{2},{\boldsymbol{\Sigma}})$, respectively. Let
$\mathbb{R}$ be as in (6) with eigenvalues
$\lambda_{1}\geq\cdots\geq\lambda_{p}\geq 0$. Assume

(a) $\lim_{p\to\infty}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}=\rho_{i}\geq 0$
for all $i\geq 1$;

(b) $\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2})\|\mathbb{R}\|_{F}}=0$ and
$\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2})^{a}}=0$ for some constant $a>0$.
If ${\boldsymbol{\mu}}_{1}={\boldsymbol{\mu}}_{2}$, then $T_{p,2}\to
b\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$ in
distribution, where $\xi_{0},\xi_{1},\xi_{2},\cdots$ are i.i.d. $N(0,1)$ and
$b=(1-\sum_{i=1}^{\infty}\rho_{i}^{2})^{1/2}$. The same conclusion also holds
for “ $T_{p,2}$”.

The above result does not require condition “$\min\\{n_{1},n_{2}\\}\to\infty$”
or “$\frac{n_{1}}{n_{1}+n_{2}}\to\gamma\in(0,1)$”, which are commonly assumed
in literature for the two-sample testing problem.

As discussed below Theorem 1, condition (8) implies that $\rho_{1}=0$ and
$\frac{p}{n\|\mathbb{R}\|_{F}}\to 0$ as $p\to\infty$. Condition
“$n_{1}+n_{2}=O(p^{\zeta})$, $\frac{1}{2}<\zeta\leq 1$” from (8) implies (b)
from Theorem 2. So our theorem is more general than the conclusion for the
two-sample test by Srivastava and Du (2008).

If $\mathbb{R}$ is $AR(1)$ model, the banded model, the Toeplitz or Hankel
matrices, please see the discussion below Theorem 1. They are also applied
here.

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
an $AR(1)$ model with $\mathbb{R}=(\gamma^{|i-j|})$ and the absolute values of
$\gamma=\gamma_{p}$ staying away from $1$. By using the Gersgorin disc theorem
[see, e.g., p. 344 from Horn and Johnson (2012)], the largest eigenvalue or
$\mathbb{R}$ is of order $O(1)$. Hence condition (a) of Theorem 1 holds with
$\rho_{1}=0$. If condition (b) also holds, then both $T_{SD}$ and $T_{p,1}$ go
to the standard normal distribution. The same conclusion is also valid for a
banded correlation matrix $\mathbb{R}=(r_{ij})$ with $r_{ij}=0$ for $|j-i|\geq
t$ where $t=t_{p}=o(\sqrt{p})$. In this case, the largest eigenvalue or
$\mathbb{R}$ is of order $o(\sqrt{p})$. Similar results can be obtained for
other patterned matrices including Toeplitz matrices, Hankel matrices and
symmetric circulant matrices; see, e.g., Brockwell et al. (2016).

Recall the definition of $T_{p,2}$ from (16). Let $\mathbb{D}$ be the diagonal
matrix of ${\boldsymbol{\Sigma}}$. Under some conditions including
$\|\mathbb{R}\|_{F}$ is of the order $p$, Zhang et al. (2020b) obtain the
limiting distribution of a normalized
$(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})^{T}{\bf
D}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})$ as
$\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$. Although this
quantity is not directly applicable because the unknown parameter matrix ${\bf
D}$ is involved, it is indeed suggestive. Their result seems to be a special
case of Theorem 2 with $\sum_{i=1}^{\infty}\rho_{i}^{2}=1$. In fact, under
their assumptions, we have confirmed that $\sum_{i=1}^{\infty}\rho_{i}^{2}=1$
in Lemma 6.37 from Section 6.4. Notice $\sqrt{p}\leq\|\mathbb{R}\|_{F}\leq p$
for any $\mathbb{R}$. Srivastava and Du (2008) consider the extreme case that
$\|\mathbb{R}\|_{F}$ is of the scale of $\sqrt{p}$. This is essentially a
weakly dependent situation with $\rho_{i}=0$ for each $i\geq 1$. According to
Theorem 2, the limiting distribution is a normal. On the other hand, roughly
speaking, Zhang et al. (2020b) intend to study another extreme case in which
$\|\mathbb{R}\|_{F}$ is of the scale $p$, and hence $\mathbb{R}$ is completely
a spiked model. As explained above, the corresponding limiting distribution is
a mixing chi-squared distribution. Our Theorem 2 handles the case for any
$\mathbb{R}$ with $\|\mathbb{R}\|_{F}$ running everywhere between $\sqrt{p}$
and $p$, and the limiting distribution turns out to be, interestingly enough,
a convolution of both.

Now we demonstrate our results by giving some examples.

###### Example 3.1

Recall Example 2.1 and $\mathbb{R}=\mathbb{A}_{p}$. For this example, Zhang et
al. (2020b) observe that “When $r=0.01$, the histograms are quite symmetric
and bell-shaped, indicating that a normal approximation as suggested by the
theory of Srivastava and Du (2008) may be applied for approximating the null
distribution of $T_{SD}$. However, when $r=0.5$ and $r=0.9$, the histograms
are quite skewed, indicating that a normal approximation is no longer
adequate.” In fact, our Theorem 2 explains their inspection very accurately as
follows. Assume $r=r_{p}$ and $\lim_{p\to\infty}\sqrt{p}\cdot r=c\geq 0$. By
changing “$n$” in Example 2.1 to “$n_{1}+n_{2}$”, we have

$\displaystyle T_{p}\to\begin{cases}\xi_{0},&\text{if $c=0$};\\\
\frac{1}{\sqrt{c^{2}+1}}\xi_{0}+\frac{c}{\sqrt{2(c^{2}+1)}}(\xi_{1}^{2}-1),&\text{if
$c\in(0,\infty)$};\\\ \frac{1}{\sqrt{2}}\cdot(\xi_{1}^{2}-1),&\text{if
$c=\infty$}\end{cases}$ (17)

under condition $p=o((n_{1}+n_{2})^{2})$ for the case $c\in[0,\infty)$ and
under the condition $p=o((n_{1}+n_{2})^{a})$ for some constant $a>0$ for the
case $c=\infty$, where $\xi_{0}$ and $\xi_{1}$ are i.i.d. $N(0,1).$ Here
$T_{p}$ stands for $T_{SD}^{\prime}$ or $T_{p,2}$. Obviously, there is a phase
transition between $c=0$ and $c=\infty$.

Figure 4: Example 3.1. Our theoretical curves always match the empirical ones
well no matter the former is a normal-like or chi-square-like curve, even one
of $n_{1}$ and $n_{2}$ being small. The normal curves by Srivastava and Du
(2008) are farther from the empirical ones as correlation $r$ increases.

Similar to Example 3.1, modifications of Examples 2.2 and 2.3 can also be
constructed for the two-sample case. It is straightforward, so we omit the
detail. Given a data set, the spikes of the population correlation
$\mathbb{R}$ have to be determined; see Fan et al. (2020) and Morales-Jimenez
et al. (2021).

Based on Example 3.1, the statistic $T_{p,2}$ from (16) is simulated for
$10^{5}$ times for each set of values of $(n_{1},n_{2},p,r)$. The value of $r$
is chosen at $0,1,0,5$ and $0.9$, designed for weakly dependence, dependence
and strong dependence, respectively. Our theoretical curves are normal-like,
mixing between $N(0,1)$ and chi-squared distribution and chi-squared-like
curves. For weakly dependent case, by examining the four graphs in the first
column, we see the normal approximation by Srivastava and Du (2008) becomes
less and less accurate as all of the values of $n_{1},n_{2},p$ increase.
However, our theoretical curve always match the empirical ones well. For
dependent and very dependent cases appeared in the second and third columns,
respectively, the normal curves are nowhere close to the empirical ones. It is
good to see our theoretical curves are always close to their empirical
counterparts.

##  4 A Law of Large Numbers for Sample Correlation Matrices

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with correlation matrix
$\mathbb{R}$. Let $\hat{\mathbb{R}}$ be the sample correlation matrix defined
in (7). One of the crucial steps in proving Theorems 1 and 2 is the use of an
asymptotically ratio-unbiased-estimator of $\mbox{tr}(\mathbb{R}^{2})$. In
their Lemma 3.2, Srivastava and Du (2008) state the following.

###### Result 4.1

If $n=O(p^{\zeta})$, $0<\zeta\leq 1$, under the condition (8),
$(\mbox{tr}(\mathbb{R}^{2})-p^{2}/n)/p$ converges to
$\lim_{p\to\infty}\mbox{tr}(\mathbb{R}^{2})/p$ in probability as
$(n,p)\to\infty$ and thus can be considered as a consistent estimator of
$\mbox{tr}(\mathbb{R}^{2})/p$ as $(n,p)\to\infty$.

The proof of the above result is give on pages 400-402 from Srivastava and Du
(2008). It lacks a mathematical rigor. By using a method from the random
matrix theory, we rigorously obtain a more general result than the above in
the following. Most efforts in the whole proofs of this paper are devoted to
this result, which is rather involved. The reason is that we assume no
independence at all among the entries of the population distribution, i.e.,
$\mathbb{R}$ is arbitrary.

######  THEOREM 3

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with correlation matrix
$\mathbb{R}$. Let $\hat{\mathbb{R}}$ be the sample correlation matrix defined
in (7). If $\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ and
$\lim_{p\to\infty}\frac{p}{n^{a}}=0$ for some constant $a>0$, then

$\displaystyle\frac{1}{\mbox{tr}(\mathbb{R}^{2})}\Big{[}\mbox{tr}(\hat{\mathbb{R}}^{2})-\frac{p(p-1)}{n-1}\Big{]}\to
1$

in probability as $p\to\infty.$

Unlike $\mbox{tr}(\mathbb{R}^{2})$, $\mbox{tr}(\mathbb{R}^{3})$ and
$\mbox{tr}(\mathbb{R}^{4})$ needed in Result 4.1, our conditions in the above
theorem are imposed on $\|\mathbb{R}\|_{F}=[\mbox{tr}(\mathbb{R}^{2})]^{1/2}$
only. Also, it is easy to see the difference between the assumptions on
restrictions of $n$ and $p$: Result 4.1 asks $n=O(p^{\zeta})$ for some
$0<\zeta\leq 1$ and our Theorem 3 holds as long as $n$ and $p$ are in the
order of a polynomial. We make some remarks on Theorem 3 next.

Naively, an obvious ratio-unbiased-estimator of $\mbox{tr}(\mathbb{R}^{2})$ is
its sample version $\mbox{tr}(\hat{\mathbb{R}}^{2})$. Theorem 3 indicates that
it is not true for the high dimensional setting. Instead, the modified version
“$\mbox{tr}(\hat{\mathbb{R}}^{2})-\frac{p(p-1)}{n-1}$” is a correct one. By
using the fact $\mbox{tr}(\hat{\mathbb{R}}^{2})\geq p$, we can see from
Theorem 3 that $\mbox{tr}(\hat{\mathbb{R}}^{2})$ is indeed a ratio-unbiased-
estimator of $\mbox{tr}(\mathbb{R}^{2})$ if $p$ is as small as $p=o(n)$, that
is, $\mbox{tr}(\hat{\mathbb{R}}^{2})/\mbox{tr}(\mathbb{R}^{2})\to 1$ in
probability.

Sample correlation matrices are a special type of random matrices. The
spectral distributions of eigenvalue values are studied in Jiang (2004) and
Xiao and Zhou (2010). The central limit theorems for determinants under
independent and correlated situations are obtained by Jiang and Yang (2013)
and Jiang (2009), respectively. The Tracy-Widom law is derived by Bao et al.
(2012). The central limit theorem for linear statistics of eigenvalues is
understood by Gao et al. (2017).

An interesting fact is the derivation of Theorem 3. Through the whole proof,
contrary to standard techniques to handle randomness, we do not use/assume any
independence from the population distribution
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$. In fact, in the most
dependent case, that is, all of the $p$ entries of
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ are identical, then
$\hat{\mathbb{R}}=\mathbb{R}$ and all of their $p^{2}$ entries are equal to
$1$, so $\mbox{tr}(\hat{\mathbb{R}}^{2})=\mbox{tr}(\mathbb{R}^{2})=p^{2}$. One
can see Theorem 3 trivially holds. On the other hand, for the most independent
case, i.e., if ${\boldsymbol{\Sigma}}=\mathbb{R}=\mathbb{I}_{p}$, then
$\mbox{tr}(\mathbb{R}^{2})=p.$ Assuming $p/n\to\rho\in(0,1]$, we have from
Theorem 3 that $\mbox{tr}(\hat{\mathbb{R}}^{2})/p\to 1+\rho.$ This actually
can be confirmed independently by a random matrix theory. In fact, let
$\mu_{p}$ be the empirical distribution of the eigenvalues of
$\hat{\mathbb{R}}$. By Theorem 2 from Jiang (2004), $\mu_{p}$ converges weakly
to the so-called Marčhenko-Pastur law $\mu$. The second moment of $\mu$ is
equal to $1+\rho$ [Lemma 3.1 from Bai and Silverstein (2010)]. Consequently,
$\mbox{tr}(\hat{\mathbb{R}}^{2})/p=\int x^{2}\,d\mu_{p}\to\int
x^{2}\,d\mu=1+\rho.$

The proof of Theorem 3 is a bit technical and lengthy. We first give an
accurate estimate of $E\mbox{tr}(\hat{\mathbb{R}}^{2})$ in Lemma 6.3. Then we
only need to show that its variance go to zero. To do so, since
$\hat{\mathbb{R}}^{2}=p+2\sum_{1\leq i<j\leq p}\hat{r}_{ij}^{2}$, we need to
understand the covariance between two sample correlations, say,
$\hat{r}_{ij}^{2}$ and $\hat{r}_{kl}^{2}$ for any $1\leq i,j,k,l\leq p$. The
precise result is given next.

######  THEOREM 4

Assume $m\geq 5$. Let
$\\{(X_{1j},X_{2j},X_{3j},X_{4j})^{T}\in\mathbb{R}^{4};\,1\leq j\leq m\\}$ be
i.i.d. random vectors with distribution $N_{4}(\mathbb{0},\mathbb{R})$, where
$\mathbb{R}=(r_{ij})_{4\times 4}$ and $r_{ii}=1$ for each $i$. Set

$\displaystyle\hat{r}_{ij}=\frac{\sum_{k=1}^{m}X_{ik}X_{jk}}{(\sum_{k=1}^{m}X_{ik}^{2})^{1/2}(\sum_{k=1}^{m}X_{jk}^{2})^{1/2}}$
(18)

for $1\leq i\leq j\leq 4.$ Then, for any $N\geq 1$,
$\mbox{Cov}(\hat{r}_{12}^{2},\hat{r}_{34}^{2})$ is equal to

$\displaystyle\varrho_{m,1}\cdot\sum_{1\leq i<j\leq
4}r_{ij}^{2}+\varrho_{m,2}\cdot r_{12}^{2}r_{34}^{2}$ $\displaystyle+$
$\displaystyle\varrho_{m,3}\cdot\big{(}r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31}\big{)}$
$\displaystyle+$
$\displaystyle\varrho_{m,4}\cdot\big{(}r_{13}r_{34}r_{41}+r_{23}r_{34}r_{42}\big{)}r_{12}^{2}$
$\displaystyle+$
$\displaystyle\varrho_{m,5}\cdot\big{(}r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}\big{)}r_{34}^{2}$
$\displaystyle+$ $\displaystyle\frac{\varrho_{m,6}}{m^{(N+1)/2}},$

where $\\{\varrho_{m,i};\,3\leq i\leq 5\\}$ are quantities not depending on
$\mathbb{R}$,

$\displaystyle|\varrho_{m,1}|\leq Km^{-2},\ \
|\varrho_{m,2}|\vee|\varrho_{m,3}|\vee|\varrho_{m,4}|\vee|\varrho_{m,5}|\leq
Km^{-1},\ |\varrho_{m,6}|\leq K$

and $K$ is a constant depending on $N$ but not on $m$ or $\mathbb{R}.$

The proof of Theorem 4 is rather involved. There are two reasons, one of them
is that no assumption is imposed on $\mathbb{R}$, as a consequence, we are not
able to use any techniques/theory related to independence. The second one is
that sample correlation coefficients are more difficult to be handled than
sample covariances. Look at (18) for $\hat{r}_{12}$ by taking $i=1$ and $j=2$.
From the law of large numbers, $(1/m)\sum_{k=1}^{m}X_{1k}^{2}\to 1$ and
$(1/m)\sum_{k=1}^{m}X_{2k}^{2}\to 1$. To understand $\hat{r}_{12}$, a naive
idea is to replace the denominator in the expression of $\hat{r}_{12}$ by $m$.
Interestingly, this works for the derivation of the Marčhenko-Pastur law
(Jiang, 2004) when $\mathbb{R}=\mathbb{I}$. Now we elaborate this more for a
further discussion. By the Taylor expansion, write

$\displaystyle\frac{1}{\sqrt{1+x}}=1-\frac{x}{2}+\frac{3}{8}x^{2}+\cdots+a_{N}x^{N}+o(x^{N+1})$
(19)

as $x$ is small. Jiang (2004) uses the above expansion with $N=0$ and $x$
taking $(1/m)\sum_{k=1}^{m}(X_{1k}^{2}-1)$ and
$(1/m)\sum_{k=1}^{m}(X_{2k}^{2}-1)$, respectively. If
$\mathbb{R}\neq\mathbb{I}$ but $\mathbb{R}=\mathbb{A}_{p}$ as in Example 2.1,
Fan and Jiang (2019) use the expansion (19) with $N=1$. In the proof our
Theorem 3, since no structure of $\mathbb{R}$ is assumed, we have to use (19)
for an arbitrary $N$. This is why the condition
“$\limsup_{p\to\infty}\frac{p}{n^{a}}=0$ for some constant $a>0$” is imposed
in Theorems 1 and 3. The new method here would be useful in the future for the
study on sample correlation matrices. It is worthwhile to mention that
handling sample correlation matrices is much more difficult than working on
sample covariance matrices. The difference is obvious: the entries of a sample
covariance matrix does not have the denominator as in the expression of
$\hat{r}_{ij}$ from (18). For instance, Bai and Saranadasa (1996) also
investigate a ratio-unbiased-estimator similar to Theorem 3 but for covariance
matrices, their argument is quick and short.

A byproduct of the proofs of Theorems 1 and 2 is the behavior of a quadratic
form related to the diagonal matrix of $\mathbb{R}$. It would be useful for
other research from this point forward. Assume
${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ is a random sample from
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with correlation matrix
$\mathbb{R}$. Review ${\bf D}$ is the diagonal matrix of
${\boldsymbol{\Sigma}}$ and $\hat{{\bf D}}$ is the diagonal matrix of
$\hat{{\bf S}}$ in (3).

######  LEMMA 4.1

Let $\boldsymbol{\eta}\sim N_{p}(\mathbb{0},{\boldsymbol{\Sigma}})$ and
$\boldsymbol{\eta}$ be independent of $\hat{{\bf S}}$. If
$\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$, then

$\displaystyle\frac{\boldsymbol{\eta}^{T}\hat{{\bf
D}}^{-1}\boldsymbol{\eta}}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}=\frac{\boldsymbol{\eta}^{T}{\bf
D}^{-1}\boldsymbol{\eta}}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}+o_{p}(1)$

as $p\to\infty.$

Since ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ is a random sample
from a multivariate normal distribution, $\bar{{\boldsymbol{X}}}$ is
independent of $\hat{{\bf S}}$, and hence $\bar{{\boldsymbol{X}}}$ is
independent of $\hat{{\bf D}}$. So the above conclusion holds if we take
$\boldsymbol{\eta}=\sqrt{p}\bar{{\boldsymbol{X}}}$. This is actually the way
we use this lemma in later proofs.

The limiting distributions in Theorems 1 and 2 essentially come from the
following.

######  LEMMA 4.2

For each $p\geq 1$, let $a_{p,1}\geq\cdots\geq a_{p,p}\geq 0$ be constants
with $a_{p,1}^{2}+\cdots+a_{p,p}^{2}=1$. Suppose
$\lim_{p\to\infty}a_{p,i}=\rho_{i}\geq 0$ for each $i\geq 1$. Let
$\xi_{1},\xi_{2},\cdots$ be i.i.d. with distribution $\chi^{2}(1)-1.$ Then
$a_{p,1}\xi_{1}+\cdots+a_{p,p}\xi_{p}$ converges to
$[2(1-\sum_{i=1}^{\infty}\rho_{i}^{2})]^{1/2}\eta+\sum_{i=1}^{\infty}\rho_{i}\xi_{i}$
in distribution as $p\to\infty$, where $\eta\sim N(0,1)$ and $\eta$ is
independent of $\\{\xi_{i};\,i\geq 1\\}$.

The comment after Example 2.3 is based on the following fact. It is
interesting useful in its own right.

######  LEMMA 4.3

Let $\mathbb{M}$ be a $p\times p$ correlation matrix, namely, $\mathbb{M}$ is
non-negative definite and its diagonal entries are all equal to $1$. Suppose
$\mathbb{M}$ has eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{p}$. Then
$\lambda_{1}+\cdots+\lambda_{k}\geq k$ for $k=1,\cdots,p$ and
$\lambda_{1}+\cdots+\lambda_{p}=p$. Conversely, for any
$\tau_{1}\geq\cdots\geq\tau_{p}\geq 0$ with $\tau_{1}+\cdots+\tau_{k}\geq k$
for each $k=1,\cdots,p$ and $\tau_{1}+\cdots+\tau_{p}=p$, there always exists
a correlation matrix with eigenvalues $\tau_{1},\cdots,\tau_{p}$.

##  5 Concluding Remarks and Discussion

1\. In this paper we have studied one-sample and two-sample mean tests. For
the multiple population case, it becomes the classical MANOVA problem. In the
“large $p$, small $n$” situation, Srivastava and Fujikoshi (2006) consider the
question by revising the classical $F$-test. They use functions of non-zero
eigenvalues of pseudo-MANOVA random matrix as test statistics. Zhang et al.
(2017) and Chen et al. (2019) study this case based on the $L_{2}$-type
statistics. An extension of our work to the MANOVA case under “large $p$,
small $n$” situation is meaningful. One direction is to modify the classical
tests, such as the Anderson test, the Pillal test, the Roy test, the Wilks
test and the Olson test, by a method similar to those from Srivastava and
Fujikoshi (2006). The classical statistics take the correlations among
population coordinates into account. It has kind of self-normalization, which
is favorable.

2\. All results in this paper are based on random samples from a multivariate
Gaussian distribution. Is it possible to generalize this to non-Gaussian
cases? The corresponding results are highly demanded simply because many data
are not sampled from a Gaussian distribution. Also, conclusions on non-
Gaussian scenario may help us understand robustness of our tests.

3\. We study the two-sample mean test by assuming the two populations have the
same covariance matrices. If the two covariance matrices are not identical,
this is the multivariate Behrens-Fisher problem. Various methods are proposed,
for example, by Yao (1965), Johansen (1980), Nel and Van der Merwe (1986),
Krishnamoorthy and Yu (2004), Chen and Qin (2010) and Chen et al. (2019). For
one-dimensional case ($p=1$), the method initiated by Welch (1947) is probably
the most popular one. The author basically normalize the difference of sample
means by its sample standard deviation. Srivastava et al. (2013) propose a
statistic by replacing $\hat{\mathbb{D}}$ from (16) with a hybrid of two
samples. Similar to the study in this paper, its properties are needed to be
understood, too.

4\. As far as proofs go, we spend most of our energy proving Theorem 4. We do
not use complicated technology. Instead, by we use brute force to compute
mixed moments of multiple Gaussian random variables in combination with some
random matrix theory. As a result, the argument is lengthy. It is possible to
shorten the proof. To get the law of large numbers stated in Theorem 4, one
may like to try Gaussian concentration inequalities; see, for example, Ledoux
(2001) and Boucheron et al. (2013). An alternative way is to get the joint
density function of $\hat{r}_{12}$ and $\hat{r}_{34}$ similar to that of the
marginal density of $\hat{r}_{12}$ treated at (26) through an hypergeometric
function.

5\. Theorems 1 and 2 present the null distributions for the one-sample and
two-sample tests, respectively. We actually have tried to derive the non-null
limiting distribution to explore power functions of the tests. The argument is
also very involved. The current paper is already very lengthy, so we postpone
and leave it as a future work.

## Acknowledgment

The bulk of the work of Tiefeng Jiang was conducted during his visit at Baidu
Research in 2019-2020. After completing an initial draft, we solicited
comments from several colleagues and we appreciate their feedback.

##  6 Proofs

One of the main steps of proving Theorems 1 and 2 is a weak law of large
numbers for $\mbox{tr}(\hat{\mathbb{R}}^{2})$ (Theorem 3), where
$\hat{\mathbb{R}}=\hat{\mathbb{R}}_{p}$ is the sample correlation matrix
defined in (7). To derive this, we need to study the mean and variance of
$\mbox{tr}(\hat{\mathbb{R}}^{2})$. In Section 6.1, we will get an accurate
estimate of $E[\mbox{tr}(\hat{\mathbb{R}}^{2})]$ by an argument on
hypergeometric functions. Then we will show
$\mbox{Var}(\mbox{tr}(\hat{\mathbb{R}}^{2}))$ go to zero. The proof of Theorem
3 will be completed once this step is established. This is a rather involved
step and it will be understood step by step in Sections 6.2.1-6.2.6.

###  6.1 Evaluation of the Mean of the Frobenius Norm of a Sample Correlation
matrix

In this section, we will work on the mean of
$\mbox{tr}(\hat{\mathbb{R}}^{2})$, where
$\hat{\mathbb{R}}=\hat{\mathbb{R}}_{p}$ is the sample correlation matrix
defined in (7). The critical tool is the hypergeometric function. We first
need a preliminary result as follows.

######  LEMMA 6.1

There exists a constant $C>0$ such that the following holds for all $k\geq 2$.

$\displaystyle\sum_{n=1}^{\infty}\binom{n+k}{k}^{-1}\leq Ck^{-1/4}.$

Proof of Lemma 6.1. For given $a>0$, set $f_{a}(x)=(1+ax^{-1})^{-x}$ for
$x>0$. We claim that $f_{a}(x)$ is decreasing in $x\in(0,\infty).$ In fact,
let $g(x)=\log f_{a}(x)$. Then

$\displaystyle
g^{\prime}(x)=\Big{[}-x\log\Big{(}1+\frac{a}{x}\Big{)}\Big{]}^{\prime}=-\log\Big{(}1+\frac{a}{x}\Big{)}-x\frac{-\frac{a}{x^{2}}}{1+\frac{a}{x}}.$

Then

$\displaystyle
g^{\prime}(x)=\log\Big{(}1-\frac{a}{x+a}\Big{)}+\frac{a}{x+a}\leq 0$

for all $x>0$ since $\log(1+y)\leq y$ for all $y>-1$. Hence, $f_{a}(x)$ is
decreasing in $x\in(0,\infty)$ for each $a>0.$ By the Stirling formula,
$m!=\sqrt{2\pi m}m^{m}e^{-m+\frac{\theta_{m}}{12m}}$ with $\theta_{m}\in(0,1)$
for each $m\geq 1$. Consequently,

$\displaystyle\sqrt{2\pi m}m^{m}e^{-m}\leq m!\leq e\cdot\sqrt{2\pi
m}m^{m}e^{-m}$

for all $m\geq 1.$ Therefore,

$\displaystyle\binom{n+k}{k}=\frac{(n+k)!}{n!k!}$ $\displaystyle\geq$
$\displaystyle e^{-2}\cdot\frac{\sqrt{n+k}}{\sqrt{2\pi
nk}}\cdot\frac{(n+k)^{n+k}e^{-(n+k)}}{k^{k}n^{n}e^{-(n+k)}}$
$\displaystyle\geq$ $\displaystyle e^{-2}\frac{1}{\sqrt{2\pi
k}}\cdot\frac{(n+k)^{n+k}}{k^{k}n^{n}}.$

This implies that

$\displaystyle\frac{1}{\binom{n+k}{k}}\leq(2\pi
e^{2}k^{1/2})\Big{(}1+\frac{k}{n}\Big{)}^{-n}\Big{(}1+\frac{n}{k}\Big{)}^{-k}.$
(20)

For fixed $k$, since $f_{k}(x)$ is decreasing in $x$, we have

$\displaystyle\sum_{n=1}^{k}\Big{(}1+\frac{k}{n}\Big{)}^{-n}\Big{(}1+\frac{n}{k}\Big{)}^{-k}$
$\displaystyle\leq$
$\displaystyle\sum_{n=1}^{k}\Big{(}1+\frac{k}{n}\Big{)}^{-n}$ (21)
$\displaystyle\leq$
$\displaystyle\frac{1}{k+1}+\sum_{n=2}^{k}\int_{n-1}^{n}\Big{(}1+\frac{k}{x}\Big{)}^{-x}\,dx$
$\displaystyle=$
$\displaystyle\frac{1}{k+1}+\int_{1}^{k}\Big{(}1+\frac{k}{x}\Big{)}^{-x}\,dx.$

Again, use the fact that $f_{k}(x)$ is decreasing to see

$\displaystyle\int_{1}^{k}\Big{(}1+\frac{k}{x}\Big{)}^{-x}\,dx$
$\displaystyle=$
$\displaystyle\int_{1}^{k^{1/4}}\Big{(}1+\frac{k}{x}\Big{)}^{-x}\,dx+\int_{k^{1/4}}^{k}\Big{(}1+\frac{k}{x}\Big{)}^{-x}\,dx$
$\displaystyle\leq$
$\displaystyle\frac{k^{1/4}}{k+1}+k\cdot(1+k^{3/4})^{-k^{1/4}}.$

If $k\geq 81$, then $k^{1/4}\geq 3$. Consequently,
$k\cdot(1+k^{3/4})^{-k^{1/4}}\leq k\cdot(k^{3/4})^{-3}=k^{-5/4}.$ From (21) we
see

$\displaystyle\sum_{n=1}^{k}\Big{(}1+\frac{k}{n}\Big{)}^{-n}\Big{(}1+\frac{n}{k}\Big{)}^{-k}\leq\frac{1}{k+1}+\frac{k^{1/4}}{k+1}+\frac{1}{k^{5/4}}\leq\frac{3}{k^{3/4}}.$
(22)

Now, $(1+\frac{k}{n})^{-n}\leq 2^{-k}$ as $n\geq k$ by the fact $f_{k}(x)$ is
decreasing. Also, $(1+\frac{n}{k})^{-k}\leq(1+\frac{n}{k})^{-81}$ as $k\geq
81$. Hence

$\displaystyle\sum_{n=k}^{\infty}\Big{(}1+\frac{k}{n}\Big{)}^{-n}\Big{(}1+\frac{n}{k}\Big{)}^{-k}\leq
2^{-k}\sum_{n=k}^{\infty}\Big{(}1+\frac{n}{k}\Big{)}^{-81}\leq
2^{-k}k^{81}\sum_{n=1}^{\infty}\frac{1}{n^{81}}$ (23)

as $k\geq 81.$ This joined with (20), (22) and (23) yields that

$\displaystyle\sum_{n=1}^{\infty}\frac{k^{1/4}}{\binom{n+k}{k}}\leq
k^{1/4}\cdot(2\pi
e^{2}k^{1/2})\cdot\big{(}3k^{-3/4}+\zeta(81)2^{-k}k^{81}\big{)}\leq
C^{\prime}$

for all $k\geq 81$, where
$\zeta(81)=\sum_{n=1}^{\infty}\frac{1}{n^{81}}<\infty$ is the Riemann zeta
function evaluated at $81$ and $C^{\prime}$ is another numerical constant not
depending on $k$. In summary,

$\displaystyle\sup_{k\geq
81}\sum_{n=1}^{\infty}k^{1/4}\binom{n+k}{k}^{-1}<\infty.$ (24)

Note that

$\displaystyle\frac{\binom{n+i+1}{i+1}}{\binom{n+i}{i}}=\frac{n+i+1}{i+1}>1$

for any $i\geq 1$ and $n\geq 1$. We know $\binom{n+i}{i}$ is increasing in
$i$, and hence

$\displaystyle\sup_{2\leq k\leq
80}\sum_{n=1}^{\infty}k^{1/4}\binom{n+k}{k}^{-1}\leq
80^{1/4}\sum_{n=1}^{\infty}\binom{n+2}{2}^{-1}\leq 2\cdot
80^{1/4}\sum_{n=1}^{\infty}\frac{1}{n^{2}}<\infty.$

This and (24) imply the conclusion. $\square$

The following result quantifies the second moment of the sample correlation
coefficient of a random sample of size $m$ up to an error $O(m^{-1/4}).$

######  LEMMA 6.2

Let $\\{(X_{i},Y_{i})^{T};\,1\leq i\leq m\\}$ be i.i.d. $2$-dimensional normal
random vectors with $EX_{1}=EY_{1}=0$, $EX_{1}^{2}=EY_{1}^{2}=1$ and
$\mbox{Cov}(X_{1},Y_{1})=r$. Set

$\displaystyle\hat{r}_{m}=\frac{\sum_{k=1}^{m}X_{i}Y_{i}}{(\sum_{k=1}^{m}X_{i}^{2})^{1/2}(\sum_{k=1}^{m}Y_{i}^{2})^{1/2}}.$

Write

$\displaystyle E\big{(}\hat{r}_{m}^{2})=\frac{1}{m}+r^{2}+b_{m}(r)\cdot r^{2}$

for $m\geq 4$. Then $\sup_{m\geq 4,\,|r|\leq 1}|m^{1/4}b_{m}(r)|<\infty$.

Proof of Lemma 6.2. If $r=\pm 1$, since $EX_{1}=EY_{1}=0$ and
$EX_{1}^{2}=EY_{1}^{2}=1$, we know $Y_{i}=\pm X_{i}$ for each $i$. By the
definition of $\hat{r}_{m}$, trivially, $\hat{r}_{m}=\pm 1$. This implies that

$\displaystyle b_{m}(\pm 1)=-\frac{1}{m}$ (25)

for each $m\geq 1.$ In the following we always assume $r^{2}<1.$ From Ghosh
(1966) or p. 156 in Muirhead (1982),

$\displaystyle
E\big{(}\hat{r}_{m}^{2})=1-\frac{m-1}{m}(1-r^{2})\cdot{}_{2}F_{1}\Big{(}1,1;\frac{1}{2}m+1;r^{2}\Big{)},$
(26)

where ${}_{2}F_{1}(a,b;c;z)$ is an hypergeometric function defined by

$\displaystyle{}_{2}F_{1}(a,b;c;z)=\sum_{k=0}^{\infty}\frac{(a)_{k}(b)_{k}}{(c)_{k}k!}z^{k},~{}~{}|z|<1,$

and where $(x)_{0}=1$ and $(x)_{k}=x(x+1)\cdots(x+k-1)$ for $k\geq 1$; see,
for example, p. 20 from Muirhead (1982). Notice $(1)_{k}=k!$ and

$\displaystyle\Big{(}\frac{1}{2}m+1\Big{)}_{k}=\frac{m+2}{2}\cdot\frac{m+4}{2}\cdots\frac{m+2k}{2}.$

Thus

$\displaystyle{}_{2}F_{1}\Big{(}1,1;\frac{1}{2}m+1;r^{2}\Big{)}$
$\displaystyle=$ $\displaystyle
1+\sum_{k=1}^{\infty}\frac{2^{k}k!}{(m+2)(m+4)\cdots(m+2k)}r^{2k}$ (27)
$\displaystyle=$ $\displaystyle
1+r^{2}\sum_{k=1}^{\infty}\frac{2^{k}k!}{(m+2)(m+4)\cdots(m+2k)}r^{2k-2}.$

Evidently, the last sum is bounded by

$\displaystyle\sum_{k=1}^{\infty}\frac{2^{k}k!}{(m+2)(m+4)\cdots(m+2k)}.$

Set $j=[\frac{m}{2}]$, where $[x]$ denotes the integer part of a real number
$x\geq 0$. Then $j\leq\frac{m}{2}$, or equivalently, $m\geq 2j$. It follows
that

$\displaystyle\frac{2^{k}k!}{(m+2)(m+4)\cdots(m+2k)}\leq\frac{k!}{(j+1)(j+2)\cdots(j+k)}=\frac{j!k!}{(j+k)!}.$

It follows that

$\displaystyle\sum_{k=1}^{\infty}\frac{2^{k}k!}{(m+2)(m+4)\cdots(m+2k)}\leq\sum_{k=1}^{\infty}\frac{1}{\binom{j+k}{j}}\leq\frac{K_{1}}{j^{1/4}}\leq\frac{K_{2}}{m^{1/4}}$

as $j\geq 2$, or equivalently, $m\geq 4$ by Lemma 6.1, where $K_{1}$ and
$K_{2}$ are constants not depending on $m$. From this and (27), we are able to
write

$\displaystyle{}_{2}F_{1}\Big{(}1,1;\frac{1}{2}m+1;r^{2}\Big{)}=1+a_{m}r^{2},$

where $0\leq a_{m}\leq K_{2}m^{-1/4}.$ Substitute this back to (26), we get

$\displaystyle E\big{(}\hat{r}_{m}^{2})$ $\displaystyle=$ $\displaystyle
1-\frac{m-1}{m}(1-r^{2})\cdot(1+a_{m}r^{2})$ $\displaystyle=$ $\displaystyle
1-\frac{m-1}{m}(1-r^{2})-\frac{m-1}{m}(1-r^{2})\cdot a_{m}r^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{m}+r^{2}+b_{m}(r)r^{2},$

where

$\displaystyle b_{m}(r)=-\frac{1}{m}r^{2}-\frac{m-1}{m}(1-r^{2})\cdot
a_{m}r^{2}.$

Obviously, $\sup_{|r|\leq 1}|b_{m}(r)|\leq K_{3}m^{-1/4}$ for every $m\geq 4$,
where $K_{3}$ is a constant not depending on $m$. This and (25) lead to the
desired conclusion. $\square$

######  LEMMA 6.3

Let ${\boldsymbol{X}}_{1},\cdots,{\boldsymbol{X}}_{n}$ be a random sample from
$N_{p}({\boldsymbol{\mu}},{\boldsymbol{\Sigma}})$ with correlation matrix
$\mathbb{R}$. Let $\hat{\mathbb{R}}=\hat{\mathbb{R}}_{p}$ be the sample
correlation matrix defined in (7). If
$\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ and
$\limsup_{p\to\infty}\frac{p}{n^{a}}=0$ for some constant $a>0$, then, as
$p\to\infty$,

$\displaystyle
E\,\mbox{tr}(\hat{\mathbb{R}}_{p}^{2})=\frac{p(p-1)}{n-1}+\mbox{tr}(\mathbb{R}_{p}^{2})\cdot\big{[}1+O(m^{-1/4})\big{]}.$

Proof of Lemma 6.3. Set $m=n-1$. The notation $W_{p}(m,{\boldsymbol{\Sigma}})$
represents the distribution of the Wishart matrix $\mathbb{U}^{T}\mathbb{U}$,
where $\mathbb{U}$ is an $m\times p$ matrix whose rows are i.i.d. with
distribution $N_{p}(\mathbb{0},{\boldsymbol{\Sigma}})$. Then, $n\hat{{\bf S}}$
has the Wishart distribution $W_{p}(m,{\boldsymbol{\Sigma}})$; see, for
example, Theorem 3.1.2 from Muirhead (1982). That is, $n\hat{{\bf
S}}\overset{d}{=}\mathbb{U}^{T}\mathbb{U}$. We claim

$\displaystyle\hat{\mathbb{R}}_{p}=(\hat{r}_{ij})_{p\times p}\
\overset{d}{=}\Big{(}\frac{\mathbb{v}_{i}^{T}\mathbb{v}_{j}}{\|\mathbb{v}_{i}^{T}\|\cdot\|\mathbb{v}_{j}\|}\Big{)}_{p\times
p},$ (28)

where the $m$ rows of $(\mathbb{v}_{1},\cdots,\mathbb{v}_{p})_{m\times p}$ are
i.i.d. with distribution $N_{p}(\mathbb{0},\mathbb{R})$. In fact, write
$\mathbb{U}=(u_{ij})=(\mathbb{u}_{1},\cdots,\mathbb{u}_{p})$ where
$\mathbb{u}_{i}=(u_{1i},\cdots,u_{mi})^{T}\in\mathbb{R}^{m}$ for each $i$.
Then $\mathbb{u}_{i}^{T}\mathbb{u}_{j}$ is the $(i,j)$-entry of
$\mathbb{U}^{T}\mathbb{U}$ and the diagonal entries are
$\|\mathbb{u}_{i}\|^{2}$ for $1\leq i\leq p$. We know from (7) that $\hat{{\bf
R}}_{p}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2}$, where
$\hat{{\bf D}}$ is the diagonal matrix of $\hat{{\bf S}}$. Then

$\displaystyle\hat{\mathbb{R}}_{p}=(\hat{r}_{ij})_{p\times p}\
\overset{d}{=}\Big{(}\frac{\mathbb{u}_{i}^{T}\mathbb{u}_{j}}{\|\mathbb{u}_{i}^{T}\|\cdot\|\mathbb{u}_{j}\|}\Big{)}_{p\times
p}.$

Notice $(\mathbb{u}_{i},\mathbb{u}_{j})=(u_{ki},u_{kj})_{1\leq k\leq m}$ for
any $i<j$ and the $m$ rows are i.i.d. bivariate normal random variables with
mean vector $\mathbb{0}$ and covariance matrix

$\displaystyle\begin{pmatrix}\sigma_{ii}&\sigma_{ij}\\\
\sigma_{ij}&\sigma_{jj}\end{pmatrix}$

where ${\boldsymbol{\Sigma}}:=(\sigma_{ij})_{p\times p}$. Write
$\mathbb{R}_{p}=(r_{ij})_{p\times p}$. Then
$r_{ij}=\sigma_{ij}(\sigma_{ii}\sigma_{jj})^{-1/2}$ by definition. Observe

$\displaystyle\frac{\mathbb{u}_{i}^{T}\mathbb{u}_{j}}{\|\mathbb{u}_{i}^{T}\|\cdot\|\mathbb{u}_{j}\|}=\frac{(\sigma_{ii}^{-1/2}\mathbb{u}_{i})^{T}(\sigma_{jj}^{-1/2}\mathbb{u}_{j})}{\|\sigma_{ii}^{-1/2}\mathbb{u}_{i}^{T}\|\cdot\|\sigma_{jj}^{-1/2}\mathbb{u}_{j}\|}.$

Take $\mathbb{v}_{i}=\sigma_{ii}^{-1/2}\mathbb{u}_{i}$ to obtain (28).

By Lemma 6.2,

$\displaystyle E\,\mbox{tr}(\hat{\mathbb{R}}_{p}^{2})$ $\displaystyle=$
$\displaystyle p+2\sum_{1\leq i<j\leq p}E\big{(}\hat{r}_{ij}^{2}\big{)}$
$\displaystyle=$ $\displaystyle p+2\sum_{1\leq i<j\leq
p}\Big{[}\frac{1}{m}+r_{ij}^{2}+b_{m}(r_{ij})r_{ij}^{2}\Big{]}$
$\displaystyle=$
$\displaystyle\frac{p(p-1)}{m}+\mbox{tr}(\mathbb{R}_{p}^{2})+2\sum_{1\leq
i<j\leq p}b_{m}(r_{ij})r_{ij}^{2},$

where $\sup_{|r|\leq 1}|b_{m}(r)|\leq Km^{-1/4}$ and $K$ is a constant not
depending on $m$ or $r$. Since

$\displaystyle 2\Big{|}\sum_{1\leq i<j\leq
p}b_{m}(r_{ij})r_{ij}^{2}\Big{|}\leq(2K)m^{-1/4}\cdot\mbox{tr}(\mathbb{R}_{p}^{2}),$

the above two displays show that

$\displaystyle
E\,\mbox{tr}(\hat{\mathbb{R}}_{p}^{2})=\frac{p(p-1)}{n-1}+\mbox{tr}(\mathbb{R}_{p}^{2})\cdot\big{[}1+O(m^{-1/4})\big{]}$

as $p\to\infty.$ The proof is completed. $\square$

###  6.2 The Proofs of Theorems 3 and 4

####  6.2.1 Mixing Moments of Gaussian Random Variables

The following lemma is a very useful tool to compute the mean of the product
of Gaussian random variables.

######  LEMMA 6.4

(Isserlis and Wick formula) Let $n\geq 2$ be an even integer and
$(X_{1},\cdots,X_{n})^{T}\in\mathbb{R}^{n}$ follows a normal distribution with
mean vector $\mathbb{0}$. Then

$\displaystyle E(X_{1}X_{2}\cdots X_{n})=\sum_{p}\prod_{(i,j)\in
p}E(X_{i}X_{j}),$

where the sum runs over every pairing $p$ of $\\{1,2,\cdots,n\\}$, that is,
all distinct ways of partitioning $\\{1,2,\cdots,n\\}$ into pairs $\\{i,j\\}$,
and the product is over the pairs contained in $p$. Sometimes we also use its
equivalent form:

$\displaystyle E(X_{1}X_{2}\cdots X_{n})=\sum_{i=2}^{n}E(X_{1}X_{i})\cdot
E\Big{(}\frac{X_{2}\cdots X_{n}}{X_{i}}\Big{)}.$

A seemingly more general, but actually an equivalent formula of Lemma 6.4 is
the following Lemma 6.5, which will be used only once in a comment after the
proof of Lemma 6.7.

######  LEMMA 6.5

(Guiard, 1986) Let $X=(X_{1},\cdots,X_{n})$ be n-dimensionally centralized
normally distributed with $Cov(X_{i},X_{j})=\sigma_{ij}$. Let
$\\{\alpha_{i};\,1\leq i\leq n\\}$ be positive integers. If
$\sum_{i=1}^{n}\alpha_{i}$ is odd, then
$\mathbb{E}(\prod_{i=1}^{n}X_{i}^{\alpha_{i}})=0.$ If
$\sum_{i=1}^{n}\alpha_{i}$ is even, then

$\displaystyle\mathbb{E}\left(\prod_{i=1}^{n}X_{i}^{\alpha_{i}}\right)=\sum_{2\beta_{ii}+\sum_{j:j\neq
i}\beta_{ij}=\alpha_{i},\forall
i}\dfrac{\prod_{i=1}^{n}\alpha_{i}!}{\prod_{i=1}^{n}2^{\beta_{ii}}\prod_{1\leq
i\leq j\leq n}\beta_{ij}!}\prod_{1\leq i\leq j\leq
n}\sigma_{ij}^{\beta_{ij}},$

where $\\{\beta_{ij},1\leq i\leq j\leq n\\}$ are non-negative integers.

In Lemmas 6.6 and 6.7 below, two identities on mixing moments of Gaussian
random variables will be presented.

######  LEMMA 6.6

Let $(X_{1},X_{2},X_{3},X_{4})^{T}\sim N_{4}(0,\mathbb{R})$. Assume
$\mathbb{R}=(r_{ij})_{4\times 4}$ with $r_{ii}=1$ for each $i$. Then

(i)
$E(X_{1}X_{2}X_{3}^{2}X_{4}^{2})=r_{12}+2r_{12}r_{34}^{2}+2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}$;

(ii)
$E(X_{1}^{2}X_{2}^{2}X_{3}^{2})=1+2r_{12}^{2}+2r_{13}^{2}+2r_{23}^{2}+8r_{12}r_{23}r_{31}.$

Proof of Lemma 6.6. Let $(\xi_{1},\cdots,\xi_{6})^{T}$ be a multivariate
normal with mean vector $\mathbb{0}$,
$\mbox{Cov}(\xi_{i},\xi_{j})=\sigma_{ij}$ and $\mbox{Var}(\xi_{i})=1$ for each
$i$. Then, by Lemma 6.4,

$\displaystyle E(\xi_{1}\xi_{2}\xi_{3}\xi_{4}\xi_{5}\xi_{6})$ $\displaystyle=$
$\displaystyle\sigma_{12}\sigma_{34}\sigma_{56}+\sigma_{12}\sigma_{35}\sigma_{46}+\sigma_{12}\sigma_{36}\sigma_{45}+$
$\displaystyle\sigma_{13}\sigma_{24}\sigma_{56}+\sigma_{13}\sigma_{25}\sigma_{46}+\sigma_{13}\sigma_{26}\sigma_{45}+$
$\displaystyle\sigma_{14}\sigma_{23}\sigma_{56}+\sigma_{14}\sigma_{25}\sigma_{36}+\sigma_{14}\sigma_{26}\sigma_{35}+$
$\displaystyle\sigma_{15}\sigma_{23}\sigma_{46}+\sigma_{15}\sigma_{24}\sigma_{36}+\sigma_{15}\sigma_{26}\sigma_{34}+$
$\displaystyle\sigma_{16}\sigma_{23}\sigma_{45}+\sigma_{16}\sigma_{24}\sigma_{35}+\sigma_{16}\sigma_{25}\sigma_{34}.$

Now take $\xi_{3}=\xi_{4}$ and $\xi_{5}=\xi_{6}$ to see

$\displaystyle E(\xi_{1}\xi_{2}\xi_{3}^{2}\xi_{5}^{2})$ $\displaystyle=$
$\displaystyle\sigma_{12}+\sigma_{12}\sigma_{35}^{2}+\sigma_{12}\sigma_{35}^{2}+$
$\displaystyle\sigma_{13}\sigma_{23}+\sigma_{13}\sigma_{25}\sigma_{35}+\sigma_{13}\sigma_{25}\sigma_{35}+$
$\displaystyle\sigma_{13}\sigma_{23}+\sigma_{13}\sigma_{25}\sigma_{35}+\sigma_{13}\sigma_{25}\sigma_{35}+$
$\displaystyle\sigma_{15}\sigma_{23}\sigma_{35}+\sigma_{15}\sigma_{23}\sigma_{35}+\sigma_{15}\sigma_{25}+$
$\displaystyle\sigma_{15}\sigma_{23}\sigma_{35}+\sigma_{15}\sigma_{23}\sigma_{35}+\sigma_{15}\sigma_{25}$
$\displaystyle=$
$\displaystyle\sigma_{12}+2\sigma_{12}\sigma_{35}^{2}+2\sigma_{13}\sigma_{23}+2\sigma_{15}\sigma_{25}+4\sigma_{13}\sigma_{25}\sigma_{35}+4\sigma_{15}\sigma_{23}\sigma_{35}.$

This says that

$\displaystyle
E(X_{1}X_{2}X_{3}^{2}X_{4}^{2})=r_{12}+2r_{12}r_{34}^{2}+2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}.$
(29)

We obtain (i). Now, take $X_{2}=X_{1}$ in (29). By using the fact $r_{12}=1$
and by changing “$2$” to “$1$” to the indices of $r$ we obtain

$\displaystyle E(X_{1}^{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle=$ $\displaystyle
1+2r_{34}^{2}+2r_{13}^{2}+2r_{14}^{2}+4r_{13}r_{14}r_{34}+4r_{14}r_{13}r_{34}$
$\displaystyle=$ $\displaystyle
1+2r_{13}^{2}+2r_{14}^{2}+2r_{34}^{2}+8r_{13}r_{14}r_{34}.$

In the above and change “$3$” to “$2$” and “$4$” to “$3$” to see

$\displaystyle
E(X_{1}^{2}X_{2}^{2}X_{3}^{2})=1+2r_{12}^{2}+2r_{13}^{2}+2r_{23}^{2}+8r_{12}r_{13}r_{23}.$

The proof is completed. $\square$

######  LEMMA 6.7

Suppose the 4-dimensional random vector $(X_{1},X_{2},X_{3},X_{4})^{T}\sim
N_{4}(0,\mathbb{R})$. Assume $\mathbb{R}=(r_{ij})_{4\times 4}$ with $r_{ii}=1$
for each $i$. Then

$\displaystyle E(X_{1}^{2}X_{2}^{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle=$
$\displaystyle 1+2\sum_{1\leq i<j\leq
4}r_{ij}^{2}+4(r_{12}^{2}r_{34}^{2}+r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle\ \
+8(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle\ \
+16(r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31}+r_{13}r_{32}r_{24}r_{41}).$

In particular, take $X_{3}=X_{1}$ and $X_{4}=X_{2}$ to see
$E(X_{1}^{4}X_{2}^{4})=9+72r_{12}^{2}+24r_{12}^{4}.$

Proof of Lemma 6.7. By Lemma 6.4,

$\displaystyle E(X_{1}^{2}X_{2}^{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle=$
$\displaystyle E(X_{1}^{2})\cdot
E(X_{2}^{2}X_{3}^{2}X_{4}^{2})+2E(X_{1}X_{2})\cdot
E(X_{1}X_{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle+2E(X_{1}X_{3})\cdot
E(X_{1}X_{3}X_{2}^{2}X_{4}^{2})+2E(X_{1}X_{4})\cdot
E(X_{1}X_{4}X_{2}^{2}X_{3}^{2}).$

By Lemma 6.6(i),

$\displaystyle
E(X_{1}X_{2}X_{3}^{2}X_{4}^{2})=r_{12}+2r_{12}r_{34}^{2}+2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}.$
(30)

Similarly, we obtain $E(X_{1}X_{3}X_{2}^{2}X_{4}^{2})$ and
$E(X_{1}X_{4}X_{2}^{2}X_{3}^{2})$ by exchanging the roles of “$X_{2}$” and
“$X_{3}$” and exchanging the roles of “$X_{2}$” and “$X_{4}$”, respectively,
from (30). By Lemma 6.6(ii),

$\displaystyle
E(X_{1}^{2}X_{2}^{2}X_{3}^{2})=1+2r_{12}^{2}+2r_{13}^{2}+2r_{23}^{2}+8r_{12}r_{13}r_{23}.$

Hence,

$\displaystyle E(X_{1}^{2}X_{2}^{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle=$
$\displaystyle 1+2r_{23}^{2}+2r_{24}^{2}+2r_{34}^{2}+8r_{23}r_{24}r_{34}+$
$\displaystyle
2r_{12}\big{(}r_{12}+2r_{12}r_{34}^{2}+2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}\big{)}+$
$\displaystyle
2r_{13}\big{(}r_{13}+2r_{13}r_{24}^{2}+2r_{12}r_{23}+2r_{14}r_{34}+4r_{12}r_{34}r_{24}+4r_{14}r_{23}r_{24}\big{)}+$
$\displaystyle
2r_{14}\big{(}r_{14}+2r_{14}r_{23}^{2}+2r_{13}r_{34}+2r_{12}r_{24}+4r_{13}r_{24}r_{23}+4r_{12}r_{34}r_{23}\big{)}.$

Sorting them out, we have

$\displaystyle E(X_{1}^{2}X_{2}^{2}X_{3}^{2}X_{4}^{2})$ $\displaystyle=$
$\displaystyle 1+2\sum_{1\leq i<j\leq
4}r_{ij}^{2}+4(r_{12}^{2}r_{34}^{2}+r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})+$
$\displaystyle
8(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})+$
$\displaystyle
16(r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31}+r_{13}r_{32}r_{24}r_{41}).$

The proof is completed. $\square$

Lemma 6.7 studies $E(X_{1}^{2}X_{2}^{2}\cdots X_{n}^{2})$ for $n=4$, which is
sufficient for our purpose. It is interesting to see the formula for $n\geq 5$
by pure curiosity. If we argue the same way as in the proof of Lemma 6.7
through Lemma 6.4, the sorting procedure would be messy. However, Lemma 6.5
provides a way to do so by figuring out the non-negative integer solutions of
the system equations $2\beta_{ii}+\sum_{j:j\neq i}\beta_{ij}=2$ for $1\leq
i\leq n.$

The following is a key step to study the covariance between two squared sample
covariances stated in Lemma 6.9.

######  LEMMA 6.8

Let $\\{(X_{1j},X_{2j},X_{3j},X_{4j})^{T}\in\mathbb{R}^{4};\,1\leq j\leq m\\}$
be i.i.d. random vectors with distribution $N_{4}(\mathbb{0},\mathbb{R})$,
where $\mathbb{R}=(r_{ij})_{4\times 4}$ and $r_{ii}=1$ for each $i$. Set

$\displaystyle B_{1}=\frac{1}{m}\sum_{j=1}^{m}X_{1j}X_{2j}\ \ \ \ \mbox{and}\
\ \ \ B_{2}=\frac{1}{m}\sum_{j=1}^{m}X_{3j}X_{4j}.$

Then,

(i) $\mbox{Var}(B_{1})=\frac{1}{m}(1+r_{12}^{2})$;

(ii) $\mbox{Cov}(B_{1},B_{2})=\frac{1}{m}(r_{13}r_{24}+r_{14}r_{23})$;

(iii)
$E[(B_{1}-r_{12})(B_{2}-r_{34})^{2}]=\frac{2}{m^{2}}(r_{13}r_{23}+r_{14}r_{24}+r_{13}r_{34}r_{42}+r_{14}r_{43}r_{32})$;

(iv) $E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]$ is identical to

$\displaystyle\ \
\frac{1}{m^{2}}(1+r_{12}^{2})(1+r_{34}^{2})+\frac{2}{m^{2}}\big{[}r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2}+2r_{13}r_{32}r_{24}r_{41}\big{]}$
$\displaystyle+\frac{2}{m^{3}}\big{[}(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+2(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+2r_{12}r_{23}r_{34}r_{41}+2r_{12}r_{24}r_{43}r_{31}+6r_{13}r_{32}r_{24}r_{41}\big{]}.$

Comment. We now conduct an independent check of the accuracy of Lemma 6.8(iv)
for two special cases. In the above result, take
$X_{1j}=X_{2j}=X_{3j}=X_{4j}$, then $r_{ij}=1$ for all $i,j$. Consequently, we
get from Theorem 6.8(iv) that

$\displaystyle
E\Big{[}\frac{1}{m}\sum_{j=1}^{m}(X_{1j}^{2}-1)\Big{]}^{4}=E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]=\frac{12}{m^{2}}+\frac{48}{m^{3}}.$

Notice $\\{X_{1j};\,1\leq j\leq m\\}$ are i.i.d. $N(0,1)$. Then
$\sum_{j=1}^{m}(X_{1j}^{2}-1)$ is a sum of i.i.d. random variables with
$E(X_{11}^{2}-1)^{2}=2$ and $E(X_{11}^{2}-1)^{4}=60$. By a classical formula
(see, for example, p. 69 from Durrett (2019)),

$\displaystyle E\Big{[}\frac{1}{m}\sum_{j=1}^{m}(X_{1j}^{2}-1)\Big{]}^{4}$
$\displaystyle=$ $\displaystyle\frac{1}{m^{4}}\big{\\{}m\cdot
E(X_{11}^{2}-1)^{4}+3m(m-1)\cdot\big{[}E(X_{11}^{2}-1)^{2}\big{]}^{2}\big{\\}}$
$\displaystyle=$ $\displaystyle\frac{12}{m^{2}}+\frac{48}{m^{3}}.$

So Theorem 6.8(iv) recovers the case for $r_{ij}=1$ for all $i,j$. On the
other hand, assume the two $2$-dimensional random vectors
$(X_{11}\,X_{21})^{T}$ and $(X_{31}\,X_{41})^{T}$ are independent, that is,
$r_{13}=r_{14}=r_{23}=r_{24}=0$. By Theorem 6.8(iv),

$\displaystyle
E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]=\frac{1}{m^{2}}(1+r_{12}^{2})(1+r_{34}^{2}).$

On the other hand, by independence and Theorem 6.8(i),

$\displaystyle
E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]=E[(B_{1}-r_{12})^{2}]\cdot
E[(B_{2}-r_{34})^{2}]=\frac{1}{m^{2}}(1+r_{12}^{2})(1+r_{34}^{2}).$

So this independent check indicates that Theorem 6.8(iv) holds for the case
$r_{13}=r_{14}=r_{23}=r_{24}=0$.

Proof of Lemma 6.8. First, by Lemma 6.4,

$\displaystyle E(X_{11}^{2}X_{21}^{2})=1+2r_{12}^{2},\ \
E(X_{11}^{3}X_{21})=3r_{12},$ (31) $\displaystyle
E(X_{11}^{2}X_{21}X_{31})=r_{23}+2r_{12}r_{13},\ \
E(X_{11}X_{21}X_{31}X_{41})=r_{12}r_{34}+r_{13}r_{24}+r_{14}r_{23}.\ \ \ \ \ \
\ \ $ (32)

In fact, the two middle identities from the four in the above can be deduced
immediately from the last one by taking $X_{11}=X_{21}=X_{31}$ and
$X_{11}=X_{21}$, respectively.

(i) Notice $mB_{1}$ and $mB_{2}$ are sums of i.i.d. random variables with mean
$r_{12}$ and $r_{34}$, respectively. Also, $X_{ij}\sim N(0,1)$ for each $i,j.$
Then $E[(B_{1}-r_{12})^{2}]=\frac{1}{m}\mbox{Var}(X_{11}X_{21})$. Since
$\mbox{Cov}(X_{11},X_{21})=E(X_{11}X_{21})=r_{12}$, we have from (31) that
$\mbox{Var}(X_{11}X_{21})=1+r_{12}^{2}$. So (i) follows.

(ii) By independence,
$E[(B_{1}-r_{12})(B_{2}-r_{34})]=\frac{1}{m}\mbox{Cov}(X_{11}X_{21},X_{31}X_{41}).$
From (32),

$\displaystyle\mbox{Cov}(X_{11}X_{21},X_{31}X_{41})=r_{12}r_{34}+r_{13}r_{24}+r_{14}r_{23}-r_{12}r_{34}=r_{13}r_{24}+r_{14}r_{23}.$
(33)

(iii) Write

$\displaystyle(B_{2}-r_{34})^{2}=\frac{1}{m^{2}}\sum_{j=1}^{m}(X_{3j}X_{4j}-r_{34})^{2}+\frac{2}{m^{2}}\sum_{1\leq
k<l\leq m}(X_{3k}X_{4k}-r_{34})(X_{3l}X_{4l}-r_{34}).$ (34)

By independence, the covariance between $B_{1}-r_{12}$ and any term from the
last sum is zero. This implies

$\displaystyle E[(B_{1}-r_{12})(B_{2}-r_{34})^{2}]$ $\displaystyle=$
$\displaystyle\frac{1}{m^{2}}E[(X_{11}X_{21}-r_{12})(X_{31}X_{41}-r_{34})^{2}].$
(35)

Use expansion
$(X_{31}X_{41}-r_{34})^{2}=X_{31}^{2}X_{41}^{2}-2r_{34}X_{31}X_{41}+r_{34}^{2}$
to see that the last expectation in (35) is equal to
$E[(X_{11}X_{21}-r_{12})(X_{31}^{2}X_{41}^{2}-2r_{34}X_{31}X_{41})]$, which is
again equal to

$\displaystyle
E(X_{11}X_{21}X_{31}^{2}X_{41}^{2})-2r_{34}E(X_{11}X_{21}X_{31}X_{41})-r_{12}E(X_{31}^{2}X_{41}^{2})+2r_{12}r_{34}E(X_{31}X_{41})$
(36) $\displaystyle=$ $\displaystyle
E(X_{11}X_{21}X_{31}^{2}X_{41}^{2})-2r_{34}(r_{12}r_{34}+r_{13}r_{24}+r_{14}r_{23})-r_{12}(1+2r_{34}^{2})+2r_{12}r_{34}^{2}$
$\displaystyle=$ $\displaystyle
E(X_{11}X_{21}X_{31}^{2}X_{41}^{2})-2r_{34}(r_{12}r_{34}+r_{13}r_{24}+r_{14}r_{23})-r_{12}$

by (31) and (32). From Lemma 6.6(i), we see that

$\displaystyle
E(X_{11}X_{21}X_{31}^{2}X_{41}^{2})=r_{12}+2r_{12}r_{34}^{2}+2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}.$
(37)

Plug this into the previous display,we arrive at

$\displaystyle
E[(B_{1}-r_{12})(B_{2}-r_{34})^{2}]=\frac{1}{m^{2}}\big{(}2r_{13}r_{23}+2r_{14}r_{24}+2r_{13}r_{24}r_{34}+2r_{14}r_{23}r_{34}\big{)}.$

(iv) By (i),

$\displaystyle E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]$ (38) $\displaystyle=$
$\displaystyle\,\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}+E(B_{1}-r_{12})^{2}\cdot
E(B_{2}-r_{34})^{2}$ $\displaystyle=$
$\displaystyle\,\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}+\frac{1}{m^{2}}(1+r_{12}^{2})(1+r_{34}^{2}).$

Similar to (34), we have

$\displaystyle(B_{1}-r_{12})^{2}=\frac{1}{m^{2}}\sum_{j=1}^{m}(X_{1j}X_{2j}-r_{12})^{2}+\frac{2}{m^{2}}\sum_{1\leq
k<l\leq m}(X_{1k}X_{2k}-r_{12})(X_{1l}X_{2l}-r_{12}).$ (39)

Recall the $m$ random variables $\\{(X_{1j},X_{2j},X_{3j},X_{4j})^{T};\,1\leq
j\leq m\\}$ are i.i.d., thus each term from the double sums in (34) and (39)
is a product of two independent random variables. This implies that
$E[(X_{3k}X_{4k}-r_{34})(X_{3l}X_{4l}-r_{34})]=0$ for any $k<l$, and the term
$(X_{1j}X_{2j}-r_{12})^{2}$ is uncorrelated to
$(X_{3k}X_{4k}-r_{34})(X_{3l}X_{4l}-r_{34})$ for any $j,k,l$ with $k<l.$ By
the same spirit, it is easy to check that
$(X_{1k}X_{2k}-r_{12})(X_{1l}X_{2l}-r_{12})$ and
$(X_{3k_{1}}X_{4k_{1}}-r_{34})(X_{3l_{1}}X_{4l_{1}}-r_{34})$ are uncorrelated
for any $k<l$ and $k_{1}<l_{1}$ as long as $(k,l)\neq(k_{1},l_{1})$. These
yield

$\displaystyle\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}$
(40) $\displaystyle=$
$\displaystyle\frac{1}{m^{3}}\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})^{2},(X_{31}X_{41}-r_{34})^{2}\big{)}+$
$\displaystyle\frac{4}{m^{4}}\sum_{1\leq k<l\leq
m}\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})(X_{12}X_{22}-r_{12}),(X_{31}X_{41}-r_{34})(X_{32}X_{42}-r_{34})\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{m^{3}}\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})^{2},(X_{31}X_{41}-r_{34})^{2}\big{)}+$
$\displaystyle\frac{2(m-1)}{m^{3}}\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})(X_{12}X_{22}-r_{12}),(X_{31}X_{41}-r_{34})(X_{32}X_{42}-r_{34})\big{)}.$

For brevity of notation, let $(X_{1},X_{2},X_{3},X_{4})\in\mathbb{R}^{4}$ and
$(Y_{1},Y_{2},Y_{3},Y_{4})\in\mathbb{R}^{4}$ be i.i.d. random vectors with
distribution $N_{4}(\mathbb{0},\mathbb{R}).$ Then the last covariance in (40)
is identical to

$\displaystyle
E\big{[}(X_{1}X_{2}-r_{12})(Y_{1}Y_{2}-r_{12})(X_{3}X_{4}-r_{34})(Y_{3}Y_{4}-r_{34})\big{]}$
(41) $\displaystyle=$
$\displaystyle\big{\\{}E[(X_{1}X_{2}-r_{12})(X_{3}X_{4}-r_{34})]\big{\\}}^{2}$
$\displaystyle=$ $\displaystyle
r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2}+2r_{13}r_{32}r_{24}r_{41}$

by independence and (33). Now we calculate the first covariance in (40). In
fact,

$\displaystyle\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})^{2},(X_{31}X_{41}-r_{34})^{2}\big{)}$
(42) $\displaystyle=$
$\displaystyle\mbox{Cov}\big{(}(X_{1}X_{2}-r_{12})^{2},(X_{3}X_{4}-r_{34})^{2}\big{)}$
$\displaystyle=$
$\displaystyle\mbox{Cov}\big{(}X_{1}^{2}X_{2}^{2},X_{3}^{2}X_{4}^{2})-2r_{12}\mbox{Cov}\big{(}X_{1}X_{2},X_{3}^{2}X_{4}^{2})-2r_{34}\mbox{Cov}\big{(}X_{1}^{2}X_{2}^{2},X_{3}X_{4})+$
$\displaystyle 4r_{12}r_{34}\mbox{Cov}\big{(}X_{1}X_{2},X_{3}X_{4}\big{)}$

since
$(X_{1}X_{2}-r_{12})^{2}=X_{1}^{2}X_{2}^{2}-2r_{12}X_{1}X_{2}+r_{12}^{2}$ and
$(X_{3}X_{4}-r_{34})^{2}=X_{3}^{2}X_{4}^{2}-2r_{34}X_{3}X_{4}+r_{34}^{2}$.
From (31) and Lemma 6.7,

$\displaystyle\mbox{Cov}(X_{1}^{2}X_{2}^{2},X_{3}^{2}X_{4}^{2})$
$\displaystyle=$ $\displaystyle
E(X_{1}^{2}X_{2}^{2}X_{3}^{2}X_{4}^{2})-E(X_{1}^{2}X_{2}^{2})\cdot
E(X_{3}^{2}X_{4}^{2})$ $\displaystyle=$ $\displaystyle
2(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+4(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+8(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+16(r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31}+r_{13}r_{32}r_{24}r_{41}).$

Also, by (37),

$\displaystyle\mbox{Cov}\,(X_{1}X_{2},X_{3}^{2}X_{4}^{2})$ $\displaystyle=$
$\displaystyle E(X_{1}X_{2}X_{3}^{3}X_{4}^{2})-r_{12}(1+2r_{34}^{2})$
$\displaystyle=$ $\displaystyle
2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}.$

By exchanging “$1$” to “$3$” and exchanging “$2$” and “$4$” in the above, we
obtain

$\displaystyle\mbox{Cov}\,(X_{1}^{2}X_{2}^{2},X_{3}X_{4})=2r_{13}r_{14}+2r_{23}r_{24}+4r_{13}r_{24}r_{12}+4r_{23}r_{14}r_{12}.$

Combining the above computations with (33) and (42), we arrive at

$\displaystyle\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})^{2},(X_{31}X_{41}-r_{34})^{2}\big{)}$
$\displaystyle=$ $\displaystyle
2(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+4(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+8(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+16(r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31}+r_{13}r_{32}r_{24}r_{41})$
$\displaystyle-2r_{12}\big{(}2r_{13}r_{23}+2r_{14}r_{24}+4r_{13}r_{24}r_{34}+4r_{14}r_{23}r_{34}\big{)}$
$\displaystyle-2r_{34}\big{(}2r_{13}r_{14}+2r_{23}r_{24}+4r_{13}r_{24}r_{12}+4r_{23}r_{14}r_{12}\big{)}$
$\displaystyle+4r_{12}r_{34}(r_{13}r_{24}+r_{14}r_{23}).$

A careful cancellation leads to

$\displaystyle\mbox{Cov}\big{(}(X_{11}X_{21}-r_{12})^{2},(X_{31}X_{41}-r_{34})^{2}\big{)}$
$\displaystyle=$ $\displaystyle
2(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+4(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+4(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+4r_{12}r_{23}r_{34}r_{41}+4r_{12}r_{24}r_{43}r_{31}+16r_{13}r_{32}r_{24}r_{41}.$

By rewriting $\frac{2(m-1)}{m^{3}}=\frac{2}{m^{2}}-\frac{2}{m^{3}}$, we see
from the above, (40) and (41) that

$\displaystyle\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{m^{3}}\big{[}2(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+4(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+4(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+4r_{12}r_{23}r_{34}r_{41}+4r_{12}r_{24}r_{43}r_{31}+16r_{13}r_{32}r_{24}r_{41}\big{]}$
$\displaystyle+\Big{(}\frac{2}{m^{2}}-\frac{2}{m^{3}}\Big{)}(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2}+2r_{13}r_{32}r_{24}r_{41}),$

which is again equal to

$\displaystyle\frac{2}{m^{2}}\big{[}r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2}+2r_{13}r_{32}r_{24}r_{41}\big{]}$
$\displaystyle+$
$\displaystyle\frac{1}{m^{3}}\big{[}2(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+2(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle+$ $\displaystyle
4(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle+$ $\displaystyle
4r_{12}r_{23}r_{34}r_{41}+4r_{12}r_{24}r_{43}r_{31}+12r_{13}r_{32}r_{24}r_{41}\big{]}.$

This and (38) imply the desired conclusion. $\square$

Now we compute the covariance of two squared sample covariances.

######  LEMMA 6.9

Let $\\{(X_{1j},X_{2j},X_{3j},X_{4j})^{T}\in\mathbb{R}^{4};\,1\leq j\leq m\\}$
be i.i.d. random vectors with distribution $N_{4}(\mathbb{0},\mathbb{R})$,
where $\mathbb{R}=(r_{ij})_{4\times 4}$ and $r_{ii}=1$ for each $i$. Let
$B_{1}$ and $B_{2}$ be defined as in Lemma 6.8. Then

$\displaystyle\mbox{Cov}(B_{1}^{2},B_{2}^{2})=\frac{1}{m}(r_{12}r_{23}r_{34}r_{41}+r_{12}r_{24}r_{43}r_{31})+\frac{\delta_{m}}{m^{2}}\sum_{1\leq
i<j\leq 4}r_{ij}^{2},$

where $|\delta_{m}|\leq\kappa$ and $\kappa$ is a numerical constant not
depending on $m$ or $\mathbb{R}.$

Proof of Lemma 6.9. Write

$\displaystyle B_{1}^{2}=(B_{1}-r_{12})^{2}+2r_{12}(B_{1}-r_{12})+r_{12}^{2}\
\ \mbox{and}\ \
B_{2}^{2}=(B_{2}-r_{34})^{2}+2r_{34}(B_{2}-r_{34})+r_{34}^{2}.$

It follows that

$\displaystyle\mbox{Cov}(B_{1}^{2},B_{2}^{2})$ $\displaystyle=$
$\displaystyle\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}+2r_{12}\cdot\mbox{Cov}\big{(}(B_{1}-r_{12}),(B_{2}-r_{34})^{2}\big{)}$
$\displaystyle+2r_{34}\cdot\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})\big{)}+r_{12}r_{34}\cdot\mbox{Cov}\big{(}B_{1},B_{2}).$

By Lemma 6.8(i) & (iv),

$\displaystyle\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}$
$\displaystyle=$ $\displaystyle
E[(B_{1}-r_{12})^{2}(B_{2}-r_{34})^{2}]-E[(B_{1}-r_{12})^{2}]\cdot
E[(B_{2}-r_{34})^{2}]$ $\displaystyle=$
$\displaystyle\frac{2}{m^{2}}\big{[}r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2}+2r_{13}r_{32}r_{24}r_{41}\big{]}+\frac{2}{m^{3}}\big{[}(r_{13}^{2}+r_{14}^{2}+r_{23}^{2}+r_{24}^{2})+(r_{13}^{2}r_{24}^{2}+r_{14}^{2}r_{23}^{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+2(r_{12}r_{23}r_{31}+r_{12}r_{24}r_{41}+r_{23}r_{34}r_{42}+r_{13}r_{34}r_{41})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+2r_{12}r_{23}r_{34}r_{41}+2r_{12}r_{24}r_{43}r_{31}+6r_{13}r_{32}r_{24}r_{41}\big{]}.$

Use $x^{2}y^{2}\leq x^{2}$ and $|xyzu|\leq|xyz|\leq|xy|\leq x^{2}+y^{2}$ for
any $x,y,z,u\in[-1,1]$ to see that

$\displaystyle\big{|}\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})^{2}\big{)}\big{|}\leq\frac{K_{1}}{m^{2}}\sum_{1\leq
i<j\leq 4}r_{ij}^{2},$

where $K_{1}$ is a numerical constant not depending on $m$ or $\mathbb{R}.$
Same bounds hold with another numerical constant $K_{2}$ for
$2r_{12}\cdot\mbox{Cov}\big{(}(B_{1}-r_{12}),(B_{2}-r_{34})^{2}\big{)}$ and
$2r_{34}\cdot\mbox{Cov}\big{(}(B_{1}-r_{12})^{2},(B_{2}-r_{34})\big{)}$ by
Lemma 6.8(iii). Finally, by Lemma 6.8(ii),

$\displaystyle
r_{12}r_{34}\cdot\mbox{Cov}\big{(}B_{1},B_{2})=\frac{1}{m}(r_{12}r_{24}r_{43}r_{31}+r_{12}r_{23}r_{34}r_{41}).$

The proof is finished by combining all of the estimates. $\square$

Review sample correlation coefficient $\hat{r}_{ij}$ defined in Theorem 4, the
denominator is not easy to be handled. By the heuristic of the law of large
numbers and the Taylor expansion, we are able to express $\hat{r}_{ij}^{2}$ as
a polynomial of Gaussian random variables plus an error. The following result
provides basic computations for this step. It will be used in the derivation
of Lemma 6.12.

######  LEMMA 6.10

Let $\\{(X_{i},Y_{i})^{T};\,1\leq i\leq m\\}$ be i.i.d. $2$-dimensional normal
random vectors with $EX_{1}=EY_{1}=0$, $EX_{1}^{2}=EY_{1}^{2}=1$ and
$\mbox{Cov}(X_{1},Y_{1})=r$. Define

$\displaystyle A=\frac{X_{1}^{2}+\cdots+X_{m}^{2}}{m}\ \ \ \ \mbox{and}\ \ \ \
B=\frac{Y_{1}^{2}+\cdots+Y_{m}^{2}}{m}.$

The following statements hold for all $m\geq 2$ with $|\delta_{m}|<\kappa$
where $\kappa$ is a constant not depending on $m$ or $r$.

(i) $E[(1-A)(1-B)]=\frac{2}{m}r^{2}$.

(ii) $E[(1-A)(1-B)^{2}]=-\frac{8}{m^{2}}r^{2}$.

(iii) $E[(1-A)(1-B)^{3}]=\frac{12r^{2}}{m^{2}}+\frac{\delta_{m}}{m^{3}}$.

(iv) $E[(1-A)^{2}(1-B)^{2}]=\frac{8r^{4}+4}{m^{2}}+\frac{\delta_{m}}{m^{3}}$.

(v) $E[(1-A)^{i}(1-B)^{j}]=\frac{\delta_{m}}{m^{3}}$ for
$(i,j)=(0,5),(1,4),(2,3).$

Proof of Lemma 6.10. Write

$\displaystyle 1-A=\frac{1}{m}\sum_{i=1}^{m}(1-X_{i}^{2})\ \ \ \mbox{and}\ \ \
1-B=\frac{1}{m}\sum_{i=1}^{m}(1-Y_{i}^{2}).$ (43)

(i) Easily,

$\displaystyle\mbox{Var}(X_{1}Y_{1})=E\big{[}(X_{1}Y_{1}-r)^{2}\big{]}=E\big{(}X_{1}^{2}Y_{1}^{2}\big{)}-r^{2}=r^{2}+1;$
(44) $\displaystyle
E[(1-X_{1}^{2})(1-Y_{1}^{2})]=1-E(X_{1}^{2})-E(Y_{1}^{2})+E(X_{1}^{2}Y_{1}^{2})=2r^{2}$
(45)

by using (31). Due to independence,

$\displaystyle
E[(1-A)(1-B)]=\frac{1}{m}E\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})\big{]}=2r^{2}.$

(ii) Write

$\displaystyle
m^{2}\cdot(1-B)^{2}=\sum_{i=1}^{m}(1-Y_{i}^{2})^{2}+2\sum_{1\leq i<j\leq
m}(1-Y_{i}^{2})(1-Y_{j}^{2}).$ (46)

Evidently, $1-X_{k}^{2}$ is uncorrelated to $(1-Y_{i}^{2})(1-Y_{j}^{2})$ for
any $i,j,k$ with $i\neq j.$ Consequently,

$\displaystyle
E[(1-A)(1-B)^{2}]=\frac{1}{m^{2}}E\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})^{2}\big{]}.$
(47)

Use the fact $(1-Y_{1}^{2})^{2}=1-2Y_{1}^{2}+Y_{1}^{4}$ to see

$\displaystyle E\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})^{2}\big{]}$ $\displaystyle=$
$\displaystyle-2E\big{[}(1-X_{1}^{2})Y_{1}^{2}\big{]}+E\big{[}(1-X_{1}^{2})Y_{1}^{4}\big{]}$
(48) $\displaystyle=$ $\displaystyle
2E\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})\big{]}+E\big{(}Y_{1}^{4}\big{)}-E\big{(}X_{1}^{2}Y_{1}^{4}\big{)}.$

Write $X_{1}=rY_{1}+\sqrt{1-r^{2}}Y_{1}^{\prime}$, where $Y_{1}^{\prime}\sim
N(0,1)$ and $Y_{1}^{\prime}$ is independent of $X_{1}$. Then

$\displaystyle
E(X_{1}^{2}Y_{1}^{4})=r^{2}E(Y_{1}^{6})+(1-r^{2})E(Y_{1}^{\prime
2})E(Y_{1}^{4})=15r^{2}+3(1-r^{2})=3+12r^{2}.$

Notice $E(Y_{1}^{4})=3$ and $E(Y_{1}^{6})=15$. By (45) and (48),

$\displaystyle E\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})^{2}\big{]}=2\cdot
2r^{2}+3-(3+12r^{2})=-8r^{2}.$

Then (ii) follows from (47).

(iii) Write

$\displaystyle m^{3}(1-B)^{3}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}(1-Y_{i}^{2})^{3}+3\sum_{1\leq i\neq j\leq
m}(1-Y_{i}^{2})^{2}(1-Y_{j}^{2})$ (49) $\displaystyle+\sum_{1\leq i\neq j\neq
k\leq m}(1-Y_{i}^{2})(1-Y_{j}^{2})(1-Y_{k}^{2})$ $\displaystyle:=$
$\displaystyle U_{1}+U_{2}+U_{3}.$

Notice $1-X_{k}^{2}$ is uncorrelated to $U_{3}$ and $1-X_{i}^{2}$ is
uncorrelated to $(1-Y_{i}^{2})^{2}(1-Y_{j}^{2})$ if $k=i\neq j$. This gives

$\displaystyle m^{4}\cdot E[(1-A)(1-B)^{3}]$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}E\big{[}(1-X_{i}^{2})(1-Y_{i}^{2})^{3}\big{]}+3\sum_{1\leq
i\neq j\leq m}E\big{[}(1-X_{j}^{2})(1-Y_{i}^{2})^{2}(1-Y_{j}^{2})\big{]}$
$\displaystyle=$ $\displaystyle
mE\big{[}(1-X_{1}^{2})(1-Y_{1}^{2})^{3}\big{]}+3m(m-1)E\big{[}(1-Y_{1}^{2})^{2}\big{]}\cdot
E\big{[}(1-X_{2}^{2})(1-Y_{2}^{2})\big{]}.$

The product of the last two expectations is equal to $2\cdot(2r^{2})=4r^{2}$
by (45). This implies

$\displaystyle
E[(1-A)(1-B)^{3}]=\frac{12r^{2}}{m^{2}}+\frac{\delta_{m}}{m^{3}}$

with $|\delta_{m}|<\kappa$ for all $m\geq 2$, where $\kappa$ here and later
represents a constant not depending on $m$ or $r$, and can be different from
line to line.

(iv) Similar to (46),

$\displaystyle m^{2}\cdot(1-A)^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}(1-X_{i}^{2})^{2}+2\sum_{1\leq i<j\leq
m}(1-X_{i}^{2})(1-X_{j}^{2})$ (50) $\displaystyle:=$ $\displaystyle
V_{1}+V_{2}.$

Trivially, $(1-X_{i}^{2})^{2}$ is uncorrelated to $(1-Y_{j}^{2})(1-Y_{k}^{2})$
for any $i,j,k$ with $j\neq k$. The same is true if $X$ and $Y$ are switched.
Thus

$\displaystyle m^{4}\cdot\mbox{Cov}\big{(}(1-A)^{2},(1-B)^{2}\big{)}$
$\displaystyle=$ $\displaystyle
m\,\mbox{Cov}\big{(}(1-X_{1}^{2})^{2},(1-Y_{1}^{2})^{2}\big{)}+2m(m-1)\mbox{Cov}\big{(}(1-X_{1}^{2})(1-X_{2}^{2}),(1-Y_{1}^{2})(1-Y_{2}^{2})\big{)}.$

Easily, by independence and (45),

$\displaystyle\mbox{Cov}\big{(}(1-X_{1}^{2})(1-X_{2}^{2}),(1-Y_{1}^{2})(1-Y_{2}^{2})\big{)}=\big{[}E(1-X_{1}^{2})(1-Y_{1}^{2})\big{]}^{2}=4r^{4}.$

Evidently, $E[(1-A)^{2}]=\frac{2}{m}$. It follows that

$\displaystyle E\big{[}(1-A)^{2}(1-B)^{2}\big{]}$ $\displaystyle=$
$\displaystyle\mbox{Cov}\big{(}(1-A)^{2},(1-B)^{2}\big{)}+E\big{[}(1-A)^{2}\big{]}\cdot
E\big{[}(1-B)^{2}\big{]}$ $\displaystyle=$
$\displaystyle\frac{1}{m^{4}}\big{[}m\,\mbox{Cov}\big{(}(1-X_{1}^{2})^{2},(1-Y_{1}^{2})^{2}\big{)}+8m(m-1)r^{4}\big{]}+\frac{4}{m^{2}}$
$\displaystyle=$
$\displaystyle\frac{8r^{4}+4}{m^{2}}+\frac{\delta_{m}}{m^{3}}.$

(v) We will study the three cases one by one.

Case 1: $(i,j)=(0,5)$. Let $Z_{1},\cdots,Z_{m}$ be i.i.d. random variables
with mean zero. Write

$\displaystyle(Z_{1}+\cdots+Z_{m})^{5}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}Z_{i}^{5}+c_{1}\sum_{i\neq
j}Z_{i}^{4}Z_{j}+c_{2}\sum_{i\neq j}Z_{i}^{3}Z_{j}^{2}+c_{3}\sum_{i\neq j\neq
k}Z_{i}^{3}Z_{j}Z_{k}$ $\displaystyle+c_{4}\sum_{i\neq j\neq
k}Z_{i}^{2}Z_{j}^{2}Z_{k}+c_{5}\sum_{i\neq j\neq
k}Z_{i}^{2}Z_{j}Z_{k}Z_{l}+c_{6}\sum_{a\neq i\neq j\neq k\neq
l}Z_{a}Z_{i}Z_{j}Z_{k}Z_{l},$

where $c_{1},\cdots,c_{6}$ are numerical coefficients not depending on $m$ or
$r$. By independence, only the first and third sums of the right hand side
above are non-zero if the expectation is taken on both sides. The total number
of the terms appeared in the two sums is $m+2\cdot\binom{m}{2}=m^{2}$. Take
$Z_{i}=1-Y_{i}^{2}$ for each $i$ to have

$\displaystyle E[(1-B)^{5}]=\frac{\delta_{m}}{m^{3}}$

for all $m\geq 2$.

Case 2: $(i,j)=(1,4)$. Similar to Case 1,

$\displaystyle(Z_{1}+\cdots+Z_{m})^{4}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}Z_{i}^{4}+c_{1}\sum_{i\neq
j}Z_{i}^{3}Z_{j}+6\sum_{i<j}Z_{i}^{2}Z_{j}^{2}+c_{3}\sum_{i\neq j\neq
k}Z_{i}^{2}Z_{j}Z_{k}$ (51) $\displaystyle+c_{4}\sum_{i\neq j\neq k\neq
l}Z_{i}Z_{j}Z_{k}Z_{l}.$

Take $Z_{i}=1-Y_{i}^{2}$ for each $i$. Recall (43), each term $1-X_{i}^{2}$ is
uncorrelated to any term in the last two sums. By the same argument as in Case
1, we know $E[(1-A)(1-B)^{4}]=\frac{1}{m^{3}}\delta_{m}$ for all $m\geq 2.$

Case 3: $(i,j)=(2,3)$. Review (49) and (50). Easily, $U_{3}$ is uncorrelated
to both $V_{1}$ and $V_{2}$; $U_{1}$ is uncorrelated to $V_{2}$. It follows
that

$\displaystyle
m^{5}\cdot\mbox{Cov}\big{(}(1-A)^{2},(1-B)^{3}\big{)}=\mbox{Cov}\big{(}U_{1},V_{1}\big{)}+\mbox{Cov}\big{(}U_{2},V_{1}\big{)}+\mbox{Cov}\big{(}U_{2},V_{2}\big{)}.$

By independence and the same argument as before, it is easy to see that
$\mbox{Cov}\big{(}U_{1},V_{1}\big{)}=m\delta_{m}$,
$\mbox{Cov}\big{(}U_{2},V_{1}\big{)}=m^{2}\delta_{m}^{\prime}$ and
$\mbox{Cov}\big{(}U_{2},V_{2}\big{)}=m^{2}\delta_{m}^{\prime\prime}$, where
$|\delta_{m}|+|\delta_{m}^{\prime}|+|\delta_{m}^{\prime\prime}|\leq\kappa$ and
$\kappa$ is a constant not depending on $m$ or $r.$ Consequently,

$\displaystyle\mbox{Cov}\big{(}(1-A)^{2},(1-B)^{3}\big{)}=\frac{\delta_{m}}{m^{3}}.$
(52)

From (50) and (49), it is readily seen $E[(1-A)^{2}]=\frac{2}{m}$ and
$E[(1-B)^{3}]=\frac{1}{m^{2}}\cdot E(1-Y_{1}^{2})^{3}$. The conclusion then
follows from these facts, (52) and the formula

$\displaystyle
E\big{[}(1-A)^{2}(1-B)^{3}\big{]}=\mbox{Cov}\big{(}(1-A)^{2},(1-B)^{3}\big{)}+E\big{[}(1-A)^{2}\big{]}\cdot
E\big{[}(1-B)^{3}].$

The proof is finished. $\square$

####  6.2.2 Evaluations of Covariances between Monomials of Gaussian Random
Variables

Let $(X_{1},X_{2},X_{3},X_{4})^{T}$ be a $4$-dimensional random vector with
distribution $N_{4}(\mathbb{0},\mathbb{R})$ where
$\mathbb{R}=(r_{ij})_{4\times 4}$ and $r_{ii}=1$ for each $i$. One of major
tasks in this section is showing
$|\mbox{Cov}(X_{1}^{\alpha}X_{2}^{\beta},X_{3}^{\gamma}X_{4}^{\delta})|\leq
C\sum_{1\leq i<j\leq 4}r_{ij}^{2}$ for some constant $C$ depending only on
non-negative integers $\alpha,\beta,\gamma,\delta$. Review the $t$-norm
$\|U\|_{t}=[E(|U|^{t})]^{1/t}$ for any random variable $U$ and $t\geq 1$. In
the following lemma, by convention we regard $\prod_{i=1}^{0}a_{i}=1$ for any
$a_{i}$.

######  LEMMA 6.11

Let $\\{X_{1},\cdots,X_{k},Y_{1},\cdots,Y_{l}\\}$ be random variables with
$k\geq 0$, $l\geq 0$ and $k+l\geq 1$. Let
$p_{1},\cdots,p_{k},q_{1},\cdots,q_{l}$ be positive integers. Assume, for each
$i$, (a) $X_{i}=\frac{1}{\sqrt{p_{i}}}\sum_{j=1}^{p_{i}}\xi_{ij}$ where
$\\{\xi_{i1},\cdots,\xi_{ip_{i}}\\}$ are i.i.d. with mean $0$, and (b)
$Y_{i}=\frac{1}{q_{i}}\chi^{2}(q_{i})$. Then

$\displaystyle
E\Big{[}\prod_{i=1}^{k}|X_{i}|^{\alpha_{i}}\cdot\prod_{j=1}^{l}Y_{j}^{\beta_{j}}\Big{]}<C\cdot\prod_{i=1}^{k}\|\,|\xi_{i1}|^{\alpha_{i}}\|_{k+l}$

for any $\alpha_{i}\geq 2/(k+l)$, $\beta_{j}>(2-q_{j})/[2(k+l)]$, $1\leq i\leq
k$ and $1\leq j\leq l$, where $C$ is a constant depending on $k$, $l$,
$\alpha_{i}$ and $\beta_{j}$ but not depending on $p_{i}$ or $q_{j}$.

Proof of Lemma 6.11. First we assume $k\geq 1$ and $l\geq 1$. By applying the
Hölder inequality to the product of $k+l$ terms, we see

$\displaystyle
E\Big{[}\prod_{i=1}^{k}|X_{i}|^{\alpha_{i}}\cdot\prod_{j=1}^{l}Y_{j}^{\beta_{j}}\Big{]}\leq\prod_{i=1}^{k}\big{(}E|X_{i}|^{\alpha_{i}(k+l)}\big{)}^{1/(k+l)}\cdot\prod_{j=1}^{l}\big{(}EY_{j}^{\beta_{j}(k+l)}\big{)}^{1/(k+l)}.$
(53)

By assumption, $\alpha_{i}(k+l)\geq 2$ for each $i$. According to the
Marcinkiewicz-Zygmund inequality [see, for example, the proof of Corollary 2
on p. 387 from Chow and Teicher (1997)],

$\displaystyle E|X_{i}|^{\alpha_{i}(k+l)}\leq C(\alpha_{i},k,l)\cdot
E\big{(}|\xi_{i1}|^{\alpha_{i}(k+l)}\big{)},$ (54)

where $C(\alpha_{i},k,l)$ is a constant depending on $\alpha_{i},k,l$ only. As
a result,

$\displaystyle\prod_{i=1}^{k}\big{(}E|X_{i}|^{\alpha_{i}(k+l)}\big{)}^{1/(k+l)}\leq
C(\alpha_{1},\cdots,\alpha_{k},k,l)\cdot\prod_{i=1}^{k}\||\xi_{i1}|^{\alpha_{i}}\|_{k+l}.$
(55)

Second, by recalling the density of $\chi^{2}(m)$ is
$f(x):=\frac{1}{\Gamma(\frac{m}{2})2^{m/2}}x^{(m/2)-1}e^{-x/2}$ for $x>0$, we
have

$\displaystyle E[\chi^{2}(m)^{\beta}]$ $\displaystyle=$
$\displaystyle\frac{1}{\Gamma(\frac{m}{2})2^{m/2}}\int_{0}^{\infty}x^{(m+2\beta)/2-1}e^{-x/2}\,dx$
(56) $\displaystyle=$
$\displaystyle\frac{\Gamma(\frac{m}{2}+\beta)}{\Gamma(\frac{m}{2})2^{-\beta}}$

for any real number $\beta>-m/2$. From the fact that
$\lim_{x\to\infty}\frac{\Gamma(x+a)}{\Gamma(x)x^{a}}=1$ for any number
$a\in\mathbb{R}$, it is seen that $E[\chi^{2}(m)/m]^{\beta}\to 1$ as
$m\to\infty.$ This implies that

$\displaystyle\sup
E\Big{[}\frac{\chi^{2}(m)}{m}\Big{]}^{\beta}=C_{\beta}<\infty,$ (57)

where the supremum is taken over all $m$ satisfying $m\geq 1$ and $m\geq
2(1-\beta)$ and where $C_{\beta}$ depends on $\beta$ only. The reason we
choose $m$ such that $\frac{m}{2}+\beta\geq 1$ is because we need
$\frac{m}{2}+\beta$ to stay away from the only singular point $0$ of
$\Gamma(x)$ defined on $[0,\infty).$ It follows that

$\displaystyle
E\big{(}|Y_{j}|^{\beta_{j}(k+l)}\big{)}=E\Big{[}\Big{(}\frac{\chi^{2}(q_{j})}{q_{j}}\Big{)}^{\beta_{j}(k+l)}\Big{]}\leq
C_{\beta_{j}(k+l)}<\infty$

provided $q_{j}\geq 2[1-\beta_{j}(k+l)]$, or equivalently,
$\beta_{j}\geq(2-q_{j})/[2(k+l)]$. This asserts that

$\displaystyle\prod_{j=1}^{l}\big{(}EY_{j}^{\beta_{j}(k+l)}\big{)}^{1/(k+l)}=C(k,l,\beta_{1},\cdots,\beta_{l})<\infty$
(58)

under the assumption $\beta_{j}\geq(2-q_{j})/[2(k+l)]$ for each $1\leq j\leq
l$. This joined with (53) and (55) leads to the desired inequality.

If $k\geq 1$ and $l=0$, then (53)-(55) still hold. If $k=0$ and $l\geq 1$,
then (53) and (58) are also true. The desired statements are then derived.
$\square$

Recall sample correlation coefficient $\hat{r}_{ij}$ defined in Theorem 4 from
Section 6.2, the next result provides an estimate for the dependency between
two denominators.

######  LEMMA 6.12

Let $\\{(X_{i},Y_{i})^{T};\,1\leq i\leq m\\}$ be i.i.d. $2$-dimensional normal
random vectors with $EX_{1}=EY_{1}=0$, $EX_{1}^{2}=EY_{1}^{2}=1$ and
$\mbox{Cov}(X_{1},Y_{1})=r$. Write

$\displaystyle
E\Big{(}\frac{m}{X_{1}^{2}+\cdots+X_{m}^{2}}\cdot\frac{m}{Y_{1}^{2}+\cdots+Y_{m}^{2}}\Big{)}=1+\frac{4+2r^{2}}{m}+\frac{12+8r^{2}+8r^{4}}{m^{2}}+\frac{\delta_{m}}{m^{3}}.$
(59)

Then $|\delta_{m}|\leq\kappa$ for all $m\geq 11$, where $\kappa$ is a constant
not depending on $m$ or $r$.

Before we present the proof let us have a quick check for two cases with $r=0$
and $r=1$. Recall (56). We have $E\frac{1}{\chi^{2}(m)}=\frac{1}{m-2}$ and
$E\frac{1}{\chi^{2}(m)^{2}}=\frac{1}{(m-2)(m-4)}$ for any $m\geq 5$. If $r=0$,
then $\\{X_{i},Y_{i};\,1\leq i\leq m\\}$ are i.i.d. $N(0,1)$. Consequently the
left hand side of (59) is identical to

$\displaystyle
m^{2}\Big{[}E\frac{1}{\chi^{2}(m)}\Big{]}^{2}=\frac{m^{2}}{(m-2)^{2}}=1+\frac{4}{m}+\frac{12}{m^{2}}+\frac{32m-48}{m^{2}(m-2)^{2}},$

which corresponds to the right hand side of (59) with $r=0$. If $r=1$, then
$X_{i}=Y_{i}$ for each $i$. Thus the left hand side of (59) becomes

$\displaystyle m^{2}\cdot
E\frac{1}{\chi^{2}(m)^{2}}=\frac{m^{2}}{(m-2)(m-4)}=1+\frac{6}{m}+\frac{28}{m^{2}}+\frac{120m-224}{m^{2}(m-2)(m-4)},$

which is equal to the right hand side of (59) with $r=1$.

Proof of Lemma 6.12. Define

$\displaystyle A=\frac{X_{1}^{2}+\cdots+X_{m}^{2}}{m}\ \ \mbox{and}\ \
B=\frac{Y_{1}^{2}+\cdots+Y_{m}^{2}}{m}.$

Set $A_{1}=1-A$ and $B_{1}=1-B$. By the formula
$x^{-1}=1+(1-x)+(1-x)^{2}+(1-x)^{3}+(1-x)^{4}+(1-x)^{5}+x^{-1}(1-x)^{6}$, we
write

$\displaystyle\frac{1}{AB}=\Big{(}\frac{A_{1}^{6}}{A}+\sum_{i=0}^{5}A_{1}^{i}\Big{)}\Big{(}\frac{B_{1}^{6}}{B}+\sum_{i=0}^{5}B_{1}^{i}\Big{)}.$

We will expand $(\sum_{i=0}^{5}A_{1}^{i})(\sum_{i=0}^{5}B_{1}^{i})$ and write
it as the sum of $A_{1}^{a}B_{1}^{b}$, and furthermore break the sum into two
sums, the first of which is for the terms with $a+b\leq 5$ and the second of
which is for those with $a+b\geq 6$. Therefore,

$\displaystyle\frac{1}{AB}$ $\displaystyle=$ $\displaystyle
1+(A_{1}+B_{1})+(A_{1}^{2}+A_{1}B_{1}+A_{2}^{2})+\big{(}A_{1}^{3}+A_{1}^{2}B_{1}+A_{1}B_{1}^{2}+B_{1}^{3}\big{)}$
(60)
$\displaystyle+\big{(}A_{1}^{4}+A_{1}^{3}B_{1}+A_{1}^{2}B_{1}^{2}+A_{1}B_{1}^{3}+B_{1}^{4}\big{)}$
$\displaystyle+\big{(}A_{1}^{5}+A_{1}^{4}B_{1}+A_{1}^{3}B_{1}^{2}+A_{1}^{2}B_{1}^{3}+A_{1}B_{1}^{4}+B_{1}^{5}\big{)}+\epsilon_{m},$

where

$\displaystyle\epsilon_{m}:=\frac{A_{1}^{6}}{A}\sum_{i=0}^{5}B_{1}^{i}+\frac{B_{1}^{6}}{B}\sum_{i=0}^{5}A_{1}^{i}+\frac{A_{1}^{6}B_{1}^{6}}{AB}+p_{m}$

and where $p_{m}=\sum A_{1}^{a}B_{1}^{b}$ with the sum running over all $a\geq
1$ and $b\geq 1$ satisfying $6\leq a+b\leq 10$. Obviously, by counting the
number of the terms on the right hand side of (60), we know the total number
of the terms in the sum of $p_{m}$ is $36-(1+2+3+4+5+6)=15$. Notice

$\displaystyle A\sim\frac{\chi^{2}(m)}{m},\ \ \ B\sim\frac{\chi^{2}(m)}{m},\ \
\ \sqrt{m}A_{1}=\frac{1}{\sqrt{m}}\sum_{i=1}^{m}(1-X_{i}^{2}),\ \ \
\sqrt{m}B_{1}=\frac{1}{\sqrt{m}}\sum_{i=1}^{m}(1-Y_{i}^{2}).$

Review the notation in Lemma 6.11. We first consider the term
$E(A_{1}^{6}B_{1}^{i}A^{-1})$. If $i=0$, take $k=l=1$, $\alpha_{1}=6$,
$\beta_{1}=-1$ and $q_{1}=m$. Then $\alpha_{1}\geq 2/(k+l)=1$ and
$\beta_{1}>(2-q_{1})/[2(k+l)]=(2-m)/4$ as $m\geq 7.$ Similarly, if $i\geq 1$,
take $k=2$, $l=1$, $\alpha_{1}=6$, $\alpha_{2}=i$, $\beta_{1}=-1$ and
$q_{1}=m$. It is always true that $\alpha_{1}\wedge\alpha_{2}\geq
2/(k+l)=2/3.$ Also, $\beta_{1}>(2-q_{1})/[2(k+l)]=(2-m)/6$ if $m\geq 9.$ A
similar but easier check can be done for $A_{1}^{a}B_{1}^{b}$ with $a\geq 1$,
$b\geq 1$ and $a+b\geq 6$. It then follows from Lemma 6.11 that

$\displaystyle
E\Big{(}\frac{A_{1}^{6}}{A}B_{1}^{i}\Big{)}=\frac{\delta(m,i)}{m^{3}}\ \
\mbox{and}\ \ E(A_{1}^{a}B_{1}^{b})=\frac{\delta(m,a,b)^{\prime}}{m^{3}}$ (61)

as $m\geq 9$, where $|\delta(m,i)|+|\delta(m,a,b)^{\prime}|\leq\kappa$ for all
$1\leq i\leq 5$ and integers $a\geq 1$ and $b\geq 1$ with $a+b\geq 6$. Here
and later $\kappa$ represents a constant not depending on $m$ or $r$ and can
be different from line to line.

Now we turn to look at the term $A_{1}^{6}B_{1}^{6}A^{-1}B^{-1}$. Take
$k=l=2$, $\alpha_{1}=\alpha_{2}=6$, $\beta_{1}=\beta_{2}=-1$ and
$q_{1}=q_{2}=m$ to see that, for each $i=1,2$, we have $\alpha_{i}\geq
2/(k+l)=1/2$ and $\beta_{i}>(2-q_{i})/[2(k+l)]=(2-m)/8$ provided $m\geq 11.$
From Lemma 6.11 again we obtain

$\displaystyle E\frac{A_{1}^{6}B_{1}^{6}}{AB}=\frac{\delta_{m}}{m^{6}}$ (62)

with $|\delta(m)|\leq\kappa$ as $m\geq 11.$ Combining (61) and (62) and using
the symmetry of $A$ and $B$ as well as that of $A_{i}$ and $B_{i}$, we see
that

$\displaystyle E\epsilon_{m}=\frac{\delta_{m}}{m^{3}}$ (63)

with $|\delta_{m}|\leq\kappa$ for all $m\geq 11$.

Back to (60), let us examine the expectation of each monomial of $A_{1}$ and
$B_{1}$ on the right hand side. By independence,

$\displaystyle
E\big{(}A_{1}^{2}\big{)}=\frac{1}{m^{2}}\sum_{i=1}^{m}E\big{[}(1-X_{i}^{2})^{2}\big{]}=\frac{2}{m}$
(64)

since $\mbox{Var}(X_{1}^{2})=2.$ Recall

$\displaystyle(t_{1}+\cdots+t_{m})^{3}=\sum_{i=1}^{m}t_{i}^{3}+3\sum_{1\leq
i\neq j\leq m}t_{i}^{2}t_{j}+\sum_{1\leq i\neq j\neq k\leq m}t_{i}t_{j}t_{k}$

for any real numbers $t_{1},\cdots,t_{m}$. Take $t_{i}=1-X_{i}^{2}$ and use
independence to see

$\displaystyle
E\big{(}A_{1}^{3}\big{)}=\frac{1}{m^{2}}E\big{[}(1-X_{1}^{2})^{3}\big{]}=\frac{1}{m^{2}}\big{[}1-3E(X_{1}^{2})+3E(X_{1}^{4})-E(X_{1}^{6})\big{]}=-\frac{8}{m^{2}}$

since $E(X_{1}^{2})=1$, $E(X_{1}^{4})=3$ and $E(X_{1}^{6})=15$. Review

$\displaystyle(t_{1}+\cdots+t_{m})^{4}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}t_{i}^{4}+c_{1}\sum_{i\neq
j}t_{i}^{3}t_{j}+6\sum_{i<j}t_{i}^{2}t_{j}^{2}+c_{3}\sum_{i\neq j\neq
k}t_{i}^{2}t_{j}t_{k}$ $\displaystyle+c_{4}\sum_{i\neq j\neq k\neq
l}t_{i}t_{j}t_{k}t_{l}$

for any $t_{1},\cdots,t_{m}$. By taking $t_{i}=1-X_{i}^{2}$ and using
independence, we get

$\displaystyle E\big{(}A_{1}^{4}\big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{m^{3}}\cdot
E\big{[}(1-X_{1}^{2})^{4}\big{]}+\frac{6}{m^{4}}\cdot\binom{m}{2}\cdot
E\big{[}(1-X_{1}^{2})^{2}(1-X_{2}^{2})^{2}\big{]}$ $\displaystyle=$
$\displaystyle\frac{12}{m^{2}}+\frac{\delta_{m}}{m^{3}}$

with $|\delta_{m}|\leq\kappa$ for all $m\geq 1$ since
$E\big{[}(1-X_{1}^{2})^{2}(1-X_{2}^{2})^{2}\big{]}=[E(1-X_{1}^{2})^{2}]^{2}=4.$
Now, by (60),

$\displaystyle E\Big{(}\frac{1}{AB}-\epsilon_{m}\Big{)}$ $\displaystyle=$
$\displaystyle
1+E\big{(}A_{1}^{2}+A_{1}B_{1}+B_{1}^{2}\big{)}+E\big{(}A_{1}^{3}+A_{1}^{2}B_{1}+A_{1}B_{1}^{2}+B_{1}^{3}\big{)}$
$\displaystyle+E\big{(}A_{1}^{4}+A_{1}^{3}B_{1}+A_{1}^{2}B_{1}^{2}+A_{1}B_{1}^{3}+B_{1}^{4}\big{)}$
$\displaystyle+E\big{(}A_{1}^{5}+A_{1}^{4}B_{1}+A_{1}^{3}B_{1}^{2}+A_{1}^{2}B_{1}^{3}+A_{1}B_{1}^{4}+B_{1}^{5}\big{)}.$

From Lemma 6.10, the expressions of $EA_{1}^{i}$ for $i=2,3,4$ above as well
as the symmetry of $A$ and $B$, we obtain

$\displaystyle E\Big{(}\frac{1}{AB}-\epsilon_{m}\Big{)}$ $\displaystyle=$
$\displaystyle
1+\Big{(}\frac{4}{m}+\frac{2}{m}r^{2}\Big{)}+\Big{(}-\frac{16}{m^{2}}-\frac{16}{m^{2}}r^{2}\Big{)}$
$\displaystyle\ \
\,+\Big{(}\frac{24}{m^{2}}+\frac{24r^{2}}{m^{2}}+\frac{8r^{4}+4}{m^{2}}\Big{)}+\frac{\delta_{m}}{m^{3}}$
$\displaystyle=$ $\displaystyle
1+\frac{4+2r^{2}}{m}+\frac{12+8r^{2}+8r^{4}}{m^{2}}+\frac{\delta_{m}}{m^{3}}$

with $|\delta_{m}|\leq\kappa$ for all $m\geq 2$. This and (63) conclude

$\displaystyle
E\frac{1}{AB}=1+\frac{4+2r^{2}}{m}+\frac{12+8r^{2}+8r^{4}}{m^{2}}+\frac{\delta_{m}}{m^{3}}$

with $|\delta_{m}|\leq\kappa$ for all $m\geq 11$. $\square$

Recall the notation $(2i-1)!!=1\cdot 3\cdots(2i-1)$ for any integer $i\geq 1$.
We set $(-1)!!=1$ by convention.

######  LEMMA 6.13

Let $X$ and $Y$ be $N(0,1)$. Assume they are jointly normal with covariance
$r$. Then the following hold.

(i) For any integers $i\geq 1$ and $j\geq 1$, there exists a polynomial $f(x)$
depending on $i$ and $j$ such that $f(0)=1$ and
$E(X^{2i-1}Y^{2j-1})=(2i-1)!!\cdot(2j-1)!!\cdot rf(r^{2}).$

(ii) For any integers $i\geq 0$ and $j\geq 0$, there exists a polynomial
$g(x)$ depending on $i$ and $j$ such that
$E(X^{2i}Y^{2j})=(2i-1)!!\cdot(2j-1)!!+r^{2}g(r^{2})$.

Proof of Lemma 6.13. Let $Z$ be a random variable with distribution $N(0,1)$
and $Z$ be independent of $X$ and $Y$. Then we may write $Y=rX+r^{\prime}Z$,
where $r^{\prime}=\sqrt{1-r^{2}}.$

(i) Easily,

$\displaystyle Y^{2j-1}=\sum_{k=0}^{2j-1}\binom{2j-1}{k}r^{k}r^{\prime
2j-k-1}X^{k}Z^{2j-k-1}.$

By independence and the fact $E[N(0,1)^{n}]=0$ for any odd number $n\geq 1$,
we have

$\displaystyle E(X^{2i-1}Y^{2j-1})$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{2j-1}\binom{2j-1}{k}r^{k}r^{\prime
2j-k-1}E\big{(}X^{2i+k-1}Z^{2j-k-1}\big{)}$ (65) $\displaystyle=$
$\displaystyle\sum_{l=1}^{j}\binom{2j-1}{2l-1}r^{2l-1}r^{\prime
2(j-l)}E\big{(}X^{2(i+l-1)}\big{)}\cdot E\big{(}Z^{2(j-l)}\big{)}$
$\displaystyle=$ $\displaystyle
r\sum_{l=1}^{j}\binom{2j-1}{2l-1}(r^{2})^{l-1}(1-r^{2})^{j-l}\cdot(2(i+l)-3)!!\cdot(2(j-l)-1)!!$
$\displaystyle:=$ $\displaystyle rh(r^{2}),$

where we set $k=2l-1$ in the second identity and the fact
$E[N(0,1)^{n}]=(n-1)!!$ for even number $n\geq 1$ is used in the last display.
Obviously the last term in (65) is obviously a polynomial. Then

$\displaystyle h(0)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{j}\binom{2j-1}{2l-1}(r^{2})^{l-1}(1-r^{2})^{j-l}\cdot(2(i+l)-3)!!\cdot(2(j-l)-1)!!\Big{|}_{r=0}$
$\displaystyle=$ $\displaystyle(2j-1)\cdot(2(i+1)-3)!!\cdot(2(j-1)-1)!!$
$\displaystyle=$ $\displaystyle(2i-1)!!\cdot(2j-1)!!.$

Set $f(x)=g(x)[(2i-1)!!\cdot(2j-1)!!]^{-1}$. The desired conclusion follows.

(ii) Recall $E(X^{2i})=(2i-1)!!$ and $E(Y^{2j})=(2j-1)!!$. The conclusion is
obviously true if $(i,j)=(0,0)$, $(i,j)=(0,1)$ and $(i,j)=(1,0)$, in which
cases $g(x)=0$ for each $x\in\mathbb{R}$. We next assume $i\geq 1$ and $j\geq
1$. By the same argument as in (i),

$\displaystyle
Y^{2j}=(rX+r^{\prime}Z)^{2j}=\sum_{k=0}^{2j}\binom{2j}{k}r^{k}r^{\prime
2j-k}X^{k}Z^{2j-k}.$

Then $E(X^{2i}Y^{2j})$ is equal to

$\displaystyle\sum_{k=0}^{2j}\binom{2j}{k}r^{k}r^{\prime
2j-k}E\big{(}X^{2i+k}Z^{2j-k}\big{)}$ $\displaystyle=$
$\displaystyle\sum_{l=0}^{j}\binom{2j}{2l}r^{2l}r^{\prime
2(j-l)}E\big{(}X^{2(i+l)}\big{)}\cdot E\big{(}Z^{2(j-l)}\big{)}$
$\displaystyle=$ $\displaystyle(1-r^{2})^{j}E\big{(}X^{2i}\big{)}\cdot
E\big{(}Z^{2j}\big{)}+r^{2}\sum_{l=1}^{j}\binom{2j}{2l}(r^{2})^{l-1}(1-r^{2})^{j-l}E\big{(}X^{2(i+l)}\big{)}\cdot
E\big{(}Z^{2(j-l)}\big{)},$

where we set $k=2l$ in the first identity and single out the term with $l=0$
in the last step. Obviously the last sum is a polynomial of $r^{2}$, say,
$h_{1}(r^{2}).$ Write $(1-r^{2})^{j}=1+r^{2}h_{2}(r^{2})$ with the function
$h_{2}(x)$ being a polynomial. Consequently,

$\displaystyle E(X^{2i}Y^{2j})$ $\displaystyle=$
$\displaystyle[1+r^{2}h_{2}(r^{2})]E(X^{2i})\cdot E(Z^{2j})+r^{2}h_{1}(r^{2})$
$\displaystyle=$ $\displaystyle(2i-1)!!\cdot(2j-1)!!+r^{2}g(r^{2}),$

where $g(x):=(2i-1)!!\cdot(2j-1)!!h_{2}(x)+h_{1}(x)$ is a polynomial. The
proof is completed. $\square$

With the help of Lemma 6.13, we will obtain Lemmas 6.14-6.16 in the following.
They estimate the size of
$|\mbox{Cov}(X_{1}^{d_{1}}X_{2}^{d_{2}},X_{3}^{d_{3}}X_{4}^{d_{4}})|$ for non-
negative integers $d_{1},d_{2},d_{3},d_{4}$ with special requirements.

######  LEMMA 6.14

Let $(X_{1},X_{2},X_{3},X_{4})^{T}$ be a $4$-dimensional random vector with
distribution $N_{4}(\mathbb{0},\mathbb{R})$ where
$\mathbb{R}=(r_{ij})_{4\times 4}$ and $r_{ii}=1$ for each $i$. Then, there
exists a constant $C>0$ depending on $i,j,k,l$ but not on $\mathbb{R}$ such
that $|\mbox{Cov}(X_{1}^{2i}X_{2}^{2j},X_{3}^{2k}X_{4}^{2l})|\leq C\sum_{1\leq
i<j\leq 4}r_{ij}^{2}$ for all non-negative integers $i,j,k$ and $l$.

Proof of Lemma 6.14. If $i+j=0$, then $i=j=0$, and hence
$X_{1}^{2i}X_{2}^{2j}=1$. The conclusion obviously holds. The same is true if
$k+l=0$. So we assume $i+j\geq 1$ and $k+l\geq 1$ next. Write

$\displaystyle\mbox{Cov}(X_{1}^{2i}X_{2}^{2j},X_{3}^{2k}X_{4}^{2l})=E\big{(}X_{1}^{2i}X_{2}^{2j}X_{3}^{2k}X_{4}^{2l}\big{)}-E\big{(}X_{1}^{2i}X_{2}^{2j}\big{)}\cdot
E\big{(}X_{3}^{2k}X_{4}^{2l}\big{)}.$ (66)

By Lemma 6.13(ii), there exists polynomials $g_{1}(x)$ and $g_{2}(x)$ for
$x\in\mathbb{R}$ such that

$\displaystyle
E\big{(}X_{1}^{2i}X_{2}^{2j}\big{)}=(2i-1)!!\cdot(2j-1)!!+r_{12}^{2}\cdot
g_{1}(r_{12}^{2});$ (67) $\displaystyle
E\big{(}X_{3}^{2k}X_{4}^{2l}\big{)}=(2k-1)!!\cdot(2l-1)!!+r_{34}^{2}\cdot
g_{2}(r_{34}^{2}).$ (68)

Set $r=(r_{12},r_{13},r_{14},r_{23},r_{24},r_{34})^{T}\in\mathbb{R}^{6}$.
Then, by Lemma 6.4, $E(X_{1}^{2i}X_{2}^{2j}X_{3}^{2k}X_{4}^{2l})$ is a
multivariate polynomial of $r_{ij},1\leq i<j\leq 4$. So we are able to write

$\displaystyle
E\big{(}X_{1}^{2i}X_{2}^{2j}X_{3}^{2k}X_{4}^{2l}\big{)}=C_{0}+\sum_{1\leq
i<j\leq 4}C_{ij}r_{ij}+c(r),$ (69)

where $C_{0}$ and $C_{ij}$ are constants and $c(r)$ is a linear combination of
$\prod_{1\leq i<j\leq 4}r_{ij}^{\alpha_{ij}}$ with $2\leq\sum_{1\leq i<j\leq
4}\alpha_{ij}\leq i+j+k+l.$ Next we will use the fact that (69) holds for
every $r_{ij}\in[-1,1]$ and $1\leq i<j\leq 4$ to identify the values of
$C_{0}$ and every $C_{ij}$.

Set $r_{ij}=0$ for all $i\neq j$. Then $X_{1},X_{2},X_{3},X_{4}$ are i.i.d.
$N(0,1)$. It follows from (69) that

$\displaystyle C_{0}$ $\displaystyle=$ $\displaystyle E(X_{1}^{2i})\cdot
E(X_{2}^{2j})\cdot E(X_{3}^{2k})\cdot E(X_{4}^{2l})$ $\displaystyle=$
$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!.$

Review (69). We claim that $C_{ij}=0$ for all $1\leq i<j\leq 4$. Take
$r_{ij}=0$ for all $1\leq i<j\leq 4$ except $(i,j)=(1,2)$. Then the three
random variables $(X_{1},X_{2})^{T}\in\mathbb{R}^{2},\,X_{3}$ and $X_{4}$ are
independent. By (69),

$\displaystyle E\big{(}X_{1}^{2i}X_{2}^{2j}\big{)}\cdot E(X_{3}^{2k})\cdot
E(X_{4}^{2l})$ (70) $\displaystyle=$
$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!+C_{12}r_{12}+c_{1}(r_{12}),$

where $c_{1}(x)=c_{2}x^{2}+c_{3}x^{3}+\cdots+c_{p}x^{p}$ with $2\leq p\leq
i+j+k+l$. Use the fact $E(X_{3}^{2k})=(2k-1)!!$, $E(X_{4}^{2l})=(2l-1)!!$ and
(67) to see

$\displaystyle E\big{(}X_{1}^{2i}X_{2}^{2j}\big{)}\cdot E(X_{3}^{2k})\cdot
E(X_{4}^{2l})$ $\displaystyle=$
$\displaystyle\big{[}(2i-1)!!\cdot(2j-1)!!+r_{12}^{2}\cdot
g_{1}(r_{12}^{2})\big{]}\cdot(2k-1)!!\cdot(2l-1)!!.$

Observe there are no linear terms of $r_{12}$ on the right hand side. Compare
this with (70), we see $C_{12}=0$. By symmetry, $C_{ij}=0$ for all $1\leq
i<j\leq 4$. Combining this, (66)-(69), we get
$\mbox{Cov}(X_{1}^{2i}X_{2}^{2j},X_{3}^{2k}X_{4}^{2l})$ is identical to

$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!+c(r)$ (71)
$\displaystyle-$ $\displaystyle\big{[}(2i-1)!!\cdot(2j-1)!!+r_{12}^{2}\cdot
g_{1}(r_{12}^{2})\big{]}\cdot\big{[}(2k-1)!!\cdot(2l-1)!!+r_{34}^{2}\cdot
g_{2}(r_{34}^{2})\big{]}$ $\displaystyle:=$ $\displaystyle c_{2}(r).$

Observe that $\mbox{Cov}(X_{1}^{2i}X_{2}^{2j},X_{3}^{2k}X_{4}^{2l})=c_{2}(r)$
is a linear combination of $\prod_{1\leq i<j\leq 4}r_{ij}^{\alpha_{ij}}$ with
$2\leq\sum_{1\leq i<j\leq 4}\alpha_{ij}\leq i+j+k+l.$ Notice $|r_{ij}|\leq 1$,
and hence $|r_{ij}^{\alpha_{ij}}|\leq r_{ij}^{2}$ for all $\alpha_{ij}\geq 2$
and $\prod_{1\leq i<j\leq
4}r_{ij}^{\alpha_{ij}}\leq|r_{kl}|^{\alpha_{kl}}|r_{uv}|^{\alpha_{uv}}$ for
any $1\leq k<l\leq 4$, $1\leq u<v\leq 4$, $(k,l)\neq(u,v)$ with
$\alpha_{kl}+\alpha_{uv}\geq 2$. Use the formula
$|xy|\leq\frac{1}{2}(x^{2}+y^{2})$ to get

$\displaystyle|r_{kl}|^{\alpha_{kl}}|r_{uv}|^{\alpha_{uv}}\leq
r_{kl}^{2}+r_{uv}^{2}\leq\sum_{1\leq i<j\leq 4}r_{ij}^{2}.$

This says that there exists a constant $C>0$ depending on $i,j,k,l$ but not on
$\mathbb{R}$ such that

$\displaystyle|\mbox{Cov}(X_{1}^{2i}X_{2}^{2j},X_{3}^{2k}X_{4}^{2l})|=|c_{2}(r)|\leq
C\sum_{1\leq i<j\leq 4}r_{ij}^{2}.$ (72)

The proof is completed. $\square$

######  LEMMA 6.15

Let $(X_{1},X_{2},X_{3},X_{4})^{T}$ be the same as in Lemma 6.14. Then, there
exists a constant $C>0$ depending on $i,j,k,l$ but not on $\mathbb{R}$ such
that $|\mbox{Cov}(X_{1}^{2i-1}X_{2}^{2j-1},X_{3}^{2k}X_{4}^{2l})|\leq
C\sum_{1\leq i<j\leq 4}r_{ij}^{2}$ for all non-negative integers $i,j,k,l$
with $i\geq 1$ and $j\geq 1$.

Proof of Lemma 6.15. If $k+l=0$, then $k=l=0$, and hence
$X_{3}^{2k}X_{4}^{2l}=1$. The conclusion trivially holds. So we assume and
$k+l\geq 1$ next.

By Lemma 6.13(i), there exists a polynomial $f(x)$ for $x\in\mathbb{R}$ with
$f(0)=1$ and

$\displaystyle E(X_{1}^{2i-1}X_{2}^{2j-1})=(2i-1)!!\cdot(2j-1)!!\cdot
r_{12}f(r_{12}^{2})$ (73)

for all integers $i\geq 1$ and $j\geq 1$. Write

$\displaystyle\mbox{Cov}(X_{1}^{2i-1}X_{2}^{2j-1},X_{3}^{2k}X_{4}^{2l})$ (74)
$\displaystyle=$ $\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k}X_{4}^{2l}\big{)}-E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}\big{)}\cdot
E\big{(}X_{3}^{2k}X_{4}^{2l}\big{)}.$

Then, by Lemma 6.4,

$\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k}X_{4}^{2l}\big{)}=D_{0}+\sum_{1\leq
i<j\leq 4}D_{ij}r_{ij}+d(r),$ (75)

where $D_{0}$ and $D_{ij}$’s are constants and $d(r)$ is a linear combination
of $\prod_{1\leq i<j\leq 4}r_{ij}^{\alpha_{ij}}$ with $2\leq\sum_{1\leq
i<j\leq 4}\alpha_{ij}\leq i+j+k+l-1.$ We now evaluate the values of $D_{0}$
and $D_{ij}$’s.

Take $r_{ij}=0$ for all $i\neq j$. Then $X_{1},X_{2},X_{3},X_{4}$ are i.i.d.
with distribution $N(0,1)$. From (75), we see

$\displaystyle D_{0}=E(X_{1}^{2i-1})\cdot E(X_{2}^{2j-1})\cdot
E(X_{3}^{2k})\cdot E(X_{4}^{2l})=0$

since $E(X_{1}^{2i-1})=0.$ Then (75) becomes

$\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k}X_{4}^{2l}\big{)}=\sum_{1\leq
i<j\leq 4}D_{ij}r_{ij}+d(r).$ (76)

Take $r_{ij}=0$ for all $1\leq i<j\leq 4$ except $(i,j)=(1,2)$. Then the three
random variables $(X_{1},X_{2})^{T}\in\mathbb{R}^{2},\,X_{3}$ and $X_{4}$ are
independent. It follows from (73) and (76) that

$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!\cdot
f(r_{12}^{2})r_{12}=D_{12}r_{12}+d_{1}(r_{12})$

where $d_{1}(x)=c_{2}x^{2}+c_{3}x^{3}+\cdots+c_{p}x^{p}$ with $2\leq p\leq
i+j+k+l-1$. Set $d_{2}(x)=c_{2}x+c_{3}x^{2}+\cdots+c_{p}x^{p-1}$. The above
implies

$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!\cdot
f(r_{12}^{2})=D_{12}+d_{2}(r_{12}).$

Take $r_{12}=0$. Recall $f(0)=1$. Then

$\displaystyle D_{12}=(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!.$

Now we claim $D_{ij}=0$ for all $(i,j)\neq(1,2).$ In fact, set $r_{ij}=0$ for
all $1\leq i<j\leq 4$ except $(i,j)=(1,3)$. Then $X_{2}$ is independent of
$(X_{1},X_{3},X_{4})^{T}$. Use the fact $E(X_{2}^{2j-1})=0$ and independence
to see from (76) that

$\displaystyle 0=D_{13}r_{13}+d_{3}(r_{13}),$ (77)

where $d_{3}(x)=a_{2}x^{2}+\cdots+a_{p}x^{p}$ with $2\leq p\leq i+j+k+l-1$.
Divide both sides of (77) by $r_{13}$ and then set $r_{13}=0$, we get
$D_{13}=0$. Similarly, $D_{ij}=0$ for all $(i,j)\neq(1,2).$ Therefore, by
(76),

$\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k}X_{4}^{2l}\big{)}=(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!\cdot
r_{12}+d(r).$ (78)

Since $f(x)$ in (73) is a polynomial with $f(0)=1$, we are able to write
$f(x)=1+xf_{1}(x)$ for each $x\in\mathbb{R}$, where $f_{1}(x)$ ia a
polynomial. Hence,

$\displaystyle E(X_{1}^{2i-1}X_{2}^{2j-1})$ $\displaystyle=$
$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot r_{12}+(2i-1)!!\cdot(2j-1)!!\cdot
r_{12}^{3}f_{1}(r_{12}^{2}).$

Joining this by (68), (74) and (78), we have that
$\mbox{Cov}(X_{1}^{2i-1}X_{2}^{2j-1},X_{3}^{2k}X_{4}^{2l})$ is equal to

$\displaystyle(2i-1)!!\cdot(2j-1)!!\cdot(2k-1)!!\cdot(2l-1)!!\cdot
r_{12}+d(r)$ $\displaystyle-$ $\displaystyle\big{[}(2i-1)!!\cdot(2j-1)!!\cdot
r_{12}+(2i-1)!!\cdot(2j-1)!!\cdot r_{12}^{3}f_{1}(r_{12}^{2})\big{]}$
$\displaystyle\cdot\big{[}(2k-1)!!\cdot(2l-1)!!+r_{34}^{2}\cdot
g_{2}(r_{34}^{2})\big{]}$ $\displaystyle:=$ $\displaystyle d_{4}(r),$

where $g_{2}(x)$ is a polynomial. Observe there are no linear terms of
$r_{ij},\,1\leq i<j\leq 4$ in the expression of $d_{4}(r)$. Consequently,
$d_{4}(r)$ is a linear combination of $\prod_{1\leq i<j\leq
4}r_{ij}^{\alpha_{ij}}$ with $2\leq\sum_{1\leq i<j\leq 4}\alpha_{ij}\leq
i+j+k+l-1.$ The conclusion then follows from the same argument between (71)
and (72). $\square$

######  LEMMA 6.16

Let $(X_{1},X_{2},X_{3},X_{4})^{T}$ be the same as in Lemma 6.14. Then, there
exists a constant $C>0$ depending on $i,j,k,l$ but not on $\mathbb{R}$ such
that $|\mbox{Cov}(X_{1}^{2i-1}X_{2}^{2j-1},X_{3}^{2k-1}X_{4}^{2l-1})|\leq
C\sum_{1\leq i<j\leq 4}r_{ij}^{2}$ for all positive integers $i,j,k,l$.

Proof of Lemma 6.16. Write

$\displaystyle\mbox{Cov}(X_{1}^{2i-1}X_{2}^{2j-1},X_{3}^{2k-1}X_{4}^{2l-1})$
(79) $\displaystyle=$ $\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}-E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}\big{)}\cdot
E\big{(}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}.$

By Lemma 6.13(i),

$\displaystyle E(X_{3}^{2k-1}X_{4}^{2l-1})=(2k-1)!!\cdot(2l-1)!!\cdot
r_{34}g(r_{34}^{2})$

for some polynomial $g(x)$. This together with (73) implies that

$\displaystyle\big{|}E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}\big{)}\cdot
E\big{(}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}\big{|}\leq C_{1}|r_{12}r_{34}|\leq
C\sum_{1\leq i<j\leq 4}r_{ij}^{2},$ (80)

where $C_{1}$ is a constant depending on $i,j,k,l$ but not on $\mathbb{R}$. By
Lemma 6.4,

$\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}=B_{0}+\sum_{1\leq
i<j\leq 4}B_{ij}r_{ij}+b(r),$ (81)

where $B_{0}$ and $B_{ij}$’s are constants and $b(r)$ is a linear combination
of $\prod_{1\leq i<j\leq 4}r_{ij}^{\alpha_{ij}}$ with $2\leq\sum_{1\leq
i<j\leq 4}\alpha_{ij}\leq i+j+k+l-2.$ We claim that

$\displaystyle B_{0}=B_{ij}=0$ (82)

for all $1\leq i<j\leq 4$. In fact, take $r_{ij}=0$ for all $1\leq i<j\leq 4$,
then $X_{1},X_{2},X_{3},X_{4}$ are i.i.d. $N(0,1)$-distributed random
variables. The left hand side of (81) is zero because of independence. The
right hand side of (81) is equal to $B_{0}$. This concludes $B_{0}=0$. Then
(81) becomes

$\displaystyle
E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}=\sum_{1\leq
i<j\leq 4}B_{ij}r_{ij}+b(r).$ (83)

For this identity we take $r_{ij}=0$ for all $1\leq i<j\leq 4$ except
$(i,j)=(1,2)$. Then the three random variables
$(X_{1},X_{2})^{T}\in\mathbb{R}^{2},\,X_{3}$ and $X_{4}$ are independent. By
independence and the fact $E(X_{3}^{2k-1})=0$, the left hand side of (83) is
zero. Hence, (83) is reduced to

$\displaystyle 0=B_{12}r_{12}+b_{1}(r_{12}),$

where $b_{1}(x)=c_{2}x^{2}+c_{3}x^{3}+\cdots+c_{p}x^{p}$ with $2\leq p\leq
i+j+k+l-2$. Set $b_{2}(x)=c_{2}x+c_{3}x^{2}+\cdots+c_{p}x^{p-1}$. Then
$0=B_{12}+b_{2}(r_{12}).$ Take $r_{12}=0$ to get $B_{12}=0$. By symmetry, we
know $B_{ij}=0$ for all $1\leq i<j\leq 4$. So (82) has been verified. It
follows that (81) is reduced to
$E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}=b(r)$. By
employing the same argument between (71) and (72), we get

$\displaystyle\big{|}E\big{(}X_{1}^{2i-1}X_{2}^{2j-1}X_{3}^{2k-1}X_{4}^{2l-1}\big{)}\big{|}\leq
C_{2}\sum_{1\leq i<j\leq 4}r_{ij}^{2},$

where $C_{2}$ is a constant depending on $i,j,k,l$ but not on $\mathbb{R}$.
This joined (79) and (80) yields the desired inequality. $\square$

The following fact supplies a convenient tool to handle the covariance between
products of independent random variables. Recall
$\|\xi\|_{q}=(E|\xi|^{q})^{1/q}$ for any random variable $\xi$ and $q\geq 1$.

######  LEMMA 6.17

For $k\geq 2$, let $\\{(U_{i},V_{i})\in\mathbb{R}^{2};\,1\leq i\leq k\\}$ be
independent random vectors. Assume
$C:=3\prod_{i=1}^{k}(1+\|U_{i}\|_{k})(1+\|V_{i}\|_{k})<\infty$ and $K=C^{k}$.
Then

$\displaystyle\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{k}U_{i},\prod_{i=1}^{k}V_{i}\Big{)}\Big{|}\leq
K\cdot\sum_{i=1}^{k}|\mbox{Cov}(U_{i},V_{i})|.$

Proof of Lemma 6.17. The inequality holds obviously for $k=1$. Assume now
$k\geq 2$. By definition,

$\displaystyle\mbox{Cov}(U_{1}U_{2},V_{1}V_{2})$ $\displaystyle=$
$\displaystyle E(U_{1}U_{2}V_{1}V_{2})-E(U_{1}U_{2})\cdot E(V_{1}V_{2})$ (84)
$\displaystyle=$ $\displaystyle E(U_{1}V_{1})\cdot E(U_{2}V_{2})-EU_{1}\cdot
EU_{2}\cdot EV_{1}\cdot EV_{2}.$

Write $E(U_{1}V_{1})=\mbox{Cov}(U_{1},V_{1})+EU_{1}\cdot EV_{1}$ and
$E(U_{2}V_{2})=\mbox{Cov}(U_{2},V_{2})+EU_{2}\cdot EV_{2}$. Plug the two
identities into (84) to see

$\displaystyle\mbox{Cov}(U_{1}U_{2},V_{1}V_{2})$ $\displaystyle=$
$\displaystyle\mbox{Cov}(U_{1},V_{1})\cdot\mbox{Cov}(U_{2},V_{2})+EU_{2}\cdot
EV_{2}\cdot\mbox{Cov}(U_{1},V_{1})$ (85)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+EU_{1}\cdot
EV_{1}\cdot\mbox{Cov}(U_{2},V_{2}).$

Since $k\geq 2$, we have from the Hölder inequality that
$E|\xi|\leq\|\xi\|_{F}\leq\|\xi\|_{k}$ for any random variable $\xi$. By the
Cauchy-Schwartz inequality,
$E|U_{1}V_{1}|\leq\|U_{1}\|_{2}\cdot\|V_{1}\|_{2}\leq\|U_{1}\|_{k}\cdot\|V_{1}\|_{k}$.
As a result, $|\mbox{Cov}(U_{1},V_{1})|$ is bounded by

$\displaystyle E|U_{1}V_{1}|+E|U_{1}|\cdot E|V_{1}|\leq
2\|U_{1}\|_{k}\|V_{1}\|_{k}.$ (86)

Set $\tau_{j}=3\prod_{i=1}^{j}(1+\|U_{i}\|_{k})(1+\|V_{i}\|_{k})$ for
$j=1,\cdots,k$. We claim that

$\displaystyle\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{k}U_{i},\prod_{i=1}^{k}V_{i}\Big{)}\Big{|}\leq(\tau_{1}\cdots\tau_{k})\cdot\sum_{i=1}^{k}|\mbox{Cov}(U_{i},V_{i})|.$
(87)

In fact, by applying (86) to (85), we see

$\displaystyle|\mbox{Cov}(U_{1}U_{2},V_{1}V_{2})|$ (88) $\displaystyle\leq$
$\displaystyle
2\|U_{1}\|_{k}\|V_{1}\|_{k}\cdot|\mbox{Cov}(U_{2},V_{2})|+\|U_{2}\|_{k}\|V_{2}\|_{k}\cdot|\mbox{Cov}(U_{1},V_{1})|+\|U_{1}\|_{k}\|V_{1}\|_{k}\cdot|\mbox{Cov}(U_{2},V_{2})|$
$\displaystyle\leq$
$\displaystyle\Big{[}3\prod_{i=1}^{2}\big{(}1+\|U_{i}\|_{k}\big{)}\big{(}1+\|V_{i}\|_{k}\big{)}\Big{]}\cdot\big{[}|\mbox{Cov}(U_{1},V_{1})|+|\mbox{Cov}(U_{2},V_{2})|\big{]}.$

So claim (87) holds for $k=2$ due to the fact $\tau_{j}\geq 1$ for each $j$.
Now we assume $k\geq 3$ and use induction to complete the proof. Assume

$\displaystyle\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{j}U_{i},\prod_{i=1}^{j}V_{i}\Big{)}\Big{|}\leq(\tau_{1}\cdots\tau_{j})\cdot\sum_{i=1}^{j}|\mbox{Cov}(U_{i},V_{i})|$
(89)

for some $2\leq j<k.$ By assumption, $(U_{j+1},V_{j+1})$ and
$(\prod_{i=1}^{j}U_{i},\prod_{i=1}^{j}V_{i})$ are independent. We obtain from
(88) that

$\displaystyle\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{j+1}U_{i},\prod_{i=1}^{j+1}V_{i}\Big{)}\Big{|}$
(90) $\displaystyle\leq$
$\displaystyle\Big{[}3\big{(}1+\|U_{j+1}\|_{k}\big{)}\big{(}1+\|V_{j+1}\|_{k}\big{)}\big{(}1+\|U_{1}\cdots
U_{j}\|_{k}\big{)}\big{(}1+\|V_{1}\cdots V_{j}\|_{k}\big{)}\Big{]}\cdot$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Big{[}|\mbox{Cov}(U_{j+1},V_{j+1})|+\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{j}U_{i},\prod_{i=1}^{j}V_{i}\Big{)}\Big{|}\Big{]}.$

Since $2\leq j<k$, then we have from the Hölder inequality that

$\displaystyle\|U_{1}\cdots U_{j}\|_{k}=\big{[}E\big{(}U_{1}\cdots
U_{j}\cdot\underbrace{1\cdots
1}_{k-j}\big{)}\big{]}^{1/k}\leq\prod_{i=1}^{j}\|U_{i}\|_{k}.$

The above also holds if the symbol “$U$” is replaced by “$V$”. Thus,

$\displaystyle 3\big{(}1+\|U_{1}\cdots
U_{j}\|_{k}\big{)}\big{(}1+\|V_{1}\cdots
V_{j}\|_{k}\big{)}\big{(}1+\|U_{j+1}\|_{k}\big{)}\big{(}1+\|V_{j+1}\|_{k}\big{)}\leq\tau_{j+1}.$

By (89) and (90),

$\displaystyle\Big{|}\mbox{Cov}\Big{(}\prod_{i=1}^{j+1}U_{i},\prod_{i=1}^{j+1}V_{i}\Big{)}\Big{|}$
$\displaystyle\leq$