---
abstract: 'This paper is to prove the asymptotic normality of a statistic for detecting the existence of heteroscedasticity for linear regression models without assuming randomness of covariates when the sample size $n$ tends to infinity and the number of covariates $p$ is either fixed or tends to infinity. Moreover our approach indicates that its asymptotic normality holds even without homoscedasticity.'
address:
- 'KLASMOE and School of Mathematics and Statistics, Northeast Normal University, Changchun, P.R.C., 130024.'
- 'School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore, 637371'
- 'KLASMOE and School of Mathematics and Statistics, Northeast Normal University, Changchun, P.R.C., 130024.'
author:
- Zhidong Bai
- Guangming Pan
- Yanqing Yin
title: 'Homoscedasticity tests for both low and high-dimensional fixed design regressions'
---

[ emptyifnotempty[@addto@macro]{} @addto@macro]{} [^1]

[ emptyifnotempty[@addto@macro]{} @addto@macro]{} [^2]

[ emptyifnotempty[Corresponding author]{}[@addto@macro]{} @addto@macro]{} [^3]

Introduction
============

A brief review of homoscedasticity test
---------------------------------------

Consider the classical multivariate linear regression model of $p$ covariates $$\begin{aligned}
  y_i={\mathbf}x_i{\mathbf}\beta+{\mathbf}\varepsilon_i,\ \ \ \ \ i=1,2,\cdots,n,\end{aligned}$$ where $y_i$ is the response variable, ${\mathbf}x_i=(x_{i,1},x_{i,2},\cdots,x_{i,p})$ is the $p$-dimensional covariates, $\beta={\left(}\beta_1,\beta_2,\cdots,\beta_p{\right)}'$ is the $p$ dimensional regression coefficient vector and $\varepsilon_i$ is the independent random errors obey the same distribution with zero mean and variance $\sigma_i^2$. In most applications of the linear regression models the homoscedasticity is a very important assumption. Without it, the loss in efficiency in using ordinary least squares (OLS) may be substantial and even worse, the biases in estimated standard errors may lead to invalid inferences. Thus, it is very important to examine the homoscedasticity. Formally, we need to test the hypothesis $$\label{a1}
H_0: \ \sigma_1^2=\sigma_2^2=\cdots=\sigma_n^2=\sigma^2,$$ where $\sigma^2$ is a positive constant.

In the literature there are a lot of work considering this hypothesis test when the dimension $p$ is fixed. Indeed, many popular tests have been proposed. For example Breusch and Pagan [@breusch1979simple] and White [@white1980heteroskedasticity] proposed statistics to investigate the relationship between the estimated errors and the covariates in economics. While in statistics, Dette and Munk [@dette1998estimating], Glejser [@glejser1969new], Harrison and McCabe [@harrison1979test], Cook and Weisberg [@cook1983diagnostics], Azzalini and Bowman[@azzalini1993use] proposed nonparametric statistics to conduct the hypothesis. One may refer to Li and Yao [@li2015homoscedasticity] for more details in this regard. The development of computer science makes it possible for people to collect and deal with high-dimensional data. As a consequence, high-dimensional linear regression problems are becoming more and more common due to widely available covariates. Note that the above mentioned tests are all developed under the low-dimensional framework when the dimension $p$ is fixed and the sample size $n$ tends to infinity. In Li and Yao’s paper, they proposed two test statistics in the high dimensional setting by using the regression residuals. The first statistic uses the idea of likelihood ratio and the second one uses the idea that “the departure of a sequence of numbers from a constant can be efficiently assessed by its coefficient of variation", which is closely related to John’s idea [@john1971some]. By assuming that the distribution of the covariates is ${\mathbf}N({\mathbf}0, {\mathbf}I_p)$ and that the error obey the normal distribution, the “coefficient of variation" statistic turns out to be a function of residuals. But its asymptotic distribution missed some part as indicated from the proof of Lemma 1 in [@li2015homoscedasticity] even in the random design. The aim of this paper is to establish central limit theorem for the “coefficient of variation" statistic without assuming randomness of the covariates by using the information in the projection matrix (the hat matrix). This ensures that the test works when the design matrix is both fixed and random. More importantly we prove that the asymptotic normality of this statistics holds even without homoscedasticity. That assures a high power of this test.

The structure of this paper is as follows. Section 2 is to give our main theorem and some simulation results, as well as two real data analysis. Some calculations and the proof of the asymptotic normality are presented in Section 3.

Main Theorem, Simulation Results and Real Data Analysis
=======================================================

The Main Theorem
----------------

Suppose that the parameter vector $\beta$ is estimated by the OLS estimator $$\hat{ {\beta}}={\left(}{\mathbf}X'{\mathbf}X{\right)}^{-1}{\mathbf}X'{\mathbf}Y.$$ Denote the residuals by $$\hat{{\mathbf}\varepsilon}={\left(}\hat{{\mathbf}\varepsilon_1},\hat{{\mathbf}\varepsilon_2},\cdots,\hat{{\mathbf}\varepsilon_n}{\right)}'={\mathbf}Y-{\mathbf}X\hat \beta={\mathbf}P\varepsilon,$$ with ${\mathbf}P=(p_{ij})_{n\times n}={\mathbf}I_n-{\mathbf}X({\mathbf}X'{\mathbf}X)^{-1}{\mathbf}X'$ and $\varepsilon={\left(}\varepsilon_1,\varepsilon_2,\cdots,\varepsilon_n{\right)}'$. Let ${\mathbf}D$ be an $n\times n$ diagonal matrix with its $i$-th diagonal entry being $\sigma_i$, set ${\mathbf}A=(a_{ij})_{n\times n}={\mathbf}P{\mathbf}D$ and let $\xi={\left(}\xi_1,\xi_2,\cdots,\xi_n{\right)}'$ stand for a standard $n$ dimensional random vector whose entries obey the same distribution with $\varepsilon$. It follows that the distribution of $\hat \varepsilon$ is the same as that of ${\mathbf}A\xi$. In the following, we use ${{\rm Diag}}{\left(}{\mathbf}B{\right)}={\left(}b_{1,1},b_{2,2},\cdots,b_{n,n}{\right)}'$ to stand for the vector formed by the diagonal entries of ${\mathbf}B$ and ${{\rm Diag}}'{\left(}{\mathbf}B{\right)}$ as its transpose, use ${\mathbf}D_{{\mathbf}B}$ stand for the diagonal matrix of ${\mathbf}B$, and use ${\mathbf}1$ stand for the vector ${\left(}1,1,\cdots,1{\right)}'$.

Consider the following statistic $$\label{a2}
{\mathbf}T=\frac{\sum_{i=1}^n{\left(}\hat{\varepsilon_i}^2-\frac{1}{n}\sum_{i=1}^n\hat{\varepsilon_i}^2{\right)}^2}{\frac{1}{n}{\left(}\sum_{i=1}^n\hat{\varepsilon_i}^2{\right)}^2}.$$ We below use ${\mathbf}A {\circ}{\mathbf}B$ to denote the Hadamard product of two matrices ${\mathbf}A$ and ${\mathbf}B$ and use ${\mathbf}A ^{{\circ}k}$ to denote the Hadamard product of k ${\mathbf}A$.

\[th1\] Under the condition that the distribution of $\varepsilon_1$ is symmetric, ${{\rm E}}|\varepsilon_1|^8\leq \infty$ and $p/n\to y\in [0,1)$ as $n\to \infty$, we have $$\frac{{\mathbf}T-a}{\sqrt b}\stackrel{d}{\longrightarrow}{\mathbf}N(0,1)$$

