arXiv:1801.06634v1 [] 20 Jan 2018

Joint CLT for eigenvalue statistics from
several dependent large dimensional sample
covariance matrices with application
Weiming Li1,* Zeng Li2,** and Jianfeng Yao3,†
1

School of Statistics and
Management
Shanghai University of Finance
and Economics
e-mail:

*

li.weiming@shufe.edu.cn
2

Department of Statistics
Pennsylvania State University
e-mail: ** zxl278@psu.edu

3

Department of Statistics and
Actuarial Science
The University of Hong Kong
e-mail: † jeffyao@hku.hk

Abstract: Let Xn = (xij ) be a k × n data matrix with complex-valued, independent
and standardized entries satisfying a Lindeberg-type moment condition. We consider
simultaneously R sample covariance matrices Bnr = n1 Qr Xn X∗n Q>
r , 1 ≤ r ≤ R,
where the Qr ’s are nonrandom real matrices with common dimensions p × k (k ≥
p). Assuming that both the dimension p and the sample size n grow to infinity,
the limiting distributions of the eigenvalues of the matrices {Bnr } are identified,
and as the main result of the paper, we establish a joint central limit theorem for
linear spectral statistics of the R matrices {Bnr }. Next, this new CLT is applied
to the problem of testing a high dimensional white noise in time series modelling.
In experiments the derived test has a controlled size and is significantly faster than
the classical permutation test, though it does have lower power. This application
highlights the necessity of such joint CLT in the presence of several dependent sample
covariance matrices. In contrast, all the existing works on CLT for linear spectral
statistics of large sample covariance matrices deal with a single sample covariance
matrix (R = 1).
Keywords and phrases: Large sample covariance matrices, central limit theorem,
linear spectral statistics, white noise test, high-dimensional times series.

MSC 2010 Mathematics Subject Classifications: Primary 62H10; secondary 60B12, 60B20.
1. Introduction
Modern information technology tremendously accelerates computing speed and greatly enlarges the amount of data storage, which enables us to collect, store and analyze data
1

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

2

of large dimensions. Classical limit theorems in multivariate analysis, which normally assume fixed dimensions, become no longer applicable for dealing with high dimensional
problems. Random matrix theory investigates the spectral properties of random matrices when their dimensions tend to infinity and hence provides a powerful framework for
solving high dimensional problems. This theory has made systematic corrections to many
classical multivariate statistical procedures in the past decades, see the monographs of Bai
and Silverstein (2010), Yao et al. (2015) and the review papers Johnstone (2007) and Paul
and Aue (2014). It has found diverse applications in various research areas, including signal processing, network security, image processing, statistical genetics and other financial
econometrics problems.
The sample covariance matrix is of central importance in multivariate analysis. Many
fundamental statistics in multivariate analysis can be written as functionals of eigenvalues
of a sample covariance matrix Sn such as linear spectral statistics (LSSs) of the form
f (λ1 ) + · · · + f (λp ) where the λj ’s are eigenvalues of Sn and f (·) is a smooth function.
The wide range of creditable applications in high dimensional statistics triggered an uptick
in the demand for CLTs of such LSSs. Actually one of the most widely used results in
this area is Bai and Silverstein (2004), which considers a sample covariance matrix of the
form Bn = n1 T1/2 Xn X∗n T1/2 , where Xn = (xij ) is p × n matrix consisting of i.i.d. complex
standardized entries and T is a p × p nonnegative Hermitian matrix. A CLT for LSSs
of Bn is established under the so-called Marčenko-Pastur regime, i.e. n, p → ∞, p/n →
c > 0. Further refinement and extensions can be found in Zheng et al. (2015), Chen and
Pan (2015), Zheng et al. (2017a), and Zheng et al. (2017b). Among them, Zheng et al.
(2015) studied the unbiased sample covariance matrix when the population mean vector
is unknown. Chen and Pan (2015) looked into the ultra-high dimensional case when the
dimension p is much larger than the sample size n, that is p/n → ∞ as n → ∞. Zheng et
al. (2017a) derived the CLT for LSSs of large dimensional general Fisher matrices. Zheng
et al. (2017b) attempted to relax the fourth order moment condition in Bai and Silverstein
(2004) and incorporated it into the limiting parameters.
However, this rich literature all deals with a single sample covariance matrix Bn . This
paper, on the contrary, aims at the joint limiting behaviour of functionals of several groups
of eigenvalues coming from different yet closely related sample covariance matrices. Specifically, we consider data samples {yjr }1≤j≤n,1≤r≤R of the form yjr = Qr xj where
(M1) {xj , 1 ≤ j ≤ n} is a sequence of k-dimensional independent and complex-valued
random vectors with independent standardized components (xij ), i.e. Exij = 0 and
E|xij |2 = 1, and the dimension k ≥ p;

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

3

(M2) {Qr , 1 ≤ r ≤ R} are R nonrandom real matrices with common dimensions p × k.
The R population covariance matrices {Tnr = Qr Q>
r , r = 1, . . . R} are assumed
product-commutative.
We thus consider R sample covariance matrices given by
n

Bnr

1
1X
∗
yjr yjr
= Qr Xn X∗n Q>
=
r , 1 ≤ r ≤ R,
n j=1
n

(1.1)

where Xn = (x1 , . . . , xn ) is of size k × n, ∗ denotes the conjugate transpose of matrices or
vectors, and > stands for the transpose of real ones. Let (λjr )1≤j≤p be the eigenvalues of
Bnr (1 ≤ r ≤ R), and consider L × R real-valued functions (flr )1≤l≤L, 1≤r≤R . This leads to
the family of L × R LSSs
ϕlr = flr (λ1r ) + · · · + flr (λpr ), 1 ≤ l ≤ L, 1 ≤ r ≤ R .
This paper establishes a joint CLT for these L × R statistics {ϕlr } under appropriate
conditions.
The importance of such joint CLT for LSSs is best explained and illustrated with the
following problem of testing a high dimensional white noise. Indeed, our motivation for the
joint CLT originates from this application to high-dimensional time series analysis. Testing
for white noise is a classical yet important problem in statistics, especially for diagnostic
checks in time series modelling. For high dimensional time series, current literatures focus
on estimation and dimension-reduction aspects of the modelling, including high dimensional VAR models and various factor models. Yet model diagnostics have largely been
untouched. Classical omnibus tests such as the multivariate Hosking and Li-McLeod tests
are no longer suitable for handling high dimensional time series. They become extremely
conservative, losing size and power dramatically. In a very recent work, Li et al. (2016)
looked into this high dimensional portmanteau test problem and proposed several new test
statistics based on single-lagged and multi-lagged sample auto-covariance matrices. More
precisely, let’s consider a p-dimensional time series modelled as a linear process
xt =

X

Al zt−l ,

(1.2)

l≥0

where {zt } is a sequence of independent p-dimensional random vectors with independent
components zt = (zit ) satisfying Ezit = 0, E|zit |2 = 1, E|zit |4 < ∞. Hence {xt } has
Ext = 0, and its lag-τ auto-covariance matrix Στ = Cov(xt+τ , xt ) depends on τ only. In
particular, Σ0 = Var(xt ) denotes the population covariance matrix of the series. The goal

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

4

is to test whether xt is a white noise, i.e.
H0 : Cov(xt+τ , xt ) = 0, τ = 1, · · · , q,

(1.3)

where q ≥ 1 is a prescribed constant integer. Let x1 , . . . , xn be a sample generated from
the model (1.2). The lag-τ sample auto-covariance matrix is
n

X
bτ = 1
Σ
xt x∗t−τ ,
n t=1

(1.4)

b τ . For any
where xt = xn+t when t ≤ 0. Li et al. (2016) proposed a test statistic based on Σ
eτ was designed to test the specific
given constant integer 1 ≤ τ ≤ q, their test statistic L
lag-τ autocorrelation of the sequence, i.e.
eτ =
L

p
X