where $a$, $b$ are determined by $n$, $p$ and ${\mathbf}A$. Under $H_0$, we further have $$a={\left(}\frac{n{\left(}3{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^2{\right)}}{{\left(}{\left(}n-p{\right)}^2+2{\left(}n-p{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P){\right)}}-1{\right)},\ b=\Delta'\Theta\Delta,$$ where $$\Delta'=(\frac{n}{{\left(}{\left(}n-p{\right)}^2+2{\left(}n-p{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P){\right)}},-\frac{n^2{\left(}3{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^2{\right)}}{{{\left(}{\left(}n-p{\right)}^2+2{\left(}n-p{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P){\right)}}^2})$$ and $$\Theta=\left(
 \begin{array}{cc}
  \Theta_{11} & \Theta_{12} \\
   \Theta_{21} & \Theta_{22} \\
   \end{array}
   \right),$$ where $$\begin{aligned}
\Theta_{11}=&72{{\rm Diag}}'({\mathbf}P) {\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}{{\rm Diag}}({\mathbf}P)+24{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}^2\\\notag
&+\nu_4{\left(}96{{\rm tr}}{\mathbf}P {\mathbf}{D_P} {\mathbf}P {\mathbf}P^{{\circ}3}+72{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^3+36{{\rm Diag}}'({\mathbf}P) {\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}^2{{\rm Diag}}({\mathbf}P) {\right)}\\\notag
&+\nu^2_4{\left(}18{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^4+16{{\rm tr}}({\mathbf}P^{{\circ}3}{\mathbf}P)^2){\right)}\\\notag
&+\nu_6{\left(}12{{\rm tr}}{\left(}{\left(}{\mathbf}P{\mathbf}D_{{\mathbf}P}{\mathbf}P {\right)}{\circ}{\left(}{\mathbf}P^{{\circ}2}{\mathbf}P^{{\circ}2}{\right)}{\right)}+16{{\rm tr}}{\mathbf}P {\mathbf}P^{{\circ}3}{\mathbf}P^{{\circ}3}{\right)}+\nu_8{\mathbf}1'({\mathbf}P^{{\circ}4}{\mathbf}P^{{\circ}4}){\mathbf}1,\end{aligned}$$

$$\begin{aligned}
\Theta_{22}=\frac{8{\left(}n-p{\right)}^3+4\nu_4{\left(}n-p{\right)}^2{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)}{n^2},\end{aligned}$$

$$\begin{aligned}
&\Theta_{12}=\Theta_{21}\\\notag
=&\frac{{\left(}n-p{\right)}}{n}{\left(}24{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}+16\nu_4{{\rm tr}}({\mathbf}P {\mathbf}P^{{\circ}3})+12\nu_4{{\rm tr}}{\left(}{\left(}{\mathbf}P{\mathbf}D_{{\mathbf}p}{\mathbf}P{\right)}{\circ}{\mathbf}P{\right)}+2\nu_6[{{\rm Diag}}({\mathbf}P)'({\mathbf}P^{{\circ}4}){\mathbf}1]{\right)},\end{aligned}$$

$\nu_4=M_4-3$ , $\nu_6=M_6-15 M_4+30$ and $\nu_8=M_8-28 M_6-35M_4^2+420M_4-630$ are the corresponding cumulants of random variable $\varepsilon_1$.

The existence of the 8-th moment is necessary because it determines the asymptotic variance of the statistic.

The explicit expressions of $a$ and $b$ are given in Theorem \[th1\] under $H_0$. However the explicit expressions of $a$ and $b$ are quite complicated under $H_1$. Nevertheless one may obtain them from (\[e1\])-(\[e2\]) and (\[t10\])-(\[ct12\]) below.

In Li and Yao’s paper, under the condition that the distribution of $\varepsilon$ is normal, they also did some simulations when the design matrices are non-Gaussian. Specifically speaking, they also investigated the test when the entries of design matrices are drawn from gamma distribution $G(2,2)$ and uniform distribution $U(0,1)$ respectively. There is no significant difference in terms of size and power between these two non-normal designs and the normal design. This seems that the proposed test is robust against the form of the distribution of the design matrix. But according to our main theorem, it is not always the case. In our main theorem, one can find that when the error $\varepsilon$ obey the normal distribution, under $H_0$ and given $p$ and $n$, the expectation of the statistics is only determined by ${{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)$. We conduct some simulations to investigate the influence of the distribution of the design matrix on this term when $n=1000$ and $p=200$. The simulation results are presented in table \[table1\].

                                               $N(0,1)$   $G(2,2)$   $U(0,1)$   $F(1,2)$   $exp(N(5,3))/100$
  ------------------------------------------- ---------- ---------- ---------- ---------- -------------------
   ${{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)$    640.3      640.7      640.2      712.5           708.3

  : The value of ${{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)$ corresponding to different design distributions[]{data-label="table1"}

It suggests that even if the entries of the design matrix are drawn from some common distribution, the expectation of the statistics may deviate far from that of the normal case. This will cause a wrong test result. Moreover, even in the normal case, our result is more accurate since we do not use any approximate value in the mean of the statistic $T$.

Let’s take an example to explain why this test works. For convenient, suppose that $\varepsilon_1$ obey the normal distribution. From the calculation in Section \[exp\] we know that the expectation of the statistic ${\mathbf}T$ defined in (\[a2\]) can be represented as $${{\rm E}}{\mathbf}T=\frac{3n\sum_{i=1}^np_{ii}^2\sigma_i^4}{(\sum_{i=1}^np_{ii}\sigma_i^2)^2}-1+o(1).$$ Now assume that $p_{ii}=\frac{n-p}{n}$ for all $i=1,\cdots,n$. Moreover, without loss of generality, suppose that $\sigma_1=\cdots=\sigma_n=1$ under $H_0$ so that we get ${{\rm E}}{\mathbf}T\to 2$ as $n \to \infty$. However, when $\sigma_1=\cdots=\sigma_{[n/2]}=1$ and $\sigma_{[n/2]+1}=\cdots=\sigma_{n}=2$, one may obtain ${{\rm E}}{\mathbf}T\to 3.08$ as $n\to \infty$. Since ${\rm Var}({\mathbf}T)=O(n^{-1})$ this ensures a high power as long as $n$ is large enough.

Some simulation results
-----------------------

We next conduct some simulation results to investigate the performance of our test statistics. Firstly, we consider the condition when the random error obey the normal distribution. Table \[table2\] shows the empirical size compared with Li and Yao’s result in [@li2015homoscedasticity] under four different design distributions. We use $``{\rm CVT}"$ and $``{\rm FCVT}"$ to represent their test and our test respectively. The entries of design matrices are $i.i.d$ random samples generated from $N(0,1)$, $t(1)$ ($t$ distribution with freedom degree 1), $F(3,2)$ ($F$ distribution with parameters 3 and 2) and logarithmic normal distribution respectively. The sample size $n$ is 512 and the dimension of covariates varies from 4 to 384. We also follow [@dette1998testing] and consider the following two models:

Model 1

:   $y_i={\mathbf}x_i{\mathbf}\beta+{\mathbf}\varepsilon_i(1+{\mathbf}x_i {\mathbf}h),\ \ \ \ \ i=1,2,\cdots,n$,\
    where ${\mathbf}h=(1,{\mathbf}0_{(p-1)})$,

Model 2

:   $y_i={\mathbf}x_i{\mathbf}\beta+{\mathbf}\varepsilon_i(1+{\mathbf}x_i {\mathbf}h),\ \ \ \ \ i=1,2,\cdots,n $\
    where ${\mathbf}h=({\mathbf}1_{(p/2)},{\mathbf}0_{(p/2)})$.

Tables \[table3\] and \[table4\] show the empirical power compared with Li and Yao’s results under four different regressors distributions mentioned above.

Then, we consider the condition that the random error obey the two-point distribution. Specifically speaking, we suppose $P(\varepsilon_1=-1)=P(\varepsilon_1=1)=1/2$. Since Li and Yao’s result is unapplicable in this situation, Table \[table5\] just shows the empirical size and empirical power under Model 2 of our test under four different regressors distributions mentioned above.

According to the simulation result, it is showed that when $p/n\to [0,1)$ as $n\to \infty$, our test always has good size and power under all regressors distributions.

         N(0,1)             t(1)             $F(3,2)$            $e^{(N(5,3))}$              
  ----- -------- -------- -------- -------- ---------- -------- ---------------- -------- -- --
    p     FCVT     CVT      FCVT     CVT       FCVT      CVT          FCVT         CVT       
    4    0.0582   0.0531   0.0600   0.0603    0.0594    0.0597       0.0590       0.0594     
   16    0.0621   0.0567   0.0585   0.0805    0.0585    0.0824       0.0595       0.0803     
   64    0.0574   0.0515   0.0605   0.2245    0.0586    0.2312       0.0578       0.2348     
   128   0.0597   0.0551   0.0597   0.5586    0.0568    0.5779       0.0590       0.5934     
   256   0.0551   0.0515   0.0620   0.9868    0.0576    0.9908       0.0595       0.9933     
   384   0.0580   0.0556   0.0595   1.0000    0.0600    1.0000       0.0600       1.0000     

  : empirical size under different distributions[]{data-label="table2"}

         N(0,1)             t(1)             $F(3,2)$            $e^{(N(5,3))}$              
  ----- -------- -------- -------- -------- ---------- -------- ---------------- -------- -- --
    p     FCVT     CVT      FCVT     CVT       FCVT      CVT          FCVT         CVT       
    4    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   16    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   64    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   128   1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   256   1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   384   0.8113   0.8072   0.9875   1.0000    0.9876    1.0000       0.9905       1.0000     

  : empirical power under model 1[]{data-label="table3"}

         N(0,1)             t(1)             $F(3,2)$            $e^{(N(5,3))}$              
  ----- -------- -------- -------- -------- ---------- -------- ---------------- -------- -- --
    p     FCVT     CVT      FCVT     CVT       FCVT      CVT          FCVT         CVT       
    4    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   16    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   64    1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   128   1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   256   1.0000   1.0000   1.0000   1.0000    1.0000    1.0000       1.0000       1.0000     
   384   0.9066   0.9034   0.9799   1.0000    0.9445    1.0000       0.8883       1.0000     

  : empirical power under model 2[]{data-label="table4"}

         N(0,1)             t(1)             $F(3,2)$            $e^{(N(5,3))}$              
  ----- -------- -------- -------- -------- ---------- -------- ---------------- -------- -- --
    p     Size    Power     Size    Power      Size     Power         Size        Power      
    4    0.0695   1.0000   0.0726   1.0000    0.0726    1.0000       0.0664       1.0000     
   16    0.0695   1.0000   0.0638   1.0000    0.0706    1.0000       0.0556       1.0000     
   64    0.0646   1.0000   0.0606   1.0000    0.0649    1.0000       0.0622       1.0000     
   128   0.0617   1.0000   0.0705   1.0000    0.0597    1.0000       0.0630       1.0000     
   256   0.0684   1.0000   0.0685   1.0000    0.0608    1.0000       0.0649       1.0000     
   384   0.0610   0.8529   0.0748   1.0000    0.0758    1.0000       0.0742       1.0000     

  : empirical size and power under different distributions[]{data-label="table5"}

Two Real Rata Analysis
----------------------

### The Death Rate Data Set

In [@mcdonald1973instabilities], the authors fitted a multiple linear regression of the total age adjusted mortality rate on 15 other variables (the average annual precipitation, the average January temperature, the average July temperature, the size of the population older than 65, the number of members per household, the number of years of schooling for persons over 22, the number of households with fully equipped kitchens, the population per square mile, the size of the nonwhite population, the number of office workers, the number of families with an income less than \$3000, the hydrocarbon pollution index, the nitric oxide pollution index, the sulfur dioxide pollution index and the degree of atmospheric moisture). The number of observations is 60. To investigate whether the homoscedasticity assumption in this models is justified, we applied our test and got a p-value of 0.4994, which strongly supported the assumption of constant variability in this model since we use the one side test. The data set is available at <http://people.sc.fsu.edu/~jburkardt/datasets/regression/regression.html>.

### The 30-Year Conventional Mortgage Rate Data Set

The 30-Year Conventional Mortgage Rate data [@Mortgage] contains the economic data information of USA from 01/04/1980 to 02/04/2000 on a weekly basis (1049 samples). The goal is to predict the 30-Year Conventional Mortgage Rate by other 15 features . We used a multiple linear regression to fit this data set and got a good result. The adjusted R-squared is 0.9986, the P value of the overall F-test is 0. Our homoscedasticity test reported a p-value 0.4439.

Proof Of The Main Theorem
=========================

This section is to prove the main theorem. The first step is to establish the asymptotic normality of ${\mathbf}T_1$, ${\mathbf}T_2$ and $\alpha{\mathbf}T_1+\beta{\mathbf}T_2$ with $\alpha^2+\beta^2 \neq 0$ by the moment convergence theorem. Next we will calculate the expectations, variances and covariance of the statistics ${\mathbf}T_1=\sum_{i=1}^n\hat{\varepsilon_i}^4$ and ${\mathbf}T_2=\frac{1}{n}{\left(}\sum_{i=1}^n\hat{\varepsilon_i}^2{\right)}^2$. The main theorem then follows by the delta method. Note that without loss of generality, under $H_0$, we can assume that $\sigma=1$.

The asymptotic normality of the statistics. {#clt}
-------------------------------------------

We start by giving a definition in Graph Theory.

A graph ${\mathbf}G={\left(}{\mathbf}V,{\mathbf}E,{\mathbf}F{\right)}$ is called two-edge connected, if removing any one edge from $G$, the resulting subgraph is still connected.

The next lemma is a fundamental theorem for Graph-Associated Multiple Matrices without the proof. For the details of this theorem, one can refer to the section ${\mathbf}A.4.2$ in [@bai2010spectral].

\[lm2\] Suppose that ${\mathbf}G={\left(}{\mathbf}V,{\mathbf}E, {\mathbf}F{\right)}$ is a two-edge connected graph with $t$ vertices and $k$ edges. Each vertex $i$ corresponds to an integer $m_i \geq 2$ and each edge $e_j$ corresponds to a matrix ${\mathbf}T^{(j)}={\left(}t_{\alpha,\beta}^{(j)}{\right)},\ j=1,\cdots,k$, with consistent dimensions, that is, if $F(e_j)=(f_i(e_j),f_e(e_j))=(g,h),$ then the matrix ${\mathbf}T^{{\left(}j{\right)}}$ has dimensions $m_g\times m_h$. Define ${\mathbf}v=(v_1,v_2,\cdots,v_t)$ and $$\begin{aligned}
   T'=\sum_{{\mathbf}v}\prod_{j=1}^kt_{v_{f_i(e_j)},v_{f_e(e_j)}}^{(j)},
 \end{aligned}$$ where the summation $\sum_{{\mathbf}v}$ is taken for $v_i=1,2,\cdots, m_i, \ i=1,2,\cdots,t.$ Then for any $i\leq t$, we have $$|T'|\leq m_i\prod_{j=1}^k\|{\mathbf}T^{(j)}\|.$$

Let $\mathcal{T}=({\mathbf}T^{(1)},\cdots,{\mathbf}T^{(k)})$ and define $G(\mathcal{T})=(G,\mathcal{T})$ as a Graph-Associated Multiple Matrices. Write $T'=sum(G(\mathcal{T}))$, which is referred to as the summation of the corresponding Graph-Associated Multiple Matrices.

We also need the following truncation lemma

\[lm3\] Suppose that $\xi_n={\left(}\xi_1,\cdots,\xi_n{\right)}$ is an i.i.d sequence with ${{\rm E}}|\xi_1|^r \leq \infty$, then there exists a sequence of positive numbers $(\eta_1,\cdots,\eta_n)$ satisfy that as $n \to \infty$, $\eta_n \to 0$ and $$P(\xi_n\neq\widehat \xi_n,\ {\rm i.o.})=0,$$ where $\widehat \xi_n={\left(}\xi_1 I(|\xi_1|\leq \eta_n n^{1/r}),\cdots,\xi_n I(|\xi_n|\leq \eta_n n^{1/r}){\right)}.$ And the convergence rate of $\eta_n$ can be slower than any preassigned rate.

${{\rm E}}|\xi_1|^r \leq \infty$ indicated that for any $\epsilon>0$, we have $$\sum_{m=1}^\infty 2^{2m}P(|\xi_1|\geq\epsilon2^{2m/r})\leq \infty.$$

Then there exists a sequence of positive numbers $\epsilon=(\epsilon_1,\cdots,\epsilon_m)$ such that $$\sum_{m=1}^\infty 2^{2m}P(|\xi_1|\geq\epsilon_m2^{2m/r})\leq \infty,$$ and $\epsilon_m \to 0$ as $m \to 0$. And the convergence rate of $\epsilon_m$ can be slower than any preassigned rate.

Now, define $\delta_n=2^{1/r}\epsilon_m$ for $2^{2m-1}\leq n\leq 2^{2m}$, we have as $n\to \infty$ $$\begin{aligned}
P(\xi_n\neq\widehat \xi_n,\ {\rm i.o.})\leq &\lim_{k\to \infty}\sum_{m=k}^{\infty}P\Big(\bigcup_{2^{2m-1}\leq n\leq 2^{2m}}\bigcup_{i=1}^n{\left(}|\xi_i|\geq\eta_nn^{1/r}{\right)}\Big)\\\notag
\leq&\lim_{k\to \infty}\sum_{m=k}^{\infty}P\Big(\bigcup_{2^{2m-1}\leq n\leq 2^{2m}}\bigcup_{i=1}^{2^{2m}}{\left(}|\xi_i|\geq\epsilon_m 2^{1/r}2^{\frac{{\left(}2m-1{\right)}}{r}}{\right)}\Big)\\\notag
\leq&\lim_{k\to \infty}\sum_{m=k}^{\infty}P\Big(\bigcup_{2^{2m-1}\leq n\leq 2^{2m}}\bigcup_{i=1}^{2^{2m}}{\left(}|\xi_i|\geq\epsilon_m 2^{{2m}/{r}}{\right)}\Big)\\\notag
=&\lim_{k\to \infty}\sum_{m=k}^{\infty}P\Big(\bigcup_{i=1}^{2^{2m}}{\left(}|\xi_i|\geq\epsilon_m 2^{{2m}/{r}}{\right)}\Big)\\\notag
\leq&\lim_{k\to \infty}\sum_{m=k}^{\infty}2^{2m}P\Big(|\xi_1|\geq\epsilon_m 2^{{2m}/{r}}\Big)=0.\end{aligned}$$

We note that the truncation will neither change the symmetry of the distribution of $\xi_1$ nor change the order of the variance of ${\mathbf}T$.

Now, we come to the proof of the asymptotic normality of the statistics. We below give the proof of the asymptotic normality of $\alpha{\mathbf}T_1+\beta{\mathbf}T_2$ , where $\alpha^2+\beta^2\neq 0$. The asymptotic normality of either ${\mathbf}T_1$ or ${\mathbf}T_2$ is a result of setting $\alpha=0$ or $\beta=0$ respectively.

Denote $\mu_1={{\rm E}}{\mathbf}T_1={{\rm E}}\sum_{i=1}^n\hat \varepsilon_i^4$, $\mu_2={{\rm E}}{\mathbf}T_2={{\rm E}}n^{-1}{\left(}\sum_{i=1}^n\hat \varepsilon_i^2{\right)}^2$ and $S=\sqrt{{\rm {Var}}{\left(}\alpha {\mathbf}T_1+\beta{\mathbf}T_2{\right)}}$. Below is devote to calculating the moments of ${\mathbf}T_0=\frac{\alpha {\mathbf}T_1+\beta {\mathbf}T_2-{\left(}\alpha \mu_1+\beta \mu_2{\right)}}{S}=\frac{\alpha {\left(}{\mathbf}T_1-\mu_1{\right)}+\beta {\left(}{\mathbf}T_2-\mu_2{\right)}}{S}$.

Note that by Lemma \[lm3\], we can assume that $\xi_1$ is truncated at $\eta_n n^{1/8}$. Then we have for large enough $n$ and $l>4$, $$M_{2l}\leq \eta_n M_8{\sqrt n}^{2l/4-1}.$$

Let’s take a look at the random variable $$\begin{aligned}
&\alpha T_1+\beta T_2=\alpha \sum_{i=1}^n{\left(}\sum_{j=1}^n a_{ij}\xi_j{\right)}^4+(n^{-1})\beta {\left(}\sum_{i=1}^n{\left(}\sum_{j=1}^n a_{ij}\xi_j{\right)}^2{\right)}^2\\\notag
=&\alpha \sum_{i,j_1,\cdots,j_4} a_{i,j_1}a_{i,j_2}a_{i,j_3}a_{i,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}+(n^{-1})\beta\sum_{i_1,i_2,j_1,\cdots,j_4} a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_3}a_{i_2,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}\\\notag
=&\alpha \sum_{i,j_1,\cdots,j_4} a_{i,j_1}a_{i,j_2}a_{i,j_3}a_{i,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}+(n^{-1})\beta\sum_{u_1,u_2,v_1,\cdots,v_4} a_{u_1,v_1}a_{u_1,v_2}a_{u_2,v_3}a_{u_2,v_4}\xi_{v_1}\xi_{v_2}\xi_{v_3}\xi_{v_4}.\end{aligned}$$ We next construct two type of graphs for the last two sums.

For given integers $i,j_1,j_2,j_3,j_4\in [1,n]$, draw a graph as follows: draw two parallel lines, called the $I$-line and the $J$-line respectively; plot $i$ on the $I$-line and $j_1,j_2,j_3$ and $j_4$ on the $J$-line; finally, we draw four edges from $i$ to $j_t$, $t=1,2,3,4$ marked with $\textcircled{1}$. Each edge $(i,j_t)$ represents the random variable $a_{i,j_t}\xi_{j_t}$ and the graph $G_1(i,{\mathbf}j)$ represents $\prod_{\rho=1}^{4}a_{i,j_\rho}\xi_{j_\rho}$. For any given integer $k_1$, we draw $k_1$ such graphs between the $I$-line and the $J$-line denoted by $G_1(\tau)=G_1(i_\tau,{\mathbf}j_\tau)$, and write $G_{(1,k_1)}=\cup_{\tau} G_1(\tau)$.

For given integers $u_1,u_2,v_1,v_2,v_3,v_4\in [1,n]$, draw a graph as follows: plot $u_1$ and $u_2$ on the $I$-line and $v_1,v_2,v_3$ and $v_4$ on the $J$-line; then, we draw two edges from $u_1$ to $v_1$ and $v_2$ marked with $\textcircled{2}$ , draw two edges from $u_2$ to $v_3$ and $v_4$ marked with $\textcircled{2}$. Each edge $(u_l,v_t)$ represents the random variable $a_{u_l,v_t}\xi_{v_t}$ and the graph $G_2({\mathbf}u,{\mathbf}v)$ represents $a_{u_1,v_1}a_{u_1,v_2}a_{u_2,v_3}a_{u_2,v_4}$. For any given integer $k_2$, we draw $k_2$ such graphs between the $I$-line and the $J$-line denoted by $G_2(\psi)=G_2({\mathbf}u_\psi,{\mathbf}v_\psi)$, and write $G_{(2,k_2)}=\cup_{\psi} G_2(\psi)$, $G_{k}=G_{(1,k_1)}\cup G_{(2,k_2)}$. Then the $k$-th order moment of ${\mathbf}T_0$ is $$\begin{aligned}
M_k'=&S^{-k}\sum_{k_1+k_2=k}{k\choose k_1}\alpha^{k_1}\beta^{k_2}\sum_{\substack{\{i_1,{\mathbf}j_1,\cdots,i_{k_1},{\mathbf}j_{k_1}\} \\ \{{\mathbf}u_1,{\mathbf}v_1,\cdots,{\mathbf}u_{k_2},{\mathbf}v_{k_2}\}}}\\
&n^{-k_2}{{\rm E}}\Big[\prod_{\tau=1}^{k_1}[G_1(i_\tau,{\mathbf}j_\tau)-{{\rm E}}(G_1(i_\tau,{\mathbf}j_\tau))]\prod_{\phi=1}^{k_2}[G_2({\mathbf}u_\psi,{\mathbf}v_\psi)-{{\rm E}}(G_2({\mathbf}u_\psi,{\mathbf}v_\psi))]\Big].\end{aligned}$$

We first consider a graph $G_k$ for the given set of integers $k_1,k_2$, $i_1, {\mathbf}j_1,\cdots,i_{k_1},{\mathbf}j_{k_1}$ and ${\mathbf}u_1,{\mathbf}v_1,\cdots,{\mathbf}u_{k_2},{\mathbf}v_{k_2}$. We have the following simple observations: Firstly, if $G_k$ contains a $j$ vertex of odd degree, then the term is zero because odd-ordered moments of random variable $\xi_j$ are 0. Secondly, if there is a subgraph $G_1(\tau)$ or $G_2(\psi)$ that does not have an $j$ vertex coinciding with any $j$ vertices of other subgraphs, the term is also 0 because $G_1(\tau)$ or $G_2(\psi)$ is independent of the remainder subgraphs.

Then, upon these two observations, we split the summation of non-zero terms in $M_k'$ into a sum of partial sums in accordance of isomorphic classes (two graphs are called isomorphic if one can be obtained from the other by a permutation of $(1,2,\cdots,n)$, and all the graphs are classified into isomorphic classes. For convenience, we shall choose one graph from an isomorphic class as the canonical graph of that class). That is, we may write $$\begin{aligned}
M_k'=&S^{-k}\sum_{k_1+k_2=k}{k\choose k_1}\alpha^{k_1}\beta^{k_2}n^{-k_2}\sum_{G_k'}M_{G_k'},\end{aligned}$$ where $$M_{G_k'}=\sum_{G_k\in G_k'}{{\rm E}}G_k.$$ Here $G_k'$ is a canonical graph and $\sum_{G_k\in G_k'}$ denotes the summation for all graphs $G_k$ isomorphic to $G_k'$.

In that follows, we need the fact that the variances of ${\mathbf}T_1$ and ${\mathbf}T_2$ and their covariance are all of order n. This will be proved in Section \[var\].

Since all of the vertices in the non-zero canonical graphs have even degrees, every connected component of them is a circle, of course a two-edge connected graph. For a given isomorphic class with canonical graph $G_k'$, denote by $c_{G_k'}$ the number of connected components of the canonical graph $G_k'$. For every connected component $G_0$ that has $l$ non-coincident $J$-vertices with degrees $d_1,\cdots,d_l$, let $d'=\max\{d_1-8,\cdots,d_l-8,0\}$, denote $\mathcal{T}=(\underbrace{{\mathbf}A,\cdots,{\mathbf}A}_{\sum_{t=1}^l d_t})$ and define $G_0(\mathcal{T})=(G_0,\mathcal{T})$ as a Graph-Associated Multiple Matrices. By Lemma \[lm2\] we then conclude that the contribution of this canonical class is at most ${\left(}\prod_{t=1}^l M_{d_t}{\right)}sum(G(\mathcal{T}))=O(\eta_n^{d'}n\sqrt n^{d'/4})$. Noticing that $\eta_n \to 0$, if $c_{G_k'}$ is less than $k/2+k_2$, then the contribution of this canonical class is negligible because $S^k\asymp n^{k/2}$ and $M_{G_k'}$ in $M_k'$ has a factor of $n^{-k_2}$. However one can see that $c_{G_k'}$ is at most $[k/2]+k_2$ for every ${G_k'}$ by the argument above and noticing that every $G_2(\bullet)$ has two $i$ vertices. Therefore, $M_k'\to 0$ if $k$ is odd.

Now we consider the limit of $M_k'$ when $k=2s$. We shall say that [*the given set of integers $i_1, {\mathbf}j_1,\cdots,i_{k_1},{\mathbf}j_{k_1}$ and ${\mathbf}u_1,{\mathbf}v_1,\cdots,{\mathbf}u_{k_2},{\mathbf}v_{k_2}$ (or equivalent the graph $G_k$) satisfies the condition $c(s_1,s_2,s_3)$ if in the graph $G_k$ plotted by this set of integers there are $2s_1$ $G_1{{\left(}\bullet{\right)}}$ connected pairwisely, $2s_2$ $G_2{{\left(}\bullet{\right)}}$ connected pairwisely and $s_3$ $G_1{{\left(}\bullet{\right)}}$ connected with $s_3$ $G_2{{\left(}\bullet{\right)}}$, where $2s_1+s_3=k_1$, $2s_2+s_3=k_2$ and $s_1+s_2+s_3=s$, say $G_1{{\left(}2\tau-1{\right)}}$ connects $G_1{{\left(}2\tau{\right)}}$, $\tau=1,2,\cdots,s_1$, $G_2{{\left(}2\psi-1{\right)}}$ connects $G_1{{\left(}2\psi{\right)}}$, $\psi=1,2,\cdots,s_2$ and $G_1{{\left(}2s_1+\varphi{\right)}}$ connects $G_2{{\left(}2s_2+\varphi{\right)}}$, $\varphi=1,2,\cdots,s_3$, and there are no other connections between subgraphs.*]{} Then, for any $G_k$ satisfying $c(s_1,s_2,s_3)$, we have $$\begin{aligned}
  {{\rm E}}G_k=&\prod_{\tau=1}^{s_1}{{\rm E}}[(G_1{{\left(}2\tau-1{\right)}}-{{\rm E}}(G_1{{\left(}2\tau-1{\right)}}))(G_1{{\left(}2\tau{\right)}}-{{\rm E}}(G_1{\left(}{2\tau}{\right)}))]\times\\\notag
           &\prod_{\psi=1}^{s_2}{{\rm E}}[(G_2{{\left(}2\psi-1{\right)}}-{{\rm E}}(G_2{{\left(}2\psi-1{\right)}}))(G_2{{\left(}2\psi{\right)}}-{{\rm E}}(G_2{\left(}{2\psi}{\right)}))]\times\\\notag
           &\prod_{\varphi=1}^{s_3}{{\rm E}}[(G_1{{\left(}2s_1+\varphi{\right)}}-{{\rm E}}(G_1{{\left(}2s_1+\varphi{\right)}}))(G_2{{\left(}2s_2+\varphi{\right)}}-{{\rm E}}(G_2{\left(}{2s_2+\varphi}{\right)}))].\end{aligned}$$

Now, we compare $$\begin{aligned}
&n^{-k_2}\sum_{G_k\in c(s_1,s_2,s_3)} {{\rm E}}G_k\\\notag
=&n^{-k_2}\sum_{G_k\in c(s_1,s_2,s_3)}\prod_{\tau=1}^{s_1}{{\rm E}}[(G_1{{\left(}2\tau-1{\right)}}-{{\rm E}}(G_1{{\left(}2\tau-1{\right)}}))(G_1{{\left(}2\tau{\right)}}-{{\rm E}}(G_1{\left(}{2\tau}{\right)}))]\times\\\notag
           &\prod_{\psi=1}^{s_2}{{\rm E}}[(G_2{{\left(}2\psi-1{\right)}}-{{\rm E}}(G_2{{\left(}2\psi-1{\right)}}))(G_2{{\left(}2\psi{\right)}}-{{\rm E}}(G_2{\left(}{2\psi}{\right)}))]\times \\\notag
           &\prod_{\varphi=1}^{s_3}{{\rm E}}[(G_1{{\left(}2s_1+\varphi{\right)}}-{{\rm E}}(G_1{{\left(}2s_1+\varphi{\right)}}))(G_2{{\left(}2s_2+\varphi{\right)}}-{{\rm E}}(G_2{\left(}{2s_2+\varphi}{\right)}))],
           \end{aligned}$$ with $$\begin{aligned}
&{\left(}{{\rm E}}{\left(}{\mathbf}T_1-\mu_1{\right)}^2{\right)}^{s_1}{\left(}{{\rm E}}{\left(}{\mathbf}T_2-\mu_2{\right)}^2{\right)}^{s_2}{\left(}{{\rm E}}{\left(}{\mathbf}T_1-\mu_1{\right)}{\left(}{\mathbf}T_2-\mu_2{\right)}{\right)}^{s_3}\\\notag
=&n^{-k_2}\sum_{G_k}\prod_{\tau=1}^{s_1}{{\rm E}}[(G_1{{\left(}2\tau-1{\right)}}-{{\rm E}}(G_1{{\left(}2\tau-1{\right)}}))(G_1{{\left(}2\tau{\right)}}-{{\rm E}}(G_1{\left(}{2\tau}{\right)}))]\times\\\notag
           &\prod_{\psi=1}^{s_2}{{\rm E}}[(G_2{{\left(}2\psi-1{\right)}}-{{\rm E}}(G_2{{\left(}2\psi-1{\right)}}))(G_2{{\left(}2\psi{\right)}}-{{\rm E}}(G_2{\left(}{2\psi}{\right)}))]\times \\\notag
           &\prod_{\varphi=1}^{s_3}{{\rm E}}[(G_1{{\left(}2s_1+\varphi{\right)}}-{{\rm E}}(G_1{{\left(}2s_1+\varphi{\right)}}))(G_2{{\left(}2s_2+\varphi{\right)}}-{{\rm E}}(G_2{\left(}{2s_2+\varphi}{\right)}))],\end{aligned}$$ where $\sum_{G_k\in c(s_1,s_2,s_3)}$ stands for the summation running over all graph $G_k$ satisfying the condition $c(s_1,s_2,s_3)$.

If $G_k$ satisfies the two observations mentioned before, then ${{\rm E}}G_k=0$, which does not appear in both expressions; if $G_k$ satisfies the condition $c(s_1,s_2,s_3)$, then the two expressions both contain ${{\rm E}}G_k$. Therefore, the second expression contains more terms that $G_k$ have more connections among subgraphs than the condition $c(s_1,s_2,s_3)$. Therefore, by Lemma \[lm2\], $${\left(}{{\rm E}}{\left(}{\mathbf}T_1-\mu_1{\right)}^2{\right)}^{s_1}{\left(}{{\rm E}}{\left(}{\mathbf}T_2-\mu_2{\right)}^2{\right)}^{s_2}{\left(}{{\rm E}}{\left(}{\mathbf}T_1-\mu_1{\right)}{\left(}{\mathbf}T_2-\mu_2{\right)}{\right)}^{s_3}=n^{-k_2}\sum_{G_k\in c(s_1,s_2,s_3)} {{\rm E}}G_k+o(S^k).
\label{map1}$$

If $G_k\in G_k'$ with $c_{G_k'}=s+k_2$, for any nonnegative integers $s_1,s_2,s_3$ satisfying $k_1=2s_1+s_3$, $k_2=2s_2+s_3$ and $s_1+s_2+s_3=s$, we have ${k_1\choose s_3}{k_2 \choose s_3}(2s_1-1)!!(2s_2-1)!!s_3!$ ways to pairing the subgraphs satisfying the condition $c(s_1,s_2,s_3)$. By (\[map1\]), we then have $$\begin{aligned}
&&\sum_{c_{G_k'}=s+k_2}n^{-k_2}{{\rm E}}G_k+o(S^k)\\
&=&\sum_{s_1+s_2+s_3=s\atop 2s_1+s_3=k_1, 2s_2+s_3=k_2}{k_1\choose s_3}{k_2 \choose s_3}(2s_1-1)!!(2s_2-1)!!s_3!
(Var({\mathbf}T_1))^{s_1}(Var({\mathbf}T_2))^{s_2}(Cov({\mathbf}T_1,{\mathbf}T_2))^{s_3}\end{aligned}$$

It follows that $$\begin{aligned}
M_k'=&S^{-k}\sum_{k_1+k_2=k}{k\choose k_1}\alpha^{k_1}\beta^{k_2}n^{-k_2}\sum_{c_{G_k'}=s+k_2}{{\rm E}}G_k+o(1)\\
=&\Big(S^{-2s}\sum_{k_1=0}^{2s}\sum_{s_3=0}^{\min\{k_1,k_2\}}{2s\choose k_1}{k_1 \choose s_3}{k_2 \choose s_3}{\left(}2s_1-1{\right)}!!{\left(}2s_2-1{\right)}!!s_3!\\
&{\left(}\alpha^2Var({\mathbf}T_1){\right)}^{s_1}{\left(}\beta^2Var({\mathbf}T_2){\right)}^{s_2}{\left(}\alpha\beta Cov({\mathbf}T_1,{\mathbf}T_2){\right)}^{s_3}\Big)+o(1)\\
=&\Big(S^{-2s}\sum_{s_1+s_2+s_3=s}{2s\choose 2s_1+s_3}{2s_1+s_3 \choose s_3}{2s_2+s_3 \choose s_3}{\left(}2s_1-1{\right)}!!{\left(}2s_2-1{\right)}!!s_3!\\
&{\left(}\alpha^2Var({\mathbf}T_1){\right)}^{s_1}{\left(}\beta^2Var({\mathbf}T_2){\right)}^{s_2}{\left(}\alpha\beta Cov({\mathbf}T_1,{\mathbf}T_2){\right)}^{s_3}\Big)+o(1)\\
=&\Big(S^{-2s}\sum_{s_1+s_2+s_3=s}\frac{(2s)!(2s_1+s_3)!(2s_2+s_3)!}{{\left(}2s_1+s_3{\right)}!(2s_2+s_3)!s_3!(2s_1)!s_3!(2s_2)!}{\left(}2s_1-1{\right)}!!{\left(}2s_2-1{\right)}!!s_3!\\
&{\left(}\alpha^2Var({\mathbf}T_1){\right)}^{s_1}{\left(}\beta^2Var({\mathbf}T_2){\right)}^{s_2}{\left(}\alpha\beta Cov({\mathbf}T_1,{\mathbf}T_2){\right)}^{s_3}\Big)+o(1)\\
=&\Big(S^{-2s}\sum_{s_1+s_2+s_3=s}(2s-1)!!\frac{s!}{s_1!s_2!s_3!}\\
&{\left(}\alpha^2Var({\mathbf}T_1){\right)}^{s_1}{\left(}\beta^2Var({\mathbf}T_2){\right)}^{s_2}{\left(}2\alpha\beta Cov({\mathbf}T_1,{\mathbf}T_2){\right)}^{s_3}\Big)+o(1),\end{aligned}$$ which implies that $$M'_{2s}\to (2s-1)!!.$$

Combining the arguments above and the moment convergence theorem we conclude that $$\frac{{\mathbf}T_1-{{\rm E}}{\mathbf}T_1}{\sqrt{{\rm Var}{\mathbf}T_1}}\stackrel{d}{\rightarrow} {\rm N}{\left(}0,1{\right)},\ \frac{{\mathbf}T_2-{{\rm E}}{\mathbf}T_2}{\sqrt{{\rm Var}{\mathbf}T_2}}\stackrel{d}{\rightarrow} {\rm N}{\left(}0,1{\right)}, \ \frac{ {\left(}\alpha {\mathbf}T_1+\beta {\mathbf}T_2{\right)}-{{\rm E}}{\left(}\alpha {\mathbf}T_1+\beta {\mathbf}T_2{\right)}}{\sqrt{{\rm Var}{\left(}\alpha {\mathbf}T_1+\beta {\mathbf}T_2{\right)}}}\stackrel{d}{\rightarrow} {\rm N}{\left(}0,1{\right)},$$ where $\alpha^2+\beta^2\neq 0.$ Let $$\Sigma=\left(
             \begin{array}{cc}
               {\rm {Var}}({\mathbf}T_1) & \rm {Cov}({\mathbf}T_1,{\mathbf}T_2) \\
               \rm {Cov}({\mathbf}T_1,{\mathbf}T_2) & \rm {Var}({\mathbf}T_2) \\
             \end{array}
           \right)
.$$ We conclude that $\Sigma^{-1/2}{\left(}{\mathbf}T_1-{{\rm E}}{\mathbf}T_1,{\mathbf}T_2-{{\rm E}}{\mathbf}T_2{\right)}'$ is asymptotic two dimensional gaussian vector.

The expectation {#exp}
---------------

In the following let ${\mathbf}B={\mathbf}A{\mathbf}A'$. Recall that

$${\mathbf}{T_1}=\sum_{i=1}^n\widehat{\varepsilon_i}^4=\sum_{i=1}^n{\left(}\sum_{j=1}^n a_{i,j}\xi_j{\right)}^4=\sum_{i=1}^n\sum_{j_1,j_2,j_3,j_4}a_{i,j_1}a_{i,j_2}a_{i,j_3}a_{i,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4},$$ $${\mathbf}T_2=n^{-1}{\left(}\sum_{i=1}^n{\left(}\sum_{j=1}^n a_{i,j}\xi_j{\right)}^2{\right)}^2
=n^{-1}\sum_{i_1,i_2}\sum_{j_1,j_2,j_3,j_4}a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_3}a_{i_2,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}.$$ Since all odd moments of $\xi_1,\cdots,\xi_n$ are 0, we know that ${{\rm E}}{\mathbf}T_1$ and ${{\rm E}}{\mathbf}T_2$ are only affected by terms whose multiplicities of distinct values in the sequence $(j_1,\cdots,j_4)$ are all even.

We need to evaluate the mixed moment ${{\rm E}}{\left(}{\mathbf}T_1^{\gamma}{\mathbf}T_2^{\omega}{\right)}$. For simplifying notations particularly in Section \[var\] we introduce the following notations $$\begin{aligned}
&\Omega_{\{\omega_1,\omega_2,\cdots,\omega_s\}}^{{\left(}\gamma_1,\gamma_2,\cdots,\gamma_t{\right)}}[\underbrace{{\left(}\phi_{1,1},\cdots,\phi_{1,s}{\right)},{\left(}\phi_{2,1},\cdots,\phi_{2,s}{\right)},
\cdots,{\left(}\phi_{t,1},\cdots,\phi_{t,s}{\right)}}_{t \ groups}]_0 \\
=&\sum_{i_1,\cdots,i_{t},j_1\neq\cdots\neq j_{s}}\prod_{\tau=1,\cdots, t}\prod_{\rho=1,\cdots, s} a_{i_{\tau},j_{\rho}}^{\phi_{\tau,\rho}},\end{aligned}$$ where $i_1,\cdots,i_{t}$ and $j_1,\cdots,j_{s}$ run over $1,\cdots,n$ and are subject to the restrictions that $j_1,\cdots,j_s$ are distinct; $\sum_{l=1}^t \gamma_l=\sum_{l=1}^s \omega_l=\theta,$ and for any $k=1,\cdots, s$, $\sum_{l=1}^t \phi_{l,k}=\theta$. Intuitively, $t$ is the number of distinct $i$-indices and $s$ that of distinct $j$’s; $\gamma_\tau$ is the multiplicity of the index $i_\tau$ and $\omega_\rho=\sum_{l=1}^t\phi_{l,\rho}$ that of $j_\rho$; $\phi_{\tau,\rho}$ the multiplicity of the factor $a_{i_\tau,j_\rho}$; and $\theta=4(\gamma+\omega)$.

Define $$\begin{aligned}
&\Omega_{\{\omega_1,\omega_2,\cdots,\omega_s\}}^{{\left(}\gamma_1,\gamma_2,\cdots,\gamma_t{\right)}}[\underbrace{{\left(}\phi_{1,1},\cdots,\phi_{1,s}{\right)},{\left(}\phi_{2,1},\cdots,\phi_{2,s}{\right)},
\cdots,{\left(}\phi_{t,1},\cdots,\phi_{t,s}{\right)}}_{t \ groups}] \\
=&\sum_{i_1,\cdots,i_{t},j_1,\cdots, j_{s}}\prod_{\tau=1,\cdots, t}\prod_{\rho=1,\cdots, s} a_{i_{\tau},j_{\rho}}^{\phi_{\tau,\rho}}.\end{aligned}$$ The definition above is similar to that of $$\Omega_{\{\omega_1,\omega_2,\cdots,\omega_s\}}^{{\left(}\gamma_1,\gamma_2,\cdots,\gamma_t{\right)}}[\underbrace{{\left(}\phi_{1,1},\cdots,\phi_{1,s}{\right)},{\left(}\phi_{2,1},\cdots,\phi_{2,s}{\right)},
\cdots,{\left(}\phi_{t,1},\cdots,\phi_{t,s}{\right)}}_{t \ groups}]_0$$ without the restriction that the indices $j_1,\cdots,j_s$ are distinct from each other. To help understand these notations we demonstrate some examples as follows. $$\begin{aligned}
\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,2,0,0),(0,0,2,2)]=\sum_{i_1,i_2,j_1,\cdots,j_4}a_{i_1,j_1}^2a_{i_1,j_2}^2a_{i_2,j_3}^2a_{i_2,j_3}^2,\end{aligned}$$

$$\begin{aligned}
\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]=\sum_{i_1,i_2,j_1,\cdots,j_4}a_{i_1,j_1}^2a_{i_1,j_2}a_{i_1,j_3}a_{i_2,j_2}a_{i_2,j_3}a_{i_2,j_4}^2,\end{aligned}$$

$$\begin{aligned}
\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,2,0,0),(0,0,2,2)]_0=\sum_{i_1,i_2,j_1\neq\cdots\neq j_4}a_{i_1,j_1}^2a_{i_1,j_2}^2a_{i_2,j_3}^2a_{i_2,j_3}^2,\end{aligned}$$

$$\begin{aligned}
\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]_0=\sum_{i_1,i_2,j_1\neq\cdots\neq j_4}a_{i_1,j_1}^2a_{i_1,j_2}a_{i_1,j_3}a_{i_2,j_2}a_{i_2,j_3}a_{i_2,j_4}^2.\end{aligned}$$

We further use $M_k$ to denote the $k$-th order moment of the error random variable. We also use ${\mathbf}C_{n}^k$ to denote the combinatorial number $n \choose k$. We then obtain

$$\begin{aligned}
\label{e1}
&{{\rm E}}{\mathbf}T_1={{\rm E}}\sum_{i=1}^n\sum_{j_1,j_2,j_3,j_4}a_{i,j_1}a_{i,j_2}a_{i,j_3}a_{i,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}\\\notag
=&M_4\Omega_{\{4\}}^{(4)}+M_2^2{\Omega_{\{2,2\}}^{(4)}}_0=M_4\Omega_{\{4\}}^{(4)}+\frac{{\mathbf}C_4^2}{2!}{\Omega_{\{2,2\}}^{(4)}}[{\left(}2,0{\right)},{\left(}0,2{\right)}]_0\\\notag
=&M_4\Omega_{\{4\}}^{(4)}+\frac{{\mathbf}C_4^2}{2!}{\left(}{\Omega_{\{2,2\}}^{(4)}}[{\left(}2,0{\right)},{\left(}0,2{\right)}]-\Omega_{\{4\}}^{(4)}{\right)}\\\notag
=&\frac{{\mathbf}C_4^2}{2!}{\Omega_{\{2,2\}}^{(4)}}[{\left(}2,0{\right)},{\left(}0,2{\right)}]+\nu_4 \Omega_{\{4\}}^{(4)}
=3\sum_i{\left(}\sum_{j}a_{i,j}^2{\right)}^2+\nu_4\sum_{ij}a_{ij}^4\\\notag
=&3\sum_ib_{i,i}^2+\nu_4\sum_{ij}a_{ij}^4=3{{\rm tr}}{\left(}{\mathbf}B{\circ}{\mathbf}B{\right)}+\nu_4{{\rm tr}}({\mathbf}A{\circ}{\mathbf}A)'({\mathbf}A{\circ}{\mathbf}A),\end{aligned}$$

where $\nu_4=M_4-3$ and

$$\begin{aligned}
\label{e2}
&{{\rm E}}{\mathbf}T_2=n^{-1}{{\rm E}}\sum_{i_1,i_2}\sum_{j_1,j_2,j_3,j_4}a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_3}a_{i_2,j_4}\xi_{j_1}\xi_{j_2}\xi_{j_3}\xi_{j_4}\\\notag
=&n^{-1}{\left(}M_4\Omega_{\{4\}}^{(2,2)}+M_2^2{\Omega_{\{2,2\}}^{(2,2)}}_0{\right)}\\\notag
=&n^{-1}{\left(}M_4\Omega_{\{4\}}^{(2,2)}+{\left(}{\Omega_{\{2,2\}}^{(2,2)}}[{\left(}2,0{\right)},{\left(}0,2{\right)}]_0+2\Omega_{\{2,2\}}^{(2,2)}[{\left(}1,1{\right)},{\left(}1,1{\right)}]_0{\right)}{\right)}\\\notag
=&n^{-1}{\left(}M_4\Omega_{\{4\}}^{(2,2)}+{\left(}{\Omega_{\{2,2\}}^{(2,2)}}[{\left(}2,0{\right)},{\left(}0,2{\right)}]+2\Omega_{\{2,2\}}^{(2,2)}[{\left(}1,1{\right)},{\left(}1,1{\right)}]{\right)}-3\Omega_{\{4\}}^{(2,2)}{\right)}\\\notag
=&n^{-1}{\left(}\sum_{i_1,i_2,j_1,j_2}a_{i_1,j_1}^2a_{i_2,j_2}^2+2\sum_{i_1,i_2,j_1,j_2}a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_1}a_{i_2,j_2}
+\nu_4\sum_{i_1,i_2,j}a_{i_1j}^2a_{i_2j}^2{\right)}\\\notag
=&n^{-1}{\left(}{\left(}\sum_{i,j}a_{i,j}^2{\right)}^2+2\sum_{i_1,i_2}{\left(}\sum_{j}a_{i_1,j}a_{i_2,j}{\right)}^2+\nu_4\sum_{i_1,i_2,j}a_{i_1j}^2a_{i_2j}^2{\right)}\\\notag
=&n^{-1}{\left(}{\left(}\sum_{i,j}a_{i,j}^2{\right)}^2+2\sum_{i_1,i_2}b_{i_1,i_2}^2+\nu_4\sum_{j=1}^nb_{jj}^2{\right)}\\ \notag
=&n^{-1}{\left(}{\left(}{{\rm tr}}{\mathbf}B{\right)}^2+2{{\rm tr}}{\mathbf}B^2+\nu_4{{\rm tr}}({\mathbf}B{\circ}{\mathbf}B){\right)}.\end{aligned}$$

The variances and covariance {#var}
----------------------------

We are now in the position to calculate the variances of ${\mathbf}T_1$, ${\mathbf}T_2$ and their covariance.

First, we have

$$\begin{aligned}
\label{t10}
&{\rm Var}( {\mathbf}T_1)={{\rm E}}{\left(}\sum_{i}\widehat{\varepsilon_i}^4-{{\rm E}}{\left(}\sum_{i}\widehat{\varepsilon_i}^4{\right)}{\right)}^2\\\notag
=&\sum_{i_1,i_2,j_1,\cdots,j_8}[{{\rm E}}G(i_1,{\mathbf}j_1)G(i_2,{\mathbf}j_2)-{{\rm E}}G(i_1,{\mathbf}j_1){{\rm E}}G(i_2,{\mathbf}j_2)]\\\notag
=&\Bigg(\Omega_{\{8\}}^{(4,4)}+\Omega_{\{2,6\}_0}^{(4,4)}+\Omega_{\{4,4\}_0}^{(4,4)}+\Omega_{\{2,2,4\}_0}^{(4,4)}+\Omega_{\{2,2,2,2\}_0}^{(4,4)}\Bigg),\\\notag\end{aligned}$$

where the first term comes from the graphs in which the 8 $J$-vertices coincide together; the second term comes from the graphs in which there are 6 $J$-vertices coincident and another two coincident and so on.

Because $G(i_1,{\mathbf}j_1)$ and $G(i_2,{\mathbf}j_2)$ have to connected each other, thus, we have $$\begin{aligned}
\label{t11}
&\Omega_{\{2,2,2,2\}_0}^{(4,4)}\\\notag=&\frac{{\mathbf}C_4^2{\mathbf}C_4^2{\mathbf}C_2^1{\mathbf}C_2^1}{2!}\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]_0+{{\mathbf}C_4^1{\mathbf}C_3^1{\mathbf}C_2^1}\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]_0\\\notag
=&72\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]-4\Omega_{\{2,2,4\}}^{(4,4)}[(2,1,1),(0,1,3)]_0-\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]_0\\\notag
 &-\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]_0-2\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]_0-\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
 &-2\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0-2\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]_0-\Omega_{\{8\}}^{(4,4)}\Big)\\\notag
 &+24\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]-6\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]_0-3\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
 &-4\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0-\Omega_{\{8\}}^{(4,4)}\Big).\\\notag
=&72\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]-4\Omega_{\{2,2,4\}}^{(4,4)}[(2,1,1),(0,1,3)]-\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]\\\notag
 &-\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+2\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]_0+\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
 &+4\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0+4\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]_0+5\Omega_{\{8\}}^{(4,4)}\Big)\\\notag
 &+24\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]-6\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+3\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
 &+8\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0+5\Omega_{\{8\}}^{(4,4)}\Big).\\\notag
=&72\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]-4\Omega_{\{2,2,4\}}^{(4,4)}[(2,1,1),(0,1,3)]-\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]\\\notag
 &-\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+2\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]+\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]\\\notag
 &+4\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]+4\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]-6\Omega_{\{8\}}^{(4,4)}\Big)\\\notag
 &+24\Big(\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]-6\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+3\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]\\\notag
 &+8\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]-6\Omega_{\{8\}}^{(4,4)}\Big).\end{aligned}$$

Likewise we have

$$\begin{aligned}
\label{t12}
&{\Omega_{\{2,2,4\}}^{(4,4)}}_0\\\notag
=&{{\mathbf}C_2^1{\mathbf}C_4^3{\mathbf}C_4^1{\mathbf}C_3^1}M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,2,1),(1,0,3)]_0\\\notag
&+\frac{{\mathbf}C_4^2{\mathbf}C_4^2{\mathbf}C_2^1{\mathbf}C_2^1}{2!}M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]_0+{{\mathbf}C_4^2{\mathbf}C_4^2}(M_4-1)\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]_0\\\notag
=&96M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,2,1),(1,0,3)]_0\\\notag
 &+72M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]_0+36(M_4-1)\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]_0,\\\ \notag