f∗ M
f τ ),
λ2j,τ = Tr(M
τ

j=1

where {λj,τ , j = 1, · · · , p} are the eigenvalues of
n


1 X
1 b
∗
b
f
Στ + Στ =
xt x∗t−τ + xt−τ x∗t ,
Mτ =
2
2n t=1

which is the symmetrized lag-τ sample auto-covariance matrix.
f τ = 1 Xn (Dτ + D> )X∗ , where
Notice that in matrix form M
τ
n
2n
Dτ =

0 In−τ
Iτ
0

!

where Im denotes the mth order unit matrix. They have proved that, under the null
hypothesis, in the simplest setting when xt = zt , the limiting distribution of the test
eτ is
statistic L


ne
p d
1
3c(ν4 − 1)
Lτ − →
− N
, 1+
.
p
2
2
2
Here, p, n → ∞ and p/n → c > 0 and ν4 = E|zit |4 . The null hypothesis should be rejected
eτ . Simulation results also show that this test statistic is consistently
for large values of L
more powerful than the Hosking and Li-McLeod tests even when the latter two have been
size adjusted.
f τ , which can be studied with the
eτ is an LSS of M
It can be seen that the test statistic L
CLT in Bai and Silverstein (2004). Indeed, the non-null eigenvalues of the sample covariance
1/2
1/2
matrix Sn = n1 Tp Xn X∗n Tp considered there are the same as the matrix Sn = n1 Xn Tp X∗n

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

5

f τ . However, the test statistic L
eτ can only detect serial
which resembles to the matrix M
dependence in a single lag each time. In order to capture a multi-lag dependence structure,
a naturally more effective way would be accumulating the lags and consider the statistic
Lq =

q
X
τ =1

eτ =
L

q
X

fτ M
f ∗ ).
Tr(M
τ

(1.5)

τ =1

Note that the CLT in Bai and Silverstein (2004) (or in its recent extensions) can only
f τ , while to derive
be used to study the correlations between different LSSs of a given M
the null distribution of Lq , we need to study the correlations between LSSs of several
f τ , 1 ≤ τ ≤ q. Consequently, we need to resort to the joint CLT
covariance matrices M
eτ , 1 ≤ τ ≤ q}. It is worth
studied in this paper to characterize the correlations among {L
noticing that Li et al. (2016) proposed another multi-lagged test statistic Uq by stacking
a number of consecutive observation vectors. It will be shown in this paper that this test
statistic Uq is essentially much less powerful than Lq considered here due to the data loss
caused by observation stacking. This superiority of Lq over Uq demonstrates the necessity
and significance of studying a joint CLT for LSSs of several dependent sample covariance
matrices as proposed in this paper.
The rest of the paper is organized as follows. The main results of the joint CLT of LSSs
of different sample covariance matrices are presented in Section 2. The application on high
dimensional white noise test is provided in Section 3 to demonstrate the utility of this joint
CLT. Numerical studies have also lent full support to the theoretical derivations. Technical
lemmas and proofs are left to Section 4. Finally, Matlab codes for reproducing simulations
in the paper are available at: http://web.hku.hk/~jeffyao/papersInfo.html.
2. Joint CLT for linear spectral statistics of {Bnr }1≤r≤R
2.1. Preliminary knowledge on LSDs of {Bnr }1≤r≤R
Recall that the empirical spectral distribution (ESD) of a p × p square matrix A is the
P
probability measure F A = p−1 pi=1 δλi , where the λi ’s are eigenvalues of A and δa denotes
the Dirac mass at point a. For any probability measure F on the real line, its Stieltjes
transform is defined by
Z
1
m(z) =
dF (x), z ∈ C+ ,
x−z
+
where C denotes the upper complex plane.
The assumptions needed for the existence of limiting spectral distributions (LSDs) of
{Bnr }1≤r≤R are as follows. The setup as well as the following Lemma 2.1 are established
in Zheng et al. (2017b).

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

6

Assumption (a) Both dimensions p and n tend to infinity such that cn = p/n → c > 0
as n → ∞.
Assumption (b) Samples are {yjr = Qr xj , j = 1, . . . , n, r = 1, . . . , R}, where Qr is p×k,
xj = (x1j , . . . , xkj )> is k × 1, and the dimension k (k ≥ p) is arbitrary. Moreover,
{xij , i = 1, . . . , k, j = 1, . . . , n} is a k × n array of independent random variables, not
necessarily identically distributed, with common moments
Exij = 0, E|x2ij | = 1,
and satisfying the following Lindeberg-type condition: for each η > 0,
n
R
k


√
1 XXX
2
2
kqir k E|xij |I |xij | > η n/kqir k → 0,
pnη 2 i=1 j=1 r=1

where kqir k is the Euclidean norm of the i-th column vector qir of Qr .
Assumption (c) The ESD Hnr of the population covariance matrix Tnr = Qr Q>
r converges weakly to a probability distribution Hr , r = 1, . . . , R. Also the sequence of
the spectral norm of (Tnr ) is bounded in n and r.
Lemma 2.1. [Theorem 2.1 of Zheng et al. (2017b)] Under Assumptions (a)-(c), almost
surely, the ESD Fnr of Bnr weakly converges to a nonrandom LSD F c,Hr . Moreover, the
Stieltjes transform mr (z) of F c,Hr is the unique solution to the following Marčenko-Pastur
equation
Z
1
mr (z) =
dHr (t) ,
(2.1)
t[1 − c − czmr (z)] − z
on the set {mr (z) ∈ C : −(1 − c)/z + cmr (z) ∈ C+ }.
Define the companion LSD of Bnr as
F c,Hr = (1 − c)δ0 + cF c,Hr .
It is readily checked that F c,Hr is the LSD of the companion sample covariance matrix Bnr =
n−1 X∗n Q>
r Qr Xn (which is n × n), and its Stieltjes transform mr (z) = −(1 − c)/z + cmr (z)
satisfies the so-called Silverstein equation
Z
t
1
+c
dHr (t).
(2.2)
z=−
mr (z)
1 + tmr (z)
2.2. Main Results
Let A and B be two real symmetric p × p matrices satisfying AB = BA. The two matrices
can then be diagonalized simultaneously. We define the joint spectral distribution of (A, B)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

7

as the two-dimensional spectral distribution of the complex matrix A + iB, i.e.,
1
G(x, y) = #{i ≤ p, <(si ) ≤ x, =(si ) ≤ y},
p
where (si ) are the p eigenvalues of A + iB and #E denotes the cardinality of a set E.
Recall the random vector of L × R LSSs of Bnr ’s
Z

f`r (x)dFnr (x)
,
(2.3)
1≤`≤L,1≤r≤R

where (Fnr ) are the corresponding empirical spectral distributions of (Bnr ) and (f`r ) are
L × R measurable functions on the real line. Our aim in this section is to establish the joint
distribution of (2.3) under suitable conditions. The main results are presented as follows.
Assumption (d) The variables {xij , i = 1, . . . , k, j = 1, . . . , n} are independent, with
common moments
Exij = 0, E|x2ij | = 1, βx = E|x4ij | − |Ex2ij |2 − 2, and αx = |Ex2ij |2 ,
and satisfying the following Lindeberg-type condition: for each η > 0
k
n
R


p
1 XXX
2
4
n/kq
k
→ 0.
kq
k
E|x
|I
|x
|
>
η
ir
ir
ij
ij
pnη 6 i=1 j=1 r=1

(2.4)

Assumption (e) Either βx = 0, or the mixing matrices {Qr } are such that the matrices
{Q>
r Qr } are diagonal (therefore with arbitrary βx ).
Assumption (f ) The joint spectral distribution Hnrs of Tnr and Tns converges weakly
to a probability distribution Hrs , 1 ≤ r, s ≤ R.
The framework with Assumptions (d)-(e)-(f) is inspired by the one advocated in Zheng
et al. (2017b). However, an extension is necessary here since we are dealing with several
random matrices simultaneously while only one matrix is considered in the reference.
Theorem 2.1. Under Assumptions (a)-(f ), let f11 , . . . , fLR be L × R functions analytic
on a complex domain containing
√
√ 2
[I(0<c<1) (1 − c)2 lim inf λT
c) lim sup λT
(2.5)
min , (1 +
max ]
n

n

T
with T = {Tnr }, and λT
min and λmax denoting the smallest and the largest eigenvalue of all
the matrices in T, respectively. Then, the random vector
Z

Z
cn ,Hnr
p
f`r (x)dFnr (x) − f`r (x)dF
(x)
.
(2.6)
1≤`≤L,1≤r≤R

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

8

converges to an (L × R)-dimensional Gaussian random vector (Xf11 , . . . , XfLR ). The mean
function is


I
1
βx
αx
+
dz,
EXf`r = −
f`r (z)g1 (z)
2πi C1
(1 − g2 (z))(1 − αx g2 (z)) 1 − g2 (z)
where
Z
g1 (z) =

cm3r (z)t2
dHr (t) and
(1 + tmr (z))3

cm2r (z)t2
dHr (t).
(1 + tmr (z))2

Z
g2 (z) =

The covariance function is
1
Cov(Xf`0 r , Xf`s ) =
4π 2

I I
f`0 (z1 )f` (z2 )

∂ 2 g(z1 , z2 )
dz1 dz2 ,
∂z1 ∂z2

(2.7)

C1 C2

where g(z1 , z2 ) = log(1 − a(z1 , z2 )) + log(1 − αx a(z1 , z2 )) − βx a(z1 , z2 ) with
ZZ
cmr (z1 )ms (z2 )t1 t2
dHrs (t1 , t2 ).
a(z1 , z2 ) =
(1 + t1 mr (z1 ))(1 + t2 ms (z2 ))
The contours C1 and C2 are non-overlapping, closed, positively orientated in the complex
plane, and enclosing both the supports of F c,Hr and of F c,Hs .
Remark 1. The centralization term in (2.6) is the expectation of f with respect to the
distribution F cn ,Hnr . This distribution is a finite dimensional version of the LSD F c,Hr ,
which is defined by (2.1) with the parameters (c, Hr ) replaced with (cn , Hnr ). The use of
F cn ,Hnr instead of F c,Hr aims to eliminate the effect of the convergence rate of (cn , Hnr ) to
(c, Hr ).
As an illustrative example of Theorem 2.1, we consider a simplified case where only
two sample covariance matrices are involved, i.e. Xn X0n /n and QXn X0n Q> /n, where Xn
is a p × n matrix of i.i.d. real standard Gaussian variables. The corresponding population
covariance matrices are Ip and Tn := QQ> , respectively. It’s clear that the ESD and its
limit of the identity matrix Ip are both the Dirac measure δ1 . Those of Tn are general
and denoted by Hn and H, respectively. Moreover, the joint spectral distribution function
Hn12 (t1 , t2 ) of Ip and Tn is equal to Hn (t2 ) for t1 = 1 and zero otherwise. Denote the ESDs
of the two sample covariance matrices by Fn1 and Fn2 , respectively, and let
Gn1 (x) = Fn1 (x) − F cn ,δ1 (x) and Gn2 (x) = Fn2 (x) − F cn ,Hn2 (x).
Then for any analytic function f , we have
Z
>
Z
d
p
f (x)dGn1 (x), f (x)dGn2 (x)
→
− N

!
!!
v1
ψ11 ψ12
,
.
v2
ψ12 ψ22

(2.8)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

9

The parameters (v1 , v2 , ψ11 , ψ22 ) of the marginal distributions in (2.8) have been derived by
many authors, see Bai and Silverstein (2004) and Zheng et al. (2017b) for example. While
the covariance parameter ψ12 is new and, from Theorem 2.1, it can be formulated as
I I
c
f (z1 )f (z2 )(m1 (z1 ) + z1 m01 (z1 ))(m2 (z2 ) + z2 m02 (z2 ))
dz1 dz2 ,
v12 = − 2
2π C1 C2
[c − (1 + z1 m1 (z1 )(1 + z2 m2 (z2 )]2
where m1 (z) and m2 (z) are the companion Stieltjes transforms of the LSDs F c,δ1 (x) and
F c,H (x), respectively, and m0 (z) denotes the derivative of m(z) with respect to z. For the
R
simplest function f (z) = z, one may figure out v12 = 2c tdH(t) by the residual theorem.
3. Application to high dimensional white noise test
As discussed in the introduction, a notable application of the joint CLT presented in this
paper is to the high dimensional white noise test. In particular, it is expected that testing
power could be gained by accumulating information across different lags, that is, by using
P
fτ M
f ∗ ) defined in (1.5).
the test statistic Lq = qτ =1 Tr(M
τ
Define the scaled statistic
n
qp
φq = Lq − .
(3.1)
p
2
The null hypothesis will be rejected for large values of φq . We consider high-dimensional
situations where the dimension p is large compared to the sample size n. By applying the
CLT in Theorem 2.1, the asymptotic null distribution of φq is derived as follows.
Theorem 3.1. Let q ≥ 1 be a fixed integer, and assume that
1. {zit , i = 1, · · · , p, t = 1, · · · , n} is a set of i.i.d. real-valued variables satisfying
Ezit = 0, Ezit2 = 1, Ezit4 = ν4 < ∞;
2. Relaxed Marčenko-Pastur regime: both the sample size n and the dimension p grow
to infinity such that
p
p
0 < lim inf ≤ lim sup < ∞.
n→∞ n
n→∞ n
Then in the simplest setting where xt = zt , we have
q d
s(cn )−1/2 {φq − } →
− N (0, 1),
(3.2)
2
where s(u) = q + u(ν4 − 1)(q 2 + q/2).
The proof of this theorem is given in Section 4.
Let Zα be the upper-α quantile of the standard normal distribution at level α. Based on
Theorem 3.1, we obtain a procedure for testing the null hypothesis in (1.3) as follows.
p
q
Multi-Lag-q test: Reject H0 if φq − > Zα s(cn ).
(3.3)
2

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

10

3.1. Simulation Experiments
Most of the experiments of this section are designed to compare our test procedure in (3.3)
and the procedure based on the test statistic Uq from Li et al. (2016) using Simes’ method
(Simes, 1986). In Li et al. (2016), several testing procedures are discussed and the test Uq
performs quite satisfactorily in terms of both size and power across different scenarios.
More precisely, let q ≥ 1 be a fixed integer, define p(q + 1)-dimensional vectors yj =


xj(q+1)−q
h i


..
n
. Since Ext = 0 and Στ = Cov(xt+τ , xt ), we have
,
j
=
1,
.
.
.
,
N
,
N
=


.
q+1

xj(q+1)

Σ0 Σ1 · · · Σq
. 
.
Σ1 Σ0 . . .. 

.
.. . .

..
. Σ1 
.
.
Σq · · · Σ1 Σ0 (q+1)p×(q+1)p




Cov(yj ) = 



The null hypothesis H0 : Cov(xt+τ , xt ) = 0, τ = 1, · · · , q becomes H0 : Σ1 = · · · = Σq =
0, a test for a block diagonal covariance structure of the stacked sequence {yj }.
When Σ0 = σ 2 Ip , the white noise test of {xt } reduces to a sphericity test of {yj }.
The well known John’s test statistic Uq can be adopted for this purpose. In our case, the
corresponding John’s test statistic Uq is defined as
Uq =

1
p(q+1)

Pp(q+1)
i=1

lq

li,q − lq

2

2
,

where {li,q , i = 1, . . ., p(q + 1)} are the eigenvalues of the sample covariance matrix Sq =
PN
1
∗
j=1 yj yj , and l q is their average.
N
Notice however
h ithat the use of blocks above reduces the sample size n to the number of
n
blocks N = q+1
. This may result in a certain loss of power for the test. To limit such
loss of power, we adopt Simes’ method for multiple hypothesis testing in Simes (1986). To
implement Simes’ method, we denote
(0)

yj = yj (x1 , . . . , xn )
as the previously defined stacked sample. Then we rotate the sample (xj ) and define a
(k)
series of new stacked samples yj for k = 1, . . . , q, that is,
(k)

yj = yj (xk+1 , . . . , xn , x1 , . . . , xk )

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

11

Then John’s test statistic Uq can be calculated based on the q + 1 samples, which results
(k)
in q + 1 different statistics {Uq }. Moreover, let Pk , 0 ≤ k ≤ q, denote the (asymptotic)
P-value for the John’s test with the k-th set of yj ’s, i.e.

Pk = 1 − Φ (N Uq(k) − p(q + 1) − ν4 + 2)/2 ,
where Φ(·) is the cumulative distribution function of the standard normal distribution. Let
P(1) ≤ · · · ≤ P(q+1) be a permutation of P0 , . . . , Pq . Then by the Simes method, we reject
k
α at least for one 1 ≤ k ≤ q + 1 for the nominal level α.
H0 if P(k) ≤ q+1
To compare our test statistic φq with multi-lag-q John’s test statistic Uq , we set the
significance level α = 5% and the critical regions of the two tests are
p
(1) Our test φq : {φq > 2q + Z0.95 s(cn ) };
(2) Multi-lag-q John’s test Uq (using Simes’ method): {at least for one 1 ≤ k ≤ q +
k
1, P(k) ≤ q+1
0.05}.
Data are generated following four different scenarios for comparison:
i.i.d.

(I) Test size under Gaussian white noise: xt = zt , (zt ) ∼ Np (0, Ip );
i.i.d.
(II) Test size under Non-Gaussian white noise: xt = zt − 2, (zit ) ∼ Gamma(4, 0.5),
E(zit ) = 2, Var(zit ) = 1, ν4 (zit ) = 4.5;
(III) Test power under a Gaussian spherical AR(1) process: xt = Axt−1 + zt , A = aIp ,
i.i.d.
a = 0.1, (zt ) ∼ Np (0, Ip );
(IV) Test power under a Non-Gaussian spherical AR(1) process: xt = Axt−1 + (zt − 2),
i.i.d.
A = aIp , a = 0.1, (zit ) ∼ Gamma(4, 0.5), E(zit ) = 2, Var(zit ) = 1, ν4 (zit ) = 4.5.
Various (p, n)-combinations are tested to show the suitability of our test statistic for both
low and high dimensional settings. Empirical statistics are obtained using 2000 independent
replications. Table 1 compares the empirical sizes of the two tests φq and Uq . It can be
seen that both of them have reasonable sizes compared to the 5% nominal level across
all the tested (p, n)-combinations. Still, the two tests become slightly conservative under
Non-Gaussian distributions in Scenario (II) compared to the Gaussian case in Scenario (I).
A sample display of these sizes is given in Figure 1 (left panel).
In Table 2, we compare the power of the two tests. Our test φq displays a generally
much higher power than the multi-lag-q John’s test Uq , especially when the dimensions
(p, n) become larger. On the other hand, both tests have slightly lower power under the
Non-Gaussian distribution than under the Gaussian distribution, which is consistent with
the previous observation that the two tests become more conservative with Non-Gaussian
populations. A sample display of these powers is given in Figure 1 (right panel).

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

0.08

12

1
φq(q=1)

φq(q=1)

φq(q=3)

φq(q=3)

Uq(q=1)

0.07

0.8

Uq(q=3)

Uq(q=1)
Uq(q=3)

Size

Power

0.06
0.5

0.05

0.04

0.03
50

0.2

200

300

500

(p/n=0.5) dimension p

0
50

200

300

500

(p/n=0.5) dimension p

Fig 1. Sample plots of empirical sizes (left panel) from Table 1 (Scenario II with cn = 0.5) and

empirical powers (right panel) from Table 2 (Scenario III with cn = 0.5).

To further explore the powers of the two tests, we varied the AR coefficient a in Scenario (III) and (IV) from -0.1 to 0.1 (a = 0 corresponds to testing size). Smaller values
of the AR(1) coefficient a are used here leading to a more difficult testing problem and a
generally decreased power for both tests. Three dimensional settings are considered with
p/n ∈ {0.1, 0.5, 1.5} while the sample size is fixed as n = 600. The number of independent
replications is still 2000 in each case. Results for Scenario (III) and (IV) are plotted in
Figure 1. This Figure further consolidates that our test φq dominates Uq under all tested
scenarios. A nonnegligible increase in the testing power of both test statistics as the dimension p becomes larger sheds more light on the blessings of high dimensionality. Still both
tests are more conservative with Non-Gaussian population distribution than with Gaussian
distribution.
3.2. Comparison to a permutation test
As many complex analytic tools are employed to derive the asymptotic null distributions of
the test statistic φq , it is natural to wonder about the performance of a “simple-minded” test
procedure, namely the permutation test. Under the null hypothesis of white noise, since the

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

1

1
φq(q=1)

φq(q=1)
0.8

φq(q=3)

Power Values

Power Values

0.8

Uq(q=1)
0.5

Uq(q=3)

0.2
0.05
0
−0.1

Uq(q=1)
0.5

−0.05
0
0.05
(p/n=0.1) AR coefficient a

0.05
0
−0.1

0.1

Power Values

Power Values

Uq(q=3)

0.2

−0.05
0
0.05
(p/n=0.5) AR coefficient a

0.5

0.05
0
−0.1

0.1

−0.05
0
0.05
(p/n=0.5) AR coefficient a

0.1

φq(q=1)
0.8

φq(q=3)

Power Values

Power Values

Uq(q=3)

1

Uq(q=1)
Uq(q=3)

0.2
0.05
0
−0.1

φq(q=3)
Uq(q=1)

φq(q=1)

0.5

0.1

0.2

1
0.8

−0.05
0
0.05
(p/n=0.1) AR coefficient a

φq(q=1)
0.8

φq(q=3)
Uq(q=1)

0.05
0
−0.1

Uq(q=3)

1
φq(q=1)

0.5

φq(q=3)

0.2

1
0.8

13

φq(q=3)
Uq(q=1)

0.5

Uq(q=3)

0.2

−0.05
0
0.05
(p/n=1.5) AR coefficient a

0.1

0.05
0
−0.1

−0.05
0
0.05
(p/n=1.5) AR coefficient a

0.1

Fig 2. Empirical Powers for the two tests with varying AR coefficient a from -0.1 to 0.1. Left panel:
Scenario(III) for Gaussian Distribution. Right panel: Scenario(IV) for Gamma Distribution

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

14

sample vectors x1 , . . . , xn have an i.i.d. structure, one can permute these n sample vectors
say B times to obtain an empirical upper 5% quantiles of the test statistic φq . The null
hypothesis will be rejected if the observed statistic φq from the original (non permuted)
sample vectors x1 , . . . , xn is larger than this empirical quantile.
Data are generated following the spherical AR(1) process in Scenario (III) and (IV) to
compare this straightforward test with our test statistic φq . In order to compare the power
performance of two tests, the AR coefficient a takes different values, a = [0, 0.05, 0.09, 0.1]
(a = 0 corresponds to testing size). The sample size is fixed as n = 300 yet data dimension
p varies. As for the permutation test, the permutation times is set as B = 500. The nominal
level is α = 5%. Testing size and power of two tests are shown in Tables 3 and 4 based on
500 replicates for all (p, n) configurations.
It can be seen that the sizes of both tests are well controlled. As for their power, our test
offers an acceptable level while the permutation test consistently performs better in the
tested cases. However, the permutation test is extremely time consuming compared to our
test. For instance, to run one set of (p, n) = (150, 300) combination for 500 replicates, it
takes only 25 seconds with our test, while almost 3 hours for the permutation test with
permutation times B = 500. Particularly the computation time increases greatly when the
sample size n grows. Therefore, our test statistic φq provides a very competitive choice for
testing high dimensional white noise while the classical permutation test is simpler, more
powerful though much slower.
4. Proofs of the main theorems
4.1. Proof of Theorem 2.1
The general strategy for our main Theorem 2.1 follows the methods advocated in Bai and
Silverstein (2004), with its most recent update in Zheng et al. (2017b). However, as we
are dealing with several random matrices simultaneously, all the technical steps for the
implementation of this strategy have to be carefully rewritten. They are presented in this
section.
4.1.1. Sketch of the proof of Theorem 2.1
Let v0 > 0 be arbitrary, xr be any number greater than the right end point of interval
(2.5), and xl be any negative number if the left end point of (2.5) is zero, otherwise choose
√ 2
c) ). Define a contour C as
xl ∈ (0, lim inf p→∞ λT
min (1 −
C = {x + iv : x ∈ {xr , xl }, v ∈ [−v0 , v0 ]} ∪ Cu with Cu = {x ± iv0 : x ∈ [xl , xr ]},

(4.1)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

15

and let Cn = Cu ∪ {x ± iv : x ∈ {xl , xr }, v ∈ [n−1 εn , v0 ]} with εn ≥ n−α for some α ∈ (0, 1).
By definition, the contour C encloses a rectangular region in the complex plane, which
contains the union of the support sets of all the LSDs F c,Hr , 1 ≤ r ≤ R. As a regularized
version of C, Cn excludes a small segment near the real line.
Let mnr (z), mnr (z), m0nr (z), m0nr (z) be the Stieltjes transforms of Fnr , F nr , F cn ,Hnr , and
F cn ,Hnr , respectively, where Fnr is the ESD of Bnr , F cn ,Hnr is the LSD defined in Remark
1, F and F are linked by the equation F = (1 − cn )δ0 + cn F . A major task of proving
Theorem 2.1 is to study the convergence of the empirical process
Mnr (z) := p[mnr (z) − m0nr (z)] = n[mnr (z) − m0nr (z)].
To this end, we need to truncate Mnr (z) as



z ∈ Cn ,

Mnr (z)
cnr (z) = M (x + in−1 ε ) x ∈ {x , x } and v ∈ [0, n−1 ε ],
M
nr
n
l
r
n



−1
−1
Mnr (x − in εn ) x ∈ {xl , xr } and v ∈ [−n εn , 0],
which agrees to Mnr (z) on Cn . This truncation is essential when proving the tightness
cnr (z) on C. Write
M


c
c
c
Mn (z) = Mn1 (z), . . . , MnR (z) ,
we will establish its convergence as stated in the following lemma.
c n (·) converges weakly to a Gaussian process
Lemma 4.1. Under Assumptions (a)-(f ), M
M(·) = (M1 , . . . , MR )(·) on C. The mean function is
EMr (z) =

αx g1 (z)
βx g1 (z)
+
,
(1 − g2 (z))(1 − αx g2 (z)) 1 − g2 (z)

where
Z
g1 (z) =

cm3r (z)t2
dHr (t) and
(1 + tmr (z))3

Z
g2 (z) =

cm2r (z)t2
dHr (t).
(1 + tmr (z))2

The covariance function is
Cov(Mr (z1 ), Ms (z2 )) = −

∂2
[log(1 − a(z1 , z2 )) + log(1 − αx a(z1 , z2 )) − βx a(z1 , z2 )]
∂z1 ∂z2

where
ZZ
a(z1 , z2 ) = c

t1 t2 mr (z1 )ms (z2 )
dHrs (t1 , t2 ).
(1 + t1 mr (z1 ))(1 + t2 ms (z2 ))

From this lemma, Theorem 2.1 follows by similar arguments on Pages 562 and 563 in
Bai and Silverstein (2004).

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

16

4.1.2. Proof of Lemma 4.1
Following closely the steps of truncation, centralization and rescaling in Appendix B of
Zheng et al. (2017b), one may find that it is sufficient to prove this lemma under the
assumption that
√
ηn n
|xij | <
,
(4.2)
max1≤r≤R {kqir k}
where the constant ηn → 0 as n → ∞.
Write for r ∈ {1, . . . , R} and z ∈ Cn ,
Mnr (z) = p[mnr (z) − Emnr (z)] + p[Emnr (z) − m0nr (z)]
1
2
:= Mnr
(z) + Mnr
(z).

The Lemma can be proved by verifying three conditions (Bai and Silverstein, 2004):
Condition 1: Finite dimensional convergence of Mn1 (z) in distribution;
Condition 2: Tightness of Mn1 (z) on Cn ;
Condition 3: Convergence of Mn2 (z).
Since the second and third conditions can be obtained directly from Lemma 5.1 in Zheng
et al. (2017b), we only consider the first one by showing that, for any W × R complex
1
(zjr ))1≤j≤W,1≤r≤R converges to a Gaussian
numbers z11 , . . . , zW R , the random vector (Mnr
vector. Without loss of generality, we assume max{kQr k} ≤ 1. We will also denote by K
any constants appearing in inequalities and K may take on different values for different
expressions.
√
With the notation rjr = (1/ n)Qr xj , we define some quantities:
Dr (z) = Bnr − zI, Djr (z) = Dr (z) − rjr r∗jr , Dijr (z) = Dr (z) − rir r∗ir − rjr r∗jr ,
−1
−1
jr (z) = r∗jr D−1
jr (z)rjr − n trTnr Djr (z),
−2
−1
δjr (z) = r∗jr D−2
jr (z)rjr − n trTnr Djr (z),
1
1
1
, βijr (z) =
, β̄jr (z) =
,
βjr (z) =
−1
−1
∗
∗
−1
1 + rjr Djr (z)rjr
1 + rir Dijr (z)rir
1 + n trTnr D−1
jr (z)
1
1
bnr (z) =
,
b
(z)
=
,
12r
1 + n−1 EtrTnr D−1
1 + n−1 EtrTnr D−1
r (z)
12r (z)

which will be frequently used in the sequel. Note that quantities in the last two rows are
all bounded in absolute value by |z|/=(z).
1
By martingale difference decomposition, the process Mnr
(z) can be expressed as
−1
p[mnr (z) − Emnr (z)] = tr[D−1
r (z) − EDr (z)]

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

=

n
X

17

−1
−1
−1
trEj [D−1
r (z) − Djr (z)] − trEj−1 [Dr (z) − Djr (z)]

j=1
n
X
= −
(Ej − Ej−1 )βjr (z)r∗jr D−2
jr (z)rjr
j=1
n

dX
(Ej − Ej−1 ) log βjr (z)
= −
dz j=1
n

dX
=
(Ej − Ej−1 ) log 1 + jr (z)β̄jr (z) ,
dz j=1
−1
−1
−1
∗
where the third equality is from the identity D−1
r (z) = Djr (z) − Djr (z)rjr rjr Djr (z)βjr (z)
and the last one is obtained using the identity βjr (z) = β̄jr (z)[1 + β̄jr (z)jr (z)]−1 . We next
show that
n
n

dX
dX
(Ej − Ej−1 ) log 1 + jr (z)β̄jr (z) −
(Ej − Ej−1 )jr (z)β̄jr (z) = op (1).
dz j=1
dz j=1

Considering the second moment of the above difference, by the Cauchy integral formula,
one may get
2

n



dX
E
(Ej − Ej−1 ) log 1 + jr (z)β̄jr (z) − jr (z)β̄jr (z)
dz j=1



I
n
X
log 1 + jr (ζ)β̄jr (ζ) − jr (ζ)β̄jr (ζ)
1
=E
dζ
2
2πi
(z
−
ζ)
|ζ−z|=v/2
j=1
n I
K X
E|jr (ζ)β̄jr (ζ)|4 |dζ|.
≤ 2 4
π v j=1 |ζ−z|=v/2

2

(4.3)

From Lemma A.1 and the truncation in (4.2), we have


k

2 X
K
−1
−1
4
8
> −1
E|jr (ζ)| ≤ 4 E trTnr Djr (ζ)Tnr Djr (ζ̄) +
E|xij | E|qir Djr (ζ)qir |
n
i=1
4

≤ Kn−2 + Kηn4 n−1 ,
by which the right hand side of (4.3) tends to zero. Therefore, we need only to consider
the limiting distribution of
n
n
dX
dX
(Ej − Ej−1 )jr (z)β̄jr (z) =
Ej jr (z)β̄jr (z)
dz j=1
dz j=1

(4.4)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

18

in finite dimensional situations. To verify the Lyapunov condition, one can show that

n
2 
X
d
d
E (Ej − Ej−1 ) jr (z)β̄jr (z) I (Ej − Ej−1 ) jr (z)β̄jr (z) ≥ 
dz
dz
j=1
n

1X
d
≤ 2
E Ej jr (z)β̄jr (z)
 j=1
dz

4

→ 0,

where the convergence is again from Lemma A.1 and (4.2). Hence, from the martingale
1
(Billingsley, 1995, Theorem 35.12), the random vector (Mnr
(zjr )) will tend to a Gaussian
vector (Mr (zjr )) with covariance function
n
X



∂
∂
jr (z1 )β̄jr (z1 ) · Ej
js (z2 )β̄js (z2 ) . (4.5)
Cov(Mr (z1 ), Ms (z2 )) = lim
Ej−1 Ej
n→∞
∂z1
∂z2
j=1
We note that the referenced martingale CLT applies also to multidimensional martingale
by considering arbitrary linear combination of its components.
Using the same approach of Bai and Silverstein (2004) on Page 571, one may replace
β̄jr (z) by bnr (z). Then, by (1.15) of Bai and Silverstein (2004), we have
Γnrs (z1 , z2 ):=

n
X

bnr (z1 )bns (z2 )Ej−1 [Ej jr (z1 )Ej js (z2 )]

j=1

"
n
1 X
−1
> −1
= 2
bnr (z1 )bns (z2 ) trEj Q>
r Djr (z1 )Qr Ej Qs Djs (z2 )Qs
n j=1
−1
> −1
+ αx trEj Q>
r Djr (z1 )Qr Ej Qs Djs (z2 )Qs

+ βx

k
X

#
−1
−1
>
q>
ir Ej Djr (z1 )qir qis Ej Djs (z2 )qis

i=1

:= Γ1 + αx Γ2 + βx Γ3 ,

(4.6)

where αx = |Ex211 |2 and βx = E|x411 | − |Ex211 |2 − 2.
Now we derive the limit of the first term in (4.6). The means is to replace D−1
jr (z) (and
−1
similarly Djs (z)) by a proper nonrandom matrix. For this, we introduce such a one
Lr (z) = zI −

n−1
b12r (z)Tnr ,
n

whose inverse spectral norm is bounded, that is,
||Lr (z)||−1 ≤

|b−1
1 + p/(nv)
12r (z)|
≤
.
−1
v
=(zb12r (z))

(4.7)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

19

−1
−1
∗
We will show that the major part of D−1
jr (z) is just −Lr (z). From the identity rir Djr (z) =
βijr (z)r∗ir D−1
ijr (z), we get
−1
−1
−1
D−1
jr (z) + Lr (z) = Lr (z) (Djr (z) + Lr (z)) Djr (z)

X

= L−1
r (z)

i6=j

X

= L−1
r (z)

n−1
b12r (z)Tnr
rir r∗ir −
n
βijr (z)rir r∗ir D−1
ijr (z) −

i6=j

!
D−1
jr (z)

n−1
b12r (z)Tnr D−1
jr (z)
n

= b12r (z)R1r (z) + R2r (z) + R3r (z),

!

(4.8)

where
R1r (z) =

X

−1
∗
−1
L−1
r (z)(rir rir − n Tnr )Dijr (z),

i6=j

X
−1
∗
(βijr (z) − b12r (z))L−1
R2r (z) =
r (z)rir rir Dijr (z),
i6=j

R3r (z) = n−1 b12r (z)L−1
r (z)Tnr

X


−1
D−1
ijr (z) − Djr (z) .

i6=j
−1
From this decomposition, after substituting −L−1
r (z) for Djr (z) in the first term in (4.6),
there are three remaining quantities. Let’s check which one (or ones) of them can be
omitted. From Lemma A.3, (4.7), and (4.3) of Bai and Silverstein (1998), for any p × p
matrix M, we have
2
1/2

E|trR2r (z)M| ≤ nE

2

1/2

(|β12r (z) − b12r (z)| )E

≤ n1/2 K|||M|||

−1
r∗1r D−1
12r MLr (z)r1r

|z|2 (1 + p/(nv))
,
v5

(4.9)

where |||M||| denotes a nonrandom bound on the spectral norm of M. From Lemma A.2,
|trR3r (z)M| ≤ |||M|||

|z|(1 + p/(nv))
.
v3

(4.10)

Again from Lemma A.3 and (4.7), for nonrandom M,
−1
−1
−1
−1
2
E|trR1r (z)M| ≤ nE1/2 |r∗ir D−1
ijr (z)MLr (z)rir − n trTnr Dijr (z)MLr (z)|

≤ n1/2 K||M||

(1 + p/(nv))
.
v2

(4.11)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

20

Therefore, quantities containing R2r (z) and R3r (z) are both negligible. For the quantity
−1
−1
−1
∗
involving R1r (z), applying the identity D−1
js (z) = Dijs (z) − Dijs (z)rjs rjs Dijs (z)βijs (z), it
can be divided into three parts, that is,
> −1
trQ>
r Ej (R1r (z1 ))Qr Qs Djs (z2 )Qs = R11 (z1 , z2 ) + R12 (z1 , z2 ) + R13 (z1 , z2 ),

(4.12)

where
R11 (z1 , z2 ) = −

X

−1
> −1
∗
> −1
βijs (z2 )r∗ir Ej (D−1
ijr (z1 ))Qr Qs Dijs (z2 )ris ris Dijs (z2 )Qs Qr Lr (z1 )rir ,

i<j

R12 (z1 , z2 ) = −tr

X

−1
−1
−1
−1
>
>
L−1
r (z1 )n Tnr Ej (Dijr (z1 ))Qr Qs (Djs (z2 ) − Dijs (z2 ))Qs Qr ,

i<j

R13 (z1 , z2 ) = tr

X

−1
∗
−1
> −1
>
L−1
r (z1 )(rir rir − n Tnr )Ej (Dijr (z1 ))Qr Qs Dijs (z2 )Qs Qr .

i<j

From Lemma A.2 and (4.7) we get |R12 (z1 , z2 )| ≤ (1 + p/(nv0 ))/v03 , and similar to (4.9),
E|R13 (z1 , z2 )| ≤ n1/2 (1 + p/(nv0 ))/v03 . Thus these two parts are trivial. We then turn to
dealing with R11 (z1 , z2 ). Using Lemma A.1, Lemma A.3, and (4.3) of Bai and Silverstein
(1998) we get, for i < j,
−1
> −1
∗
> −1
E βijs (z2 )r∗ir Ej (D−1
ijr (z1 ))Qr Qs Dijs (z2 )ris ris Dijs (z2 )Qs Qr Lr (z1 )rir



−1
> −1
>
> −1
>
−b12s (z2 )n−2 tr Ej (D−1
ijr (z1 ))Qr Qs Dijs (z2 )Qs Qr tr Dijs (z2 )Qs Qr Lr (z1 )Qr Qs }
≤ Kn−1/2 .
So we may simplify R11 (z1 , z2 ) by replacing βijs (z2 ) with b12s (z2 ) and remove the random
parts of rir and ris . By Lemma A.2, we have


−1
> −1
>
> −1
>
tr Ej (D−1
ijr (z1 ))Qr Qs Dijs (z2 )Qs Qr tr Dijs (z2 )Qs Qr Lr (z1 )Qr Qs }


−1
> −1
>
> −1
>
≤ Kn.
−tr Ej (D−1
jr (z1 ))Qr Qs Djs (z2 )Qs Qr tr Djs (z2 )Qs Qr Lr (z1 )Qr Qs }
−1
−1
It implies that we may further replace D−1
ijr (z1 ) and Dijs (z2 ) in R11 (z1 , z2 ) with Djr (z1 )
and D−1
js (z2 ), respectively, which yields

E R11 (z1 , z2 )
+



j−1
−1
> −1
> −1
> −1
b12s (z2 )tr Q>
r Ej (Djr (z1 ))Qr Qs Djs (z2 )Qs tr Qs Djs (z2 )Qs Qr Lr (z1 )Qr
2
n

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

≤ Kn1/2 .

21

(4.13)

Integrating the results in (4.8)-(4.13), we obtain

−1
> −1
tr Q>
r Ej (Djr (z1 ))Qr Qs Djs (z2 )Qs



j−1
> −1
> −1
b12r (z1 )b12s (z2 )tr Qs Djs (z2 )Qs Qr Lr (z1 )Qr
× 1+
n2
> −1
−1
= −trQ>
r Lr (z1 )Qr Qs Djs (z2 )Qs + R14 (z1 , z2 ),

(4.14)

where E|R14 (z1 , z2 )| ≤ Kn1/2 . Furthermore, from this and (4.8)-(4.11), we may substitute
−1
for the second and third D−1
js (z2 ) in (4.14) with −Ls (z2 ) and then get

−1
> −1
tr Q>
E
(D
(z
))Q
Q
D
(z
)Q
j
1
r
2
s
r
s
jr
js



j−1
> −1
> −1
b12r (z1 )b12s (z2 )tr Qs Ls (z2 )Qs Qr Lr (z1 )Qr
× 1−
n2
−1
> −1
= trQ>
r Lr (z1 )Qr Qs Ls (z2 )Qs + R15 (z1 , z2 ),

(4.15)

where E|R15 (z1 , z2 )| ≤ Kn1/2 .
From Lemma A.2 and (4.3) of Bai and Silverstein (1998), we have
|b12r (z) − bnr (z)| ≤ Kn−1

and |bnr (z) − Eβ1r (z)| ≤ Kn−1/2 ,

respectively. By (2.2) of Silverstein (1995) and discussions in Section 5 of Bai and Silverstein
(1998), we have
Eβ1r (z) = −zEmnr (z) and |Emnr (z) − m0nr (z)| ≤ Kn−1 ,
respectively. Therefore, we get
|b12r (z) + zm0nr (z)| ≤ Kn−1/2 .
Combining this and (4.15), it follows that

j−1 0
tr
1−
mnr (z1 )m0ns (z2 )
n2


>
0
−1
>
0
−1
×tr Qs (I + mns (z2 )Tns ) Qs Qr (I + mnr (z1 )Tnr ) Qr


j−1
−1
>
> −1
= tr Qr Ej (Djr (z1 ))Qr Qs Djs (z2 )Qs 1 −
cn m0nr (z1 )m0ns (z2 )
n

Z Z
t1 t2 dHprs (t1 , t2 )
×
(1 + t1 m0nr (z1 ))(1 + t2 m0ns (z2 ))
−1
> −1
Q>
r Ej (Djr (z1 ))Qr Qs Djs (z2 )Qs



(4.16)

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

ncn
=
z1 z2

Z Z

t1 t2 dHprs (t1 , t2 )
+ R16 (z1 , z2 ),
(1 + t1 m0nr (z1 ))(1 + t2 m0ns (z2 ))

22

(4.17)

where E|R16 (z1 , z2 )| ≤ Kn1/2 . Similar to the arguments in Bai and Silverstein (2004) (page
577), using (4.17) and letting
Z Z
t1 t2 dHprs (t1 , t2 )
0
0
an (z1 , z2 ) = cn mnr (z1 )mns (z2 )
,
(1 + t1 m0nr (z1 ))(1 + t2 m0ns (z2 ))
we get
n

1X
an (z1 , z2 )
i.p.
Γ1 =
+ op (1) −→ − log(1 − a(z1 , z2 )),
n j=1 1 − ((j − 1)/n)an (z1 , z2 )

(4.18)

where
Z Z
a(z1 , z2 ) = cmr (z1 )ms (z2 )

t1 t2 dHrs (t1 , t2 )
.
(1 + tr mr (z1 ))(1 + ts ms (z2 ))

Similar to the derivation of Γ1 , one can easily show that
n

1X
αx an (z1 , z2 )
i.p.
αx Γ2 =
+ op (1) −→ − log(1 − αx a(z1 , z2 )).
n j=1 1 − αx ((j − 1)/n)an (z1 , z2 )

(4.19)

Considering the third term of (4.6), βx Γ3 with βx 6= 0, from Assumption (e), the matrix
> −1
Q>
r Qr is diagonal, so is Qr Lr Qr . Using (4.8)-(4.11) and (4.16), we have
n
k
1 XX > >
−1
> >
βx Γ3 = βx bnr (z1 )bns (z2 ) 2
e Q Ej (D−1
jr (z1 ))Qr ei ei Qs Ej (Djs (z2 ))Qs ei
n j=1 i=1 i r
n
k
1 X X > > −1
> −1
= βx bnr (z1 )bns (z2 ) 2
ei Qr Lr (z1 )Qr ei e>
i Qs Ls (z2 )Qs ei + op (1)
n j=1 i=1
n

1 X
−1
> −1
tr Q>
= βx bnr (z1 )bns (z2 ) 2
r Lr (z1 )Qr Qs Ls (z2 )Qs + op (1)
n j=1
i.p.

−→ βx a(z1 , z2 ).

(4.20)

Collecting the results in (4.5), (4.18)-(4.20), we finally get the covariance function in the
Lemma and the proof of this Lemma is completed.

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

23

4.2. Proof of Theorem 3.1
First we show that it is it is enough to establish the following claim: under the (classical)
Marčenko-Pastur regime, i.e., n → ∞, p = pn → ∞ such that p/n → c > 0, it holds that
φq,n −

q d
→
− N (0, s(c)) ,
2

(4.21)

where recall that s(u) = q + u(ν4 − 1)(q 2 + 2q ). Here we use φq,n for φq to signify the
dependence in n. So assume this claim is true. Under the relaxed Marčenko-Pastur regime,
the sequence {pn /n} is bounded below and above. For any subsequence (pnk , nk )k of (pn , n),
we can extract a further subsequence (pnk` , nk` )` such that the ratios pnk` /nk` converge to
α > 0 when ` → ∞. On this subsequence, by Claim (4.21),
φq,nk` −

q d
→
− N (0, s(α)) ,
2

` → ∞.

By continuity of the function u → s(u), we have

q d
−1/2
→
− N (0, 1) ,
s(pnk` /nk` )
φq,nk` −
2

` → ∞.

As this limit is independent of the subsequence and it holds for all such subsequences, the
same limit holds for the whole sequence, that is,

q d
→
− N (0, 1) , ` → ∞.
s(pn /n)−1/2 φq,n −
2
The required asymptotic normality is thus established.
The remaining of the section is devoted to a proof of Claim (4.21) assuming p/n → c > 0.
Define the banded Toeplitz matrix


0 · · · 21 · · · 0
 .. . .
.. 
1
 .
. 0
. 
2


 1

.
1
.
Cn,τ = 
. 0 2 
0
2


. 
 .. 1
..
. .. 
0
 . 2
0 · · · 12 · · · 0 n×n
and

n
 1
1 X
b
Nτ =
xt x∗t−τ + xt−τ x∗t = Xn Cn,τ X∗n .
2p t=1+τ
p

Define the associated Fourier series f (λ) of the banded Toeplitz matrix Cn,τ as
f (λ) = lim

n→∞

τ
X
k=−τ

tk eikλ =


1 iτ λ
e + e−iτ λ = cos(τ λ),
2

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

24

where tk is entry on the k-diagonal of Cn,τ .
According to the fundamental eigenvalue distribution theorem of Szegö for Toeplitz
forms, see Section 5.2 in Grenander and Szegö (1958) and Theorem 4.1 in Gray (2006), we
can infer that
Lemma 4.2. Suppose {lt , t = 1, · · · , n} are eigenvalues of Cn,τ with Fourier series f (λ),
then
(1) For any positive integer s,
n

1
1X s
lt =
lim
n→∞ n
2π
t=1

Z
0

2π

1
f (λ) dλ =
2π
s

Z

2π

(cos(τ λ))s dλ,

0

(2) For any continuous function on support of {lt , t = 1, · · · , n},
n

1
1X
F (lt ) =
lim
n→∞ n
2π
t=1

Z

2π

F (cos(τ λ)) dλ,
0

(3) Sequence {lt , t = 1, · · · , n} and


 
2πτ t
, t = 1, · · · , n
cos
n
are asymptotically equally distributed.
The limiting spectral distribution of Cn,τ is also derived in Lemma 3.1 of Bai and Wang
(2015), this useful lemma is stated as follows:
Lemma 4.3. As T → ∞, the ESD of Cn,τ tends to H, which is an Arcsine distribution
with density function
1
H 0 (t) = √
, t ∈ (−1, 1).
π 1 − t2
Recall for the permutation matrices D1 and Dτ = Dτ1 defined in Introduction, it holds
that
>
>
>
D1 D>
1 = D1 D1 = Dτ Dτ = Dτ Dτ = In .
Meanwhile, from the properties of Chebyshev polynomials, we can derive the following
lemma.

has eigenvalue
Lemma 4.4. (1) 12 D1 + D>
1
 


2πt
cos
, t = 1, · · · , n ,
n

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

(2)

1
2

25


Dτ + D>
has eigenvalue
τ
  


2πt
Tτ cos
, t = 1, · · · , n ,
n

where Tτ (·) stands for the Chebyshev polynomial of order τ .

shares the same asymptotic spectral distribution with Cn,τ as n → ∞.
(3) 12 Dτ + D>
τ
Since
n

X

eτ = 1
xt x∗t−τ + xt−τ x∗t
N
2p t=1
!
q
X

1
Dτ + D>
(x1 , · · · , xn )∗
(x1 , · · · , xn )
=
τ
2p
τ =1
!
q
X

1
= Xn
Dτ + D>
X∗n ,
τ
2p
τ =1
here for t ≤ 0, xt = xn+t , by Lemmas 4.2 and 4.4, it doesn’t take too much effort to see
e τ and N
b τ share the same limiting spectral distribution.
that N
e τ , by
Consider the Stieltjes transform mτ (z) of the limiting spectral distribution of N
implementing the Silverstein equation (2.2), we can infer that mτ (z) satisfies
Z
1
1
t
z=−
+
dHτ (t)
mτ (z) c
1 + tmτ (z)
!
1
1
1
−1 + − p
=
,
mτ (z)
c c 1 − m2τ (z)
where p/n → c > 0 as n → ∞, which coincides with the results in Bai and Wang (2015)
and Li et al. (2016).
Note that our test statistic
Lq =

q
X
τ =1



eτ =
L

q
X
τ =1

fτ M
f∗ ) =
Tr(M
τ

q
 p 2 X

n

e τN
e ∗ ),
Tr(N
τ

τ =1




fτ = 1 Σ
bτ + Σ
b ∗ = 1 Pn xt x∗ + xt−τ x∗ , thus the asymptotic properties
where M
τ
t−τ
t
t=1
2
2n
f
e τ since M
fτ = p N
e .
of Mτ can be inferred from those of N
n τ
eτ
Actually in Li et al. (2016), the asymptotic behavior of the single-lag-τ test statistic L
has already been thoroughly explored and characterized. Theorem 2.1 in Li et al. (2016) is
stated as follows:
Lemma 4.5. Let τ ≥ 1 be a fixed integer, and assume that

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

26

1. {zit , i = 1, · · · , p, t = 1, · · · , n} are all independently distributed satisfying Ezit =
0, Ezit2 = 1, Ezit4 = ν4 < ∞;
2. (Marčenko-Pastur regime). The dimension p and the sample size n grow to infinity
in a related way such that cn := p/n → c > 0.
eτ
Then in the simplest setting when xt = zt , the limiting distribution of the test statistic L
is


3(ν4 − 1)
p d
1
ne
, 1+
c .
Lτ − →
− N
p
2
2
2
Pq e
Now consider the multi-lag-q test statistic Lq =
τ =1 Lτ , combining with Lemma
4.5,
 all we need is the joint distribution of any two different single-lag test statistic, i.e.
er , L
es , 1 ≤ r 6= s ≤ q.
L
For a given integer q > 0, 1 ≤ r 6= s ≤ q, let fr (x) = fs (x) = x2 ,
 ∗
 ∗
1
1
Dr + D>
Ds + D>
r Xn Xn , Bns =
s Xn Xn ,
2p
2p


e r = 1 Xn Dr + D> X∗ , B ns = N
e s = 1 Xn Ds + D> X∗ ,
=N
r
n
s
n
2p
2p

Bnr =
Bnr

then both the LSDs of Bnr and Bns , i.e. F c,Hr , F c,Hs have Stieltjes transform mr (z) and
ms (z) satisfying the equation
!
1
1
1
.
z=
−1 + − p
m(z)
c c 1 − m2 (z)
Meanwhile,
e rN
e ∗) =
Tr(N
r

Z

e sN
e ∗) =
fr (x) dFnr (x), Tr(N
s

Z
fs (x) dFns (x).

where Fnr and Fns are the ESDs of Bnr and Bns . Thus, by directly implementing our joint
CLT for linear spectral statistics ofthe samplecovariance matrices,
i.e. Theorem 2.1, we

e1 , · · · , L
eq or any pair of L
er , L
es , 1 ≤ r 6= s ≤ q.
can derive the joint distribution of L
Precisely, the covariance function in Theorem 2.1 for the present case can be calculated to
be

3

 1+ 2 c(β2 x +2) , if r = s,
c
(4.22)
Cov (Xfr , Xfs ) =

 βx +2 ,
if
r
=
6
s.
c
The details of this lengthy derivation are postponed to Appendix B. Combining with
Lemma 4.5, for any given integer q ≤ 1, it can be inferred that, under the same assumptions

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

27



e1 , · · · , L
eq is
in Lemma 4.5, the joint limiting distribution of L




e1
c(ν4 − 1)
1 + 3(ν42−1) c · · ·
L
n .  p
d

1

..
..
...
− N  1q , 
 ,
 ..  − 1q →
.
.
p
2
2
eq
c(ν4 − 1) · · · 1 + 3(ν42−1) c
L


where 1q = (1, · · · , 1)> is a q dimensional vector with q ones.
Recall that
q
q
 p 2 X
X
e ∗ ),
e τN
e
Tr(N
Lq =
Lτ =
τ
n
τ =1
τ =1
then by the Delta method, we can derive the limiting distribution of our test statistic Lq ,
i.e.,
q
n
pq d
q 
Lq −
→
− N
, q + c(ν4 − 1)(q 2 + ) .
p
2
2
2
Claim (4.21) is thus established.
5. Discussions
In this paper we have introduced, for the first time in the literature on eigenvalues of
large sample covariance matrices, a joint central limit theorem involving several population covariance matrices. This theorem is believed to provide wide applications to current
problems in high-dimensional statistics, especially for testing on structures of population
covariances. As a show-case, we treated the problem of testing for a high-dimensional
white noise in time series modelling. The derived new test shows very promising performance compared to existing competitors For future study, it would be worth investigating
other significant applications of this CLT.
Acknowledgments. We thank a Referee for the suggestion of the comparison to a
permutation test made in Section 3.1. Weiming Li’s research is partially supported by
National Natural Science Foundation of China, No. 11401037, MOE Project of Humanities
and Social Sciences, No. 17YJC790057, and Program of IRTSHUFE. Jianfeng Yao thanks
support from HKSAR GRF Grant 17305814.
References
Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis. 3rd Edition.
John Wiley & Sons.

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

28

Bai, Z. D., (1999) Methodologies in spectral analysis of large dimensional random matrices,
a review. Statistica Sinica, 9, 611-677.
Bai, Z. D., Jiang, D. D., Yao, J. F. and Zheng, S. R. (2009). Corrections to LRT on large
dimensional covariance matrix by RMT. Annals of Statistics, 37, 3822-3840.
Bai, Z. D., Jiang, D. D., Yao, J. F., and Zheng, S. R. (2013). Testing linear hypotheses in
high-dimensional regressions. Statistics, 47, 1207-1223.
Bai, Z. D. and Silverstein, J. W. (1998). No eigenvalues outside the support of the limiting
spectral distribution of large-dimensional sample covariance matrices. Ann. Probab., 26,
316-345.
Bai, Z. D. and Silverstein, J. W. (2004). CLT for linear spectral statistics of large dimensional sample covariance matrices. Ann. Probab., 32, 553-605.
Bai, Z. D. and Silverstein, J. W. (2010). Spectral Analysis of Large Dimensional Random
Matrices. Science Press: Beijing.
Bai, Z. D. and Wang, C. (2015). A note on the limiting spectral distribution of a symmetrized auto-cross covariance matrix. Statistics & Probability Letters, 96, 333-340.
Bai, Z. D. and Zhou, W. (2008). Large sample covariance matrices without independence
structures in columns. Statistica Sinica, 28, 425-442.
Chen, B. B. and Pan, G. M. (2015) CLT for linear spectral statistics of normalized sample
covariance matrices with the dimension much larger than the sample size. Bernoulli, 21,
1089-1133.
Billingsley, P. (1995). Probability and Measure. John Wiley & Sons: New York
Burkholder, D. L. (1973). Distribution function inequalities for martingales. Ann. Probab.,
1, 19-42.
Gray, R. M. (2006). Toeplitz and Circulant Matrices: a Review. Now Publishers Inc.
Grenander, U. and Szegö, G.(1958). Toeplitz Forms and Their Applications. In: California
Monographs in Mathematical Sciences. University of California Press, Berkeley.
Johnstone, I. M. (2007). High dimensional statistical inference and random matrices. Int.
Cong. Mathematicians, Vol. I, 307-333. Zürich, Switzerland: European Mathematical
Society.
Li, Z., Yao, J., Lam, C., and Yao, Q. (2016). On testing a high-dimensional white noise.
Manuscript.
Marčenko, V. A. and Pastur, L. A. (1967). Distribution of eigenvalues for some sets of
random matrices. Math. USSR-Sb, 1, 457-483.
Pan, G. M. (2012). Comparison between two types of large sample covariance matrices.
Annales de l’Institut Henri Poincare-Probabiliteset Statistiques, 50, 655-677.
Paul, D. and Aue, A. (2014). Random matrix theory in statistics: A review, Journal of

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

29

Statistical Planning and Inference, 150, 1-29.
Silverstein, J. W. (1995). Strong convergence of the empirical distribution of eigenvalues
of large dimensional random matrices. J. Multivariate Anal., 5, 331-339.
Silverstein, J. W. and Bai, Z. D. (1995). On the empirical distribution of eigenvalues of a
class of large dimensional random matrices. J. Multivariate Anal., 54, 175-192.
Silverstein, J. W. and Choi, S. I. (1995). Analysis of the limiting spectral distribution of
large dimensional random matrices. J. Multivariate Anal., 54, 295-309.
Simes, R. J. (1986). An improved Bonferroni procedure for multiple tests of significance.
Biometrika, 73, 751-754.
Srivastava, M. S. (2005). Some tests concerning the covariance matrix in high dimensional
data. J. Japan Statist. Soc., 35, 251–272.
Yao, J. F., Bai, Z. D., and Zheng, S. R. (2015). Large Sample Covariance Matrices and
High-Dimensional Data Analysis, Cambridge University Press.
Zheng, S. R., Bai, Z. D., and Yao, J. F. (2015). Substitution principle for CLT of linear
spectral statistics of high-dimensional sample covariance matrices with applications to
hypothesis testing. Ann. Statist., 43, 546-591.
Zheng, S. R., Bai, Z. D., and Yao, J. F. (2017a). CLT for large dimensional general Fisher
matrices and its applications in high-dimensional data analysis. Bernoulli, 23, 1130-1178.
Zheng, S. R., Bai, Z. D., Yao, J. F., and Zhu, H. T. (2017b). CLT for linear spectral statistics
of large dimensional sample covariance matrices with dependent data. Manuscript.

Appendix A: Mathematical Tools
Lemma A.1. Let X = (X1 , . . . , Xn )> be a (complex) random vector with independent and
standardized entries having finite fourth moment and C = (cij ) be n × n (complex) matrix.
We have
!
n
X
E|X∗ CX − trC|4 ≤ K [tr(CC∗ )]2 +
E|Xii8 ||cii |4 .
i=1

The proof of the lemma follows easily by simple calculus and thus omitted.
Lemma A.2. (Lemma 2.6 of Silverstein and Bai (1995)). Let z ∈ C+ with v = = z, A
and B being n × n with B Hermitian, and r ∈ Cn . Then

kAk
r∗ (B − zI)−1 A(B − zI)−1 r
≤
.
tr (B − zI)−1 − (B + rr∗ − zI)−1 A =
∗
−1
1 + r (B − zI) r
v

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

30

Lemma A.3. [Formula 2.3 of Bai and Silverstein (2004)] For any nonrandom p × p
e ` , ` = 1, . . . q2 .
matrices Ck , k = 1, . . . , q1 and C
!
q2
q1
Y
Y
e ` r1r − n−1 trTnr C
e `)
E
(r∗ C
r∗ Ck r1r
1r

1r
`=1
q1
q2
Y
Y
−(1∧q2 ) (q2 −2)∨0
e ` ||,
Kn
δn
||Ck ||
||C
k=1
`=1
k=1

≤

q1 , q2 ≥ 0,

(A.1)

where K is a positive constant depending on q1 and q2 .
Appendix B: Derivation of the covariance (4.22)
Applying Theorem 2.1 to the functions fr (x) = fs (x) where 1 ≤ r 6= s ≤ q, the
corresponding covariance function is
I I
∂ 2 g(z1 , z2 )
1
dz1 dz2 ,
(B.1)
f
(z
)f
(z
)
Cov(Xfr , Xfs ) =
r 1 s 2
4π 2
∂z1 ∂z2
C1 C 2

where g(z1 , z2 ) = log(1 − a(z1 , z2 )) + log(1 − αx a(z2 , z2 )) − βx a(z1 , z2 ) with
ZZ
cmr (z1 )ms (z2 )t1 t2
a(z1 , z2 ) =
dHrs (t1 , t2 ).
(1 + t1 mr (z1 ))(1 + t2 ms (z2 ))
Mapping into our case, we have
1
p ↔ n, n ↔ p, c ↔ , αx = 1, βx = ν4 − 3, fr (x) = fs (x) = x2 ,
c
1
a(z1 , z2 ) =
c
1
=
c

ZZ

mr (z1 )ms (z2 )t1 t2
dHrs (t1 , t2 )
(1 + t1 mr (z1 ))(1 + t2 ms (z2 ))

1

m1 m2 Tr (t)Ts (t)

−1

(1 + Tr (t)m1 )(1 + Ts (t)m2 )

Z

dH(t),

where m1 , mr (z1 ), m2 , ms (z2 ), Tr (t), Ts (t) are Chebyshev polynomial of order r and
s respectively. Furthermore, we have
∂a (z1 , z2 )
1
=
∂z1
c

Z

∂a (z1 , z2 )
1
=
∂z2
c

Z

1

−1

1

−1

Ts (t) m2
Tr (t)
∂m1
·
dH(t),
2 ·
1 + Ts (t) m2 (1 + Tr (t) m1 )
∂z1
Tr (t) m1
Ts (t)
∂m2
·
dH(t),
2 ·
1 + Tr (t) m1 (1 + Ts (t) m2 )
∂z2

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

1
∂ 2 a (z1 , z2 )
=
∂z1 ∂z2
c

Z

1

Tr (t)
Ts (t)
∂m1 ∂m2
·
dH(t).
2 ·
2 ·
∂z1 ∂z2
(1 + Tr (t) m1 ) (1 + Ts (t) m2 )

−1

Since g(z1 , z2 ) = 2 log (1 − a(z1 , z2 )) − βx a(z1 , z2 ),
2

∂ a
∂a
· ∂a
∂ 2 log (1 − a(z1 , z2 ))
∂z1 ∂z2
∂z1 ∂z2
,
=−
−
∂z1 ∂z2
1 − a (1 − a)2

Thus
1
Cov (Xfr , Xfs ) = 2
4π

I I
C1

z12 z22

C2

∂ 2 g (z1 , z2 )
dz1 dz2
∂z1 ∂z2

I I
∂2a
∂a
· ∂a
1
1
2 2 ∂z1 ∂z2
=− 2
dz1 dz2 − 2
z1 z2
z12 z22 ∂z1 ∂z22 dz1 dz2
2π C1 C2
1−a
2π C1 C2
(1 − a)
I I
2
βx
∂ a
− 2
dz1 dz2 , M1 + M2 + M3 .
z12 z22
4π C1 C2
∂z1 ∂z2
I I

Consider M1 first, by Cauchy’s residue theorem, we have
I
1
z12
∂ 2a
·
dz1
2πi C1 1 − a ∂z1 ∂z2
I
Z
1
Ts (t)
∂m1 ∂m2
z12
1 1
Tr (t)
=
·
dz1
·
2 ·
2 dH(t) ·
2πi C1 1 − a c −1 (1 + Tr (t) m1 ) (1 + Ts (t) m2 )
∂z1 ∂z2

Z 
I
1
∂m2 1 1
Tr (t)Ts (t)
z12
·
dm1 dH(t)
=−
·
∂z2 c −1 (1 + Ts (t) m2 )2 2πi C1 (1 − a) (1 + Tr (t) m1 )2



2
1
1
Z
I
−1 + c − √

c 1−m21
1
∂m2 1 1 
 Tr (t)Ts (t)

=−
·
·
dm

1  dH(t)
∂z2 c −1  (1 + Ts (t) m2 )2 2πi C1 m21 (1 − a) (1 + Tr (t) m1 )2


∂m2 1
=−
·
∂z2 c




 T (t)T (t)

r
s

2

−1 (1 + Ts (t) m2 )


Z

1



2 (1)

1
1
 −1 + c − c√1−m2

1


·

 (1 − a) (1 + Tr (t) m1 )2 




 dH(t).


m1 =0

Note that
∂a
1
=
∂m1
c

Z

1

−1

Ts (t)m2
Tr (t)
·
dH(t),
1 + Ts (t) m2 (1 + Tr (t) m1 )2

31

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

32

then
2 (1)



1
1

 −1 + c − c√1−m2
1



2
 (1 − a) (1 + Tr (t) m1 ) 

1
= −2Tr (t) +
c

Z

1

−1

Ts (u)Tr (u)m2
dH(u),
1 + Ts (u) m2

m1 =0

therefore,
I I
∂ 2a
1
1
2 2
·
dz1 dz2
z z ·
M1 = − 2
2π C1 C2 1 2 1 − a ∂z1 ∂z2


I
Z 1
Z
2
Tr (t)Ts (t)
1 1 Ts (u)Tr (u)m2
2 1
=
z ·
−2Tr (t) +
dH(u) dH(t)dm2
2πi C2 2 c −1 (1 + Ts (t) m2 )2
c −1 1 + Ts (u) m2


2
I  −1 + 1c − √ 1 2
Z

c 1−m2
1 1 −2Tr2 (t)Ts (t)
2


·
dH(t)
=
 dm2

2
2
2πi C2 
m2
c −1 (1 + Ts (t) m2 )


1
+
c

1
=
c

Z

1

 2 I


−1  2πi C2

Z

1

−1


−1 + 1c − √ 1

2




Tr (t)Ts (t)  Z 1

1
Ts (u)Tr (u)

·
dH(u)
dm
2  dH(t)
2
c
1
+
T
(u)
m

m2 (1 + Ts (t) m2 )
s
−1
2
c

1−m22

 
2 (1)
1
1
2 −1 + c − √

 
c 1−m22


2
−2Tr (t)Ts (t) · 

2


(1 + Ts (t) m2 )

dH(t)
m =0

2

 Z 1
Z 1
2
1
Ts (u)Tr (u)dH(u) ·
Tr (t)Ts (t)dH(t)
+
c −1
c −1
Z 1
2
Z
8 1 2
2
2
=
T (t)Ts (t)dH(t) + 2
Ts (u)Tr (u)dH(u) .
c −1 r
c
−1

Similarly, for M2 , considering
I
1
1
∂a ∂a
z12 ·
·
dz1
2·
2πi C1
(1 − a) ∂z1 ∂z2

2
1
1
 Z 1

I
−1 + c − √
c 1−m21
1
1
Ts (t) m2
Tr (t)
∂m1
=
·
·
·
dH(t)
2πi C1
c −1 1 + Ts (t) m2 (1 + Tr (t) m1 )2 ∂z1
m21 (1 − a)2
 Z 1

1
Tr (t) m1
Ts (t)
∂m2
·
·
·
dH(t) dz1
c −1 1 + Tr (t) m1 (1 + Ts (t) m2 )2 ∂z2

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics


∂m2 1
=−
·
∂z2 c

1

Z

 1 I


−1  2πi C1

2


−1 + 1c − √ 1
c

33

1−m21
2

m1 (1 − a)

·

Tr (t)Ts (t) m2
(1 + Ts (t) m2 ) (1 + Tr (t) m1 )2

 Z 1


1
Tr (u) Ts (u)
·
dH(u) dm1 dH(t)
c −1 (1 + Tr (u) m1 ) (1 + Ts (u) m2 )2
Z
Z
1 1 Tr (u) Ts (u)
∂m2 1 1 Tr (t)Ts (t) m2
·
dH(t) ·
dH(u).
=−
∂z2 c −1 1 + Ts (t) m2
c −1 (1 + Ts (u) m2 )2
Then,
I I
1
1
∂a ∂a
M2 = − 2
z12 z22 ·
·
dz1 dz2
2·
2π C1 C2
(1 − a) ∂z1 ∂z2


I
Z
Z
1 1 Tr (u) Ts (u)
2
∂m2 1 1 Tr (t)Ts (t) m2
2
dH(t) ·
dH(u) dz2
=
z −
·
2πi C2 2
∂z2 c −1 1 + Ts (t) m2
c −1 (1 + Ts (u) m2 )2



2
1
1
 Z 1

Z
I Tr (t)Ts (t) −1 + c − √

c 1−m22
1 1
Tr (u) Ts (u)
1
 2

=
dH(u)
dm

2  dH(t)
c −1  2πi C2
m2 (1 + Ts (t) m2 )
c −1 (1 + Ts (u) m2 )2

2
= 2
c

Z

2

1

Tr (t)Ts (t) dH(t)

.

−1

As for M3 ,
I
1
∂ 2a
z12
dz1
2πi C1 ∂z1 ∂z2

2
1
1

 Z 1
I
−1 + c − √
c 1−m21
Tr (t)Ts (t)
∂m1 ∂m2
1
1
=
·
dH(t)
·
dz1
2
2
2πi C1
m21
c −1 (1 + Tr (t) m1 ) (1 + Ts (t) m2 )
∂z1 ∂z2



2
1
1
Z
−1 + c − √
 1 I

c 1−m21
Tr (t)Ts (t)
∂m2 1 1


=−
·
·
dm
1  dH(t)
2 
2
2
∂z2 c −1 (1 + Ts (t) m2 )  2πi C1 m1 (1 + Tr (t) m1 )


∂m2 1
= −
·
∂z2 c

Z

1

−1

Tr (t)Ts (t)
(1 + Ts (t) m2 )2

1
c

√1
c 1−m21

 −1 + −

·

(1 + Tr (t) m1 )2

2 (1)





dH(t)
m1 =0

=

∂m2 1
·
∂z2 c

Z

1

−1

2Tr2 (t)Ts (t)
(1 + Ts (t) m2 )2

dH(t),

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

thus
βx
M3 =
2πi

I

z22



C2

∂m2 1
·
∂z2 c


Z

1

 1 I


−1  2πi C2

βx
= −
c

1


2Tr2 (t)Ts (t)
2 dH(t) dz2
−1 (1 + Ts (t) m2 )


2
1
1
−1 + c − √

c 1−m22
2Tr2 (t)Ts (t)

·
dm
2  dH(t)
2
m22

(1 + Ts (t) m2 )
Z


βx
= −
c

Z

Z

1

1
c

√1
c 1−m22

2 (1)

 −1 + −



(1 + Ts (t) m2 )2

1

2Tr2 (t)Ts (t)

−1






dH(t)
m2 =0

=

4βx
c

Tr2 (t)Ts2 (t)dH(t).

−1

Combining the three terms gives
4(βx + 2)
Cov (Xfr , Xfs ) =
c

Z

1

Tr2 (t)Ts2 (t)dH(t)

−1

1
Note that for dH(t) = √
dt, we have
π 1 − t2
Z 1
Tr (t)Ts (t) dH(t) =
−1
Z 1
Tr2 (t)Ts2 (t) dH(t) =
−1

4
+ 2
c

Z

1

2
Ts (u)Tr (u)dH(u) .

−1

1
1{r=s} ,
2
3
1
1{r=s} + 1{r6=s} .
8
4

Therefore,
Cov (Xfr , Xfs ) =
The required formula is established.





1+ 23 c(βx +2)
,
c2

if r = s,




βx +2
,
c

if r 6= s.

34

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

Table 1
Test Size under Scenario (I) and (II)
p
5
10
25
40
10
20
50
80
50
100
250
400
10
40
60
100
50
200
300
500
90
360
540
900
150
600
900
1500
200
800
1200
2000
500
2000
3000
5000

n
1000
2000
5000
8000
1000
2000
5000
8000
1000
2000
5000
8000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000

cn
0.005
0.005
0.005
0.005
0.01
0.01
0.01
0.01
0.05
0.05
0.05
0.05
0.1
0.1
0.1
0.1
0.5
0.5
0.5
0.5
0.9
0.9
0.9
0.9
1.5
1.5
1.5
1.5
2
2
2
2
5
5
5
5

φq (I)
q=1 q=3
0.081 0.078
0.059 0.060
0.051 0.054
0.050 0.049
0.072 0.067
0.066 0.059
0.046 0.053
0.056 0.049
0.063 0.056
0.058 0.052
0.056 0.055
0.049 0.047
0.073 0.075
0.053 0.062
0.049 0.043
0.062 0.058
0.066 0.067
0.053 0.052
0.053 0.054
0.052 0.050
0.051 0.055
0.056 0.051
0.058 0.050
0.048 0.053
0.042 0.047
0.048 0.054
0.056 0.055
0.051 0.052
0.051 0.051
0.047 0.056
0.055 0.050
0.054 0.053
0.063 0.053
0.049 0.054
0.052 0.056
0.052 0.052

Uq (I)
q=1 q=3
0.065 0.049
0.052 0.050
0.047 0.043
0.054 0.036
0.057 0.048
0.050 0.040
0.045 0.044
0.051 0.046
0.051 0.048
0.061 0.047
0.050 0.047
0.047 0.040
0.062 0.061
0.050 0.043
0.055 0.047
0.053 0.045
0.066 0.050
0.051 0.045
0.044 0.046
0.045 0.051
0.050 0.057
0.050 0.050
0.061 0.049
0.058 0.046
0.065 0.064
0.049 0.040
0.040 0.046
0.041 0.043
0.057 0.049
0.052 0.049
0.051 0.052
0.050 0.051
0.049 0.044
0.042 0.044
0.052 0.042
0.047 0.049

φq (II)
q=1 q=3
0.081 0.074
0.058 0.055
0.054 0.059
0.062 0.055
0.063 0.060
0.057 0.058
0.048 0.053
0.046 0.046
0.046 0.058
0.048 0.052
0.046 0.048
0.055 0.042
0.072 0.074
0.051 0.059
0.052 0.053
0.051 0.054
0.051 0.059
0.048 0.045
0.040 0.048
0.041 0.045
0.048 0.051
0.047 0.048
0.037 0.046
0.045 0.045
0.039 0.048
0.036 0.045
0.039 0.043
0.042 0.049
0.045 0.045
0.043 0.052
0.043 0.052
0.048 0.048
0.051 0.052
0.040 0.048
0.034 0.044
0.043 0.045

Uq (II)
q=1 q=3
0.090 0.081
0.068 0.067
0.055 0.050
0.057 0.051
0.068 0.052
0.065 0.049
0.055 0.044
0.052 0.050
0.055 0.058
0.049 0.047
0.047 0.037
0.046 0.044
0.079 0.083
0.056 0.055
0.053 0.049
0.054 0.044
0.069 0.070
0.053 0.059
0.043 0.038
0.048 0.038
0.069 0.078
0.047 0.047
0.054 0.057
0.048 0.045
0.059 0.070
0.048 0.049
0.048 0.046
0.049 0.047
0.067 0.066
0.046 0.048
0.053 0.040
0.049 0.040
0.051 0.068
0.055 0.048
0.049 0.051
0.049 0.057

35

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

Table 2
Test Power under Scenario (III) and (IV)
p
5
10
25
40
10
20
50
80
50
100
250
400
10
40
60
100
50
200
300
500
90
360
540
900
150
600
900
1500
200
800
1200
2000
500
2000
3000
5000

n
1000
2000
5000
8000
1000
2000
5000
8000
1000
2000
5000
8000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000
100
400
600
1000

cn
0.005
0.005
0.005
0.005
0.01
0.01
0.01
0.01
0.05
0.05
0.05
0.05
0.1
0.1
0.1
0.1
0.5
0.5
0.5
0.5
0.9
0.9
0.9
0.9
1.5
1.5
1.5
1.5
2
2
2
2
5
5
5
5

φq (III)
q=1
q=3
0.999
0.990
1
1
1
1
1
1
1
0.998
1
1
1
1
1
1
1
0.9995
1
1
1
1
1
1
0.335
0.213
0.983
0.764
1
0.973
1
1
0.416
0.245
1
0.957
1
1
1
1
0.506
0.311
1
0.991
1
1
1
1
0.607
0.365
1
1
1
1
1
1
0.694
0.428
1
1
1
1
1
1
0.9375 0.704
1
1
1
1
1
1

Uq (III)
q=1 q=3
0.807 0.635
0.9995 0.986
1
1
1
1
0.824 0.647
1
0.992
1
1
1
1
0.859 0.681
1
0.993
1
1
1
1
0.085 0.077
0.294 0.207
0.482 0.338
0.851 0.651
0.093 0.093
0.287 0.190
0.497 0.327
0.878 0.669
0.090 0.079
0.323 0.208
0.527 0.338
0.897 0.655
0.089 0.078
0.337 0.203
0.573 0.346
0.922 0.674
0.092 0.082
0.351 0.207
0.609 0.351
0.923 0.657
0.102 0.082
0.452 0.190
0.734 0.353
0.986 0.659

φq (IV)
q=1
q=3
0.999
0.986
1
1
1
1
1
1
1
0.997
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0.3055 0.1955
0.968 0.6935
1
0.94
1
1
0.3075 0.1805
0.9995 0.8585
1
0.995
1
1
0.3785 0.217
1
0.952
1
0.9995
1
1
0.4545 0.2745
1
0.992
1
1
1
1
0.5165 0.2985
1
0.9985
1
1
1
1
0.809
0.519
1
1
1
1
1
1

Uq (IV)
q=1 q=3
0.761 0.633
0.999 0.984
1
1
1
1
0.79 0.636
1
0.989
1
1
1
1
0.835 0.660
1
0.996
1
1
1
1
0.129 0.128
0.283 0.225
0.483 0.331
0.852 0.652
0.094 0.091
0.310 0.223
0.491 0.322
0.886 0.646
0.096 0.097
0.315 0.214
0.504 0.331
0.894 0.671
0.102 0.098
0.319 0.212
0.576 0.326
0.906 0.656
0.108 0.098
0.352 0.204
0.592 0.342
0.933 0.656
0.108 0.097
0.430 0.202
0.712 0.329
0.984 0.679

36

Weiming Li, Zeng Li and Jianfeng Yao/Joint CLT for linear spectral statistics

37

Table 3
Size and power comparison of permutation test and φq for Gaussian distribution under Scenario (III),
nominal level α = 5%, testing size is shown in the first row of each (p, n) configuration block.
Permutation

φq

p

n

cn

a

q=1

q=3

q=1

q=3

150

300

0.5

0

0.056

0.062

0.052

0.054

150

300

0.5

0.05

0.360

0.256

0.254

0.146

150

300

0.5

0.09

0.992

0.948

0.934

0.694

150

300

0.5

0.1

1

0.988

0.976

0.820

270

300

0.9

0

0.068

0.074

0.060

0.050

270

300

0.9

0.05

0.510

0.408

0.376

0.216

270

300

0.9

0.09

1

1

0.990

0.820

270

300

0.9

0.1

1

1

0.998

0.914

600

300

2

0

0.056

0.066

0.062

0.062

600

300

2

0.05

0.808

0.708

0.536

0.294

600

300

2

0.09

1

1

1

0.974

600

300

2

0.1

1

1

1

0.994

Table 4
Size and power comparison of permutation test and φq for Gamma distribution under Scenario (IV),
nominal level α = 5%, testing size is shown in the first row of each (p, n) configuration block.
Permutation

φq

p

n

cn

a

q=1

q=3

q=1

q=3

150

300

0.5

0

0.052

0.044

0.050

0.058

150

300

0.5

0.05

0.344

0.276

0.228

0.142

150

300

0.5

0.09

0.996

0.962

0.892

0.534

150

300

0.5

0.1

1

0.988

0.956

0.692

270

300

0.9

0

0.058

0.048

0.046

0.052

270

300

0.9

0.05

0.470

0.372

0.244

0.144

270

300

0.9

0.09

1

1

0.960

0.662

270

300

0.9

0.1

1

1

0.988

0.802

600

300

2

0

0.058

0.054

0.042

0.056

600

300

2

0.05

0.766

0.692

0.366

0.202

600

300

2

0.09

1

1

0.998

0.884

600

300

2

0.1

1

1

1

0.944