=&96M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,2,1),(1,0,3)]\\\notag
 &+72M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+36(M_4-1)\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]\\\ \notag
 &-96M_4\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]_0-(108 M_4-36)\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
 &-240M_4\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0-(168M_4-72)\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]_0-(204M_4-36)\Omega_{\{8\}}^{(4,4)}\\\notag
 =&96M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,2,1),(1,0,3)]\\\notag
 &+72M_4\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+36(M_4-1)\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]\\\ \notag
 &-96M_4\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]-(108 M_4-36)\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]\\\notag
 &-240M_4\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]-(168M_4-72)\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]+(408M_4-72)\Omega_{\{8\}}^{(4,4)}\end{aligned}$$

$$\begin{aligned}
\label{t13}
\Omega_{\{4,4\}_0}^{(4,4)}=&{{\mathbf}C_2^1{\mathbf}C_4^1{\mathbf}C_4^3}M_4^2\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]_0
+\frac{{\mathbf}C_4^2{\mathbf}C_4^2}{2!}(M_4^2-1)\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]_0\\\notag
=&16M_4^2\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]+18(M_4^2-1)\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]\\\notag
&-(34M_4^2-18)\Omega_{\{8\}}^{(4,4)},\end{aligned}$$

$$\begin{aligned}
\label{t14}
\Omega_{\{2,6\}_0}^{(4,4)}=&{\mathbf}C_2^1{\mathbf}C_4^2(M_6-M_4)\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]_0+{\mathbf}C_4^1{\mathbf}C_4^1M_6\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]_0\\\notag
=&12(M_6-M_4)\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]+16M_6\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]\\\notag
& -(28 M_6-12 M_4)\Omega_{\{8\}}^{(4,4)}.\end{aligned}$$

and $$\begin{aligned}
\label{t15}
\Omega_{\{8\}_0}^{(4,4)}=&(M_8-M_4^2)\Omega_{\{8\}}^{(4,4)}[(4),(4)].\end{aligned}$$

Combining (\[t10\]), (\[t11\]), (\[t12\]), (\[t13\]), (\[t14\]) and (\[t15\]), we obtain

$$\begin{aligned}
\label{vt1}
&{\rm Var}( {\mathbf}T_1)=72\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]+24\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]\\\notag
&+96(M_4-3)\Omega_{\{2,2,4\}}^{(4,4)}[(2,1,1),(0,1,3)]
+36(M_4-3)\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]\\\notag
&+72(M_4-3)\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]+16(M_4^2-6M_4+9)\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]\\\notag
&+18(M_4^2-6M_4+9)\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]+16(M_6-15M_4+30)\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]\\\notag
&+12(M_6-15M_4+30)\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]
+(M_8-28M_6-35M_4^2+420M_4-630)\Omega_{\{8\}}^{(4,4)}[(4),(4)],\end{aligned}$$

where

$$\begin{aligned}
&\Omega_{\{2,2,2,2\}}^{(4,4)}[(2,1,1,0),(0,1,1,2)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_1,j_2}a_{i_1,j_3}a_{i_2,j_2}a_{i_2,j_3}a_{i_2,j_4}^2\\\notag
&={{\rm Diag}}'({\mathbf}B) {\left(}{\mathbf}B{\circ}{\mathbf}B{\right)}{{\rm Diag}}({\mathbf}B),\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,2,2,2\}}^{(4,4)}[(1,1,1,1),(1,1,1,1)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}a_{i_1,j_2}a_{i_1,j_3}a_{i_1,j_4}a_{i_2,j_1}a_{i_2,j_2}a_{i_2,j_3}a_{i_2,j_4}\\\notag
&={{\rm tr}}{\left(}{\mathbf}B{\circ}{\mathbf}B{\right)}^2,\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,2,4\}}^{(4,4)}[(2,1,1),(0,1,3)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_1,j_2}a_{i_1,j_3}a_{i_2,j_2}a_{i_2,j_3}^3={{\rm tr}}{\mathbf}B {\mathbf}{D_B} {\mathbf}A {\mathbf}A'^{{\circ}3},\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,2,4\}}^{(4,4)}[(2,0,2),(0,2,2)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_1,j_3}^2a_{i_2,j_2}^2a_{i_2,j_3}^2\\\notag
&={{\rm Diag}}'({\mathbf}B) {\left(}{\mathbf}A{\circ}{\mathbf}A{\right)}{\left(}{\mathbf}A{\circ}{\mathbf}A{\right)}'{{\rm Diag}}({\mathbf}B) ,\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,2,4\}}^{(4,4)}[(1,1,2),(1,1,2)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}a_{i_1,j_2}a_{i_1,j_3}^2a_{i_2,j_1}a_{i_2,j_2}a_{i_2,j_3}^2\\\notag
&={{\rm tr}}{\left(}{\left(}{\mathbf}B{\circ}{\mathbf}B{\right)}{\left(}{\mathbf}A{\circ}{\mathbf}A{\right)}{\left(}{\mathbf}A{\circ}{\mathbf}A{\right)}'{\right)},\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{4,4\}}^{(4,4)}[(3,1),(1,3)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^3a_{i_1,j_2}a_{i_2,j_1}a_{i_2,j_2}^3={{\rm tr}}{\left(}{\left(}{\mathbf}A^{{\circ}3}{\mathbf}A'{\right)}{\left(}{\mathbf}A^{{\circ}3}{\mathbf}A'{\right)}'{\right)},\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{4,4\}}^{(4,4)}[(2,2),(2,2)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_1,j_2}^2a_{i_2,j_1}^2a_{i_2,j_2}^2={{\rm tr}}{\left(}{\left(}{\mathbf}A {\circ}{\mathbf}A{\right)}{\left(}{\mathbf}A {\circ}{\mathbf}A{\right)}'{\right)}^2 ,\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,6\}}^{(4,4)}[(1,3),(1,3)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}a_{i_1,j_2}^3a_{i_2,j_1}a_{i_2,j_2}^3={{\rm tr}}{\left(}{\mathbf}B {\mathbf}A^{{\circ}3} {\mathbf}A'^{{\circ}3}{\right)},\end{aligned}$$

$$\begin{aligned}
&\Omega_{\{2,6\}}^{(4,4)}[(2,2),(0,4)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_1,j_2}^2a_{i_2,j_2}^4={{\rm tr}}{\left(}{\left(}{\mathbf}A' {\mathbf}D_{{\mathbf}B}{\mathbf}A {\right)}{\circ}{\left(}{\mathbf}A'^{{\circ}2}{\mathbf}A^{{\circ}2}{\right)}{\right)},\end{aligned}$$

and

$$\begin{aligned}
&\Omega_{\{8\}}^{(4,4)}[(4),(4)]=\sum_{i_1,\cdots,i_2,j_1,\cdots, j_4}a_{i_1,j_1}^42a_{i_2,j_1}^4={\mathbf}1'{\mathbf}A^{{\circ}4}{\mathbf}A'^{{\circ}4}{\mathbf}1,\end{aligned}$$

Using the same procedure, we have

$$\begin{aligned}
&{\rm Var}({\mathbf}T_2)=n^{-2}{\left(}{{\rm E}}{\left(}\sum_{i}\widehat{\varepsilon_i}^2{\right)}^4-{{\rm E}}^2{\left(}\sum_{i}\widehat{\varepsilon_i}^2{\right)}^2{\right)}\\\notag
=&n^{-2}\sum_{i_1,\cdots,i_4,j_1,\cdots,j_8}a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_3}a_{i_2,j_4}a_{i_3,j_5}a_{i_3,j_6}a_{i_4,j_7}a_{i_4,j_8}
{\left(}{{\rm E}}\prod_{t=1}^8\xi_{j_t}-{{\rm E}}\prod_{t=1}^4\xi_{j_t}{{\rm E}}\prod_{t=5}^8\xi_{j_t}{\right)}\\\notag
=&n^{-2}(P_{2,1}+P_{2,2})+O(1),\end{aligned}$$

where

$$\begin{aligned}
P_{2,1}= {\mathbf}C_2^1{\mathbf}C_2^1{\mathbf}C_2^1\sum_{i_1,\cdots,i_4,j_1,\cdots, j_4}a^2_{i_1,j_1}a_{i_2,j_2}a_{i_3,j_2}a_{i_2,j_3}a_{i_3,j_3}a_{i_4,j_4}^2
=8{\left(}{{\rm tr}}{\mathbf}B{\right)}^2{{\rm tr}}{\mathbf}B^2,\end{aligned}$$

$$\begin{aligned}
P_{2,2}=\nu_4{\mathbf}C_2^1{\mathbf}C_2^1\sum_{i_1,\cdots,i_4,j_1,j_2, j_3}a^2_{i_1,j_1}a^2_{i_2,j_2}a^2_{i_3,j_2}a_{i_4,j_3}^2
=4\nu_4{{\rm tr}}({\mathbf}B'{\circ}{\mathbf}B'){\left(}{{\rm tr}}{\mathbf}B{\right)}^2,\end{aligned}$$

Similarly, we have

$$\begin{aligned}
&{\rm Cov}({\mathbf}T_1,{\mathbf}T_2)=n^{-1}{\left(}{{\rm E}}{\left(}\sum_{i}\widehat{\varepsilon_i}^2{\right)}^2\sum_{i}\widehat{\varepsilon_i}^4-{{\rm E}}{\left(}\sum_{i}\widehat{\varepsilon_i}^2{\right)}^2\sum_{i}\widehat{\varepsilon_i}^4{\right)}\\\notag
=&n^{-1}\sum_{i_1,\cdots,i_3,j_1,\cdots,j_8}a_{i_1,j_1}a_{i_1,j_2}a_{i_2,j_3}a_{i_2,j_4}a_{i_3,j_5}a_{i_3,j_6}a_{i_3,j_7}a_{i_3,j_8}
{\left(}{{\rm E}}\prod_{t=1}^8\xi_{j_t}-{{\rm E}}\prod_{t=1}^4\xi_{j_t}{{\rm E}}\prod_{t=5}^8\xi_{j_t}{\right)}\\
=&n^{-1}(P_{3.1}+P_{3,2}+P_{3,3}+P_{3,4})+O(1),\end{aligned}$$

where

$$\begin{aligned}
P_{3,1}={\mathbf}C_4^2{\mathbf}C_2^1{\mathbf}C_2^1\sum_{i_1,\cdots,i_3,j_1,\cdots, j_4}a_{i_1,j_1}^2a_{i_2,j_2}a_{i_2,j_3}a_{i_3,j_2}a_{i_3,j_3}a_{i_3,j_4}^2
=24{{\rm tr}}{\left(}{\mathbf}B^2{\circ}{\mathbf}B{\right)}{{\rm tr}}{\mathbf}B,\end{aligned}$$

$$\begin{aligned}
P_{3,2}=\nu_4{\mathbf}C_4^1{\mathbf}C_2^1{\mathbf}C_2^1\sum_{i_1,\cdots,i_3,j_1,\cdots, j_3}a_{i_1,j_1}^2a_{i_2,j_2}a_{i_3,j_2}a_{i_2,j_3}a_{i_3,j_3}^3
=16\nu_4{{\rm tr}}({\mathbf}B{\mathbf}A {\mathbf}A'^{{\circ}3}){{\rm tr}}{\mathbf}B,\end{aligned}$$

$$\begin{aligned}
P_{3,3}=\nu_4{\mathbf}C_4^2{\mathbf}C_2^1\sum_{i_1,\cdots,i_3,j_1,\cdots, j_3}a_{i_1,j_1}^2a^2_{i_2,j_2}a^2_{i_3,j_2}a_{i_3,j_3}^2
=12\nu_4{{\rm tr}}{\left(}{\left(}{\mathbf}A'{\mathbf}D_{{\mathbf}B}{\mathbf}A{\right)}{\circ}{\left(}{\mathbf}A'{\mathbf}A{\right)}{\right)}{{\rm tr}}{\mathbf}B,\end{aligned}$$

$$\begin{aligned}
\label{ct12}
P_{3,3}=\nu_6{\mathbf}C_2^1\sum_{i_1,i_2,i_3,j_1, j_2}a_{i_1,j_1}^2a^2_{i_2,j_2}a^4_{i_3,j_2}
=2\nu_6[{{\rm Diag}}({\mathbf}A'{\mathbf}A)'({\mathbf}A'^{{\circ}4}){\mathbf}1]{{\rm tr}}{\mathbf}B,\end{aligned}$$

We would like to point out that we do not need the assumption that $H_0$ holds up to now. From now on, in order to simplify the above formulas we assume $H_0$ holds. Summarizing the calculations above, we obtain under $H_0$

$$\begin{aligned}
{{\rm E}}{\mathbf}T_1=3\sum_ib_{i,i}^2+\nu_4\sum_{ij}a_{ij}^4=3{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^2,\end{aligned}$$

$$\begin{aligned}
\label{e2}
{{\rm E}}{\mathbf}T_2=n^{-1}{\left(}{\left(}n-p{\right)}^2+2{\left(}n-p{\right)}+\nu_4{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P){\right)},\end{aligned}$$

$$\begin{aligned}
{\rm Var}{\mathbf}T_1=&72{{\rm Diag}}'({\mathbf}P) {\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}{{\rm Diag}}({\mathbf}P)+24{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}^2\\\notag
&+\nu_4{\left(}96{{\rm tr}}{\mathbf}P {\mathbf}{D_P} {\mathbf}P {\mathbf}P^{{\circ}3}+72{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^3+36{{\rm Diag}}'({\mathbf}P) {\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}^2{{\rm Diag}}({\mathbf}P) {\right)}\\\notag
&+\nu^2_4{\left(}18{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)^4+16{{\rm tr}}({\mathbf}P^{{\circ}3}{\mathbf}P)^2){\right)}\\\notag
&+\nu_6{\left(}12{{\rm tr}}{\left(}{\left(}{\mathbf}P{\mathbf}D_{{\mathbf}P}{\mathbf}P {\right)}{\circ}{\left(}{\mathbf}P^{{\circ}2}{\mathbf}P^{{\circ}2}{\right)}{\right)}+16{{\rm tr}}{\mathbf}P {\mathbf}P^{{\circ}3}{\mathbf}P^{{\circ}3}{\right)}+\nu_8{\mathbf}1'({\mathbf}P^{{\circ}4}{\mathbf}P^{{\circ}4}){\mathbf}1,\end{aligned}$$

$$\begin{aligned}
{\rm Var}({\mathbf}T_2)=\frac{8{\left(}n-p{\right)}^3+4\nu_4{\left(}n-p{\right)}^2{{\rm tr}}({\mathbf}P{\circ}{\mathbf}P)}{n^2}+O(1),\end{aligned}$$

and $$\begin{aligned}
&{\rm Cov}({\mathbf}T_1,{\mathbf}T_2)\\\notag
=&\frac{{\left(}n-p{\right)}}{n}{\left(}24{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}+16\nu_4{{\rm tr}}({\mathbf}P {\mathbf}P^{{\circ}3})+12\nu_4{{\rm tr}}{\left(}{\left(}{\mathbf}P{\mathbf}D_{{\mathbf}p}{\mathbf}P{\right)}{\circ}{\mathbf}P{\right)}+2\nu_6[{{\rm Diag}}({\mathbf}P)'({\mathbf}P^{{\circ}4}){\mathbf}1]{\right)}.\end{aligned}$$

The proof of the main theorem
-----------------------------

Define a function $f(x,y)=\frac{x}{y}-1$. One may verify that $f_x(x,y)=\frac{1}{y},$ $f_y(x,y)=-\frac{x}{y^2}$, where $f_x(x,y)$ and $f_y(x,y)$ are the first order partial derivative. Since ${\mathbf}T=\frac{{\mathbf}T_1}{{\mathbf}T_2}-1,$ using the delta method, we have under $H_0$, $${{\rm E}}{\mathbf}T=f({{\rm E}}{\mathbf}T_1,{{\rm E}}{\mathbf}T_2)={\left(}\frac{3n{{\rm tr}}{\left(}{\mathbf}P{\circ}{\mathbf}P{\right)}}{(n-p)^2+2{\left(}n-p{\right)}}-1{\right)},$$ $$\begin{aligned}
  {\rm {Var}} {\mathbf}T=(f_x({{\rm E}}{\mathbf}T_1,{{\rm E}}{\mathbf}T_2),f_y({{\rm E}}{\mathbf}T_1,{{\rm E}}{\mathbf}T_2))\Sigma(f_x({{\rm E}}{\mathbf}T_1,{{\rm E}}{\mathbf}T_2),f_y({{\rm E}}{\mathbf}T_1,{{\rm E}}{\mathbf}T_2))'.\end{aligned}$$ The proof of the main theorem is complete.

[10]{}

Adelchi Azzalini and Adrian Bowman. On the use of nonparametric regression for checking linear relationships. , pages 549–557, 1993.

Zhidong Bai and Jack W Silverstein. , volume 20. Springer, 2010.

Trevor S Breusch and Adrian R Pagan. A simple test for heteroscedasticity and random coefficient variation. , pages 1287–1294, 1979.

R Dennis Cook and Sanford Weisberg. Diagnostics for heteroscedasticity in regression. , 70(1):1–10, 1983.

Holger Dette and Axel Munk. Testing heteroscedasticity in nonparametric regression. , 60(4):693–708, 1998.

Holger Dette, Axel Munk, and Thorsten Wagner. Estimating the variance in nonparametric regression—what is a reasonable choice? , 60(4):751–764, 1998.

Herbert Glejser. A new test for heteroskedasticity. , 64(325):316–323, 1969.

Michael J Harrison and Brendan PM McCabe. A test for heteroscedasticity based on ordinary least squares residuals. , 74(366a):494–499, 1979.

S John. Some optimal multivariate tests. , 58(1):123–127, 1971.

Zhaoyuan Li and Jianfeng Yao. Homoscedasticity tests valid in both low and high-dimensional regressions. , 2015.

Gary C McDonald and Richard C Schwing. Instabilities of regression estimates relating air pollution to mortality. , 15(3):463–481, 1973.

Halbert White. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. , pages 817–838, 1980.

H.Altay Guvenir and I.Uysal. Bilkent University Function Approximation Repository. , 2000

[^1]: Zhidong Bai is partially supported by a grant NSF China 11571067

[^2]: G. M. Pan was partially supported by a MOE Tier 2 grant 2014-T2-2-060 and by a MOE Tier 1 Grant RG25/14 at the Nanyang Technological University, Singapore.

[^3]: Yanqing Yin was partially supported by a project of China Scholarship Council
