1 Introduction
###### Abstract

Testing heteroscedasticity of the errors is a major challenge in high-dimensional regressions where the number of covariates is large compared to the sample size. Traditional procedures such as the White and the Breusch-Pagan tests typically suffer from low sizes and powers. This paper proposes two new test procedures based on standard OLS residuals. Using the theory of random Haar orthogonal matrices, the asymptotic normality of both test statistics is obtained under the null when the degrees of freedom tend to infinity. This encompasses both the classical low-dimensional setting where the number of variables is fixed while the sample size tends to infinity, and the proportional high-dimensional setting where these dimensions grow to infinity proportionally. These procedures thus offer a wide coverage of dimensions in applications. To our best knowledge, this is the first procedures in the literature for testing heteroscedasticity which are valid for medium and high-dimensional regressions. The superiority of our proposed tests over the existing methods are demonstrated by extensive simulations and by several real data analyses as well.

Keywords. Breusch and Pagan test, White’s test, heteroscedasticity, high-dimensional regression, hypothesis testing, Haar matrix.

Testing for Heteroscedasticity in High-dimensional Regressions

[2mm] Zhaoyuan Li and Jianfeng Yao

Department of Statistics and Actuarial Science

The University of Hong Kong

## 1 Introduction

Consider the linear regression model

 yi=Xiβ+εi,i=1,…,n, (1)

where is the dependent variable, is a vector of regressors, is the -dimensional coefficient vector, and the error , where could depend on covariates , and are independent standard normal distributed. A significant part of the inference theory for the model is based on the assumption that the errors are homoscedastic, i.e. under the hypothesis

 H0:σ21=…=σ2n=σ2, (2)

for some constant , that is, the unconditional and conditional variances of the noise coincide and are independent of the covariates. However, this assumption cannot be always guaranteed in practice, and it is well known that heteroscedasticity of the error variance leads to inefficient parameter estimates and inconsistent covariance estimates. We consider testing the hypothesis in (2) when the number of covariates goes to infinity together with the sample size .

Studying this testing problem is also motivated by recent advances in the estimation of high-dimensional regressions. In this paper, we consider testing the hypothesis in (2) when the number of covariates can be large with respect to the sample size . High-dimensional regressions become vital due to the increasingly wide availability of data sets with a large number of variables in empirical economics, finance (Belloni et al., 2014a) and biology (Daye et al., 2012). For example, the American Housing Survey records prices as well as a multitude of features of the house sold; scanner datasets record prices and numerous characteristics of products sold at a store or on the internet (Belloni et al., 2011); and text data frequently lead to counts of words in documents from a large dictionary (Taddy, 2013). Importantly, heteroscedasticity is possible in such data sets. Among the existing methods for high-dimensional regressions, one approach assumes sparsity where the number of important regressors is much smaller than or . For example, Belloni et al. (2012) and Belloni et al. (2014b) studied the estimation problem with heteroscedasticity and proposed the heteroscedastic form of Lasso and square-root Lasso methods, respectively. However, if the errors are homoscedastic, these heteroscedasticity-consistent methods may lose efficiency as suggested by the phenomenon arising in low-dimensional regressions. Here, we conduct a small simulation study with 5000 replications to illustrate this point by using data generated according to Model 1 in Section 3 with . The ratio of the heteroscedastic form of Lasso estimator (root mean squared error) over the OLS estimator (root mean squared error) is 0.8791 when the errors are related to 10 regressors. But the ratio is 12.08 for homoscedastic errors, and 7.53 when the errors are related to only one regressor. El Karoui et al. (2013) and Bean et al. (2013) stated that the Lasso-type of methods result in biased estimates of the coefficients, and the least squares method is preferable to other M-estimators in high-dimensional regression under homoscedasticity. Bean et al. (2013) proposed an optimal least square algorithm with the assumption that the error is homoscedastic with a known distribution. However, the performance of this optimal algorithm is largely unknown if the error is in fact heteroscedastic. In summary, the discussion above on two recent high-dimensional estimation methodologies highlights the importance of conducting heteroscedasticity detection as a preliminary step in practice in order to select a suitable estimation method for high-dimensional regressions.

Heteroscedastic testing has been extensively studied for classical low-dimensional regressions in the literature. Many popular tests examine whether the estimated residuals are correlated with some covariates or any auxiliary variables that would be useful in explaining the departure from homoscedasticity, see for example Breusch and Pagan (1979), White (1980), Cook and Weisberg (1983), Azzalini and Bowman (1993), Diblasi and Bowman (1997), and Su and Ullah (2013). These tests, however, will not have much power if the existing heteroscedasticity is not strongly related to either the chosen auxiliary variables or covariates. In consequence, many nonparametric test procedures are thus proposed to avoid such potential model misspecification, see for example, Eubank and Thomas (1993) and Dette and Munk (1998). Koenker and Bassett (1982) and Newey and Powell (1987) proposed to test heteroscedasticity by comparing different quantile or expectile estimates. Their approach is much preferable to many other tests for heavy tailed errors (Lee, 1992). However, there is some difficulty in applying this approach because no clear criterion exists for selecting the used quantiles.

Testing the homoscedasticity hypothesis (2) becomes very challenging for high-dimensional regressions. The large sample theory of all the existing tests discussed above is developed under the low-dimensional framework where the dimension should be fixed while the sample size tends to infinity. By referring to recent advances in high-dimensional statistics (Paul and Aue, 2014; Yao et al., 2015), it clearly appears that these test methods are not suitable for analysing data sets where the number of variables is not small enough” compared to the sample size. For example, the limiting approximation for White’s test statistic is typically misleading even for a moderate dimension while the sample size is (see Table 1 for more details). As an additional illustration, many published Monte Carlo studies of tests for heteroscedasticity have used very low-dimensional designs and the error variances are determined by a single variable in the alternative model, see for example Dette and Munk (1998). Godfrey and Orme (1999) and Godfrey (1996) showed that the results obtained from very simple experimental designs (for example ) may be an unreliable guide to finite sample performance with a moderately large number of variables. Another illustration of high-dimensional effect is an interesting phenomenon shown in Ferrari and Cribari-Neto (2002) and Godfrey and Orme (1999) where the actual size of many popular tests stays far from the nominal level for the moderately large sample size . Therefore, accurate and powerful test procedure is an urgent need for detecting heteroscedasticity in a high-dimensional regression.

In this paper, we propose two new procedures for testing heteroscedasticity, which are dimension-proof in the sense that they are valid for a wide range of dimension (covering both low and high-dimensional settings). More precisely, our procedures are theoretically valid once the degree of freedom is large enough (precisely when ). This includes for instance the low-dimensional setting where and the high-dimensional situation where and grow to infinity proportionally such that with . Simulation experiments reported in Section 3 show that the proposed tests outperform the popular existing methods for medium or high-dimensional regressions. More surprisingly, even in low-dimensional setting, our procedures perform better than these classical procedures.

The paper is organised as follows. The main results of the paper are reported in Section 2. Two new tests are here proposed using the residuals of a least squares fit. Section 3 reports several simulation experiments to assess the finite sample performance of the proposed tests and compare them to the existing ones. In Section 4 we apply the suggested procedures to analyse four real data sets. All technical proofs of the results presented in Section 2 are relegated to the Appendix.

## 2 Main results

The following assumptions will be used in our set-up of the regression model (1):

• Assumption (a): The errors are independent and normal distributed: ;

• Assumption (b): In the design matrix , are independent normal distributed vectors with mean 0 and covariance matrix ;

• Assumption (c): As , the degree of freedom ;

• Assumption (d): In addition to Assumption (c), , where .

Both Assumptions (a) and (b) are classical in a regression model. Assumptions (c) and (d) define the asymptotic setting of the paper which is quite general. In particular, the setting includes the situation where both and are large while remaining comparable, i.e. for some , and . Meanwhile, the setting encompasses the classical low-dimensional situation where is a constant and . Therefore, the procedure derived under this setting will be applicable to both the high and low-dimensional settings. It is however noted that since our methods will use the OLS residuals, it is required that although both dimensions can grow to infinity.

In the regression model (1) and under homoscedasticity, the parameter vector is estimated by the OLS estimator where and . Then, the vector of residuals is

 (3)

Here and throughout of the paper, denotes the -th order identity matrix. Notice that is a projection matrix of rank . In the following, two test statistics are proposed based on the residuals .

Note that each covariate vector so we have where . Let be the corresponding design” matrix. Then we have . Therefore the projection matrix is independent of the covariance structure . In what follows we can assume and the coordinates of are i.i.d standard normals.

### 2.1 An approximate likelihood-ratio test

We first derive a test statistic from the concept of likelihood ratio test. For the regression model (1) and under Assumption (a), the likelihood function is simply

 L(β,σ21,…,σ2n)=(2π)−n/2(σ21⋯σ2n)−1/2exp{−12n∑i=1(yi−Xiβ)2σ2i}.

Without assuming the homoscedasticity, the likelihood is maximised by solving the system of equations

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∂logL∂σ2i=−12σ2i+12σ4i(yi−Xiβ)2=0,∂logL∂β=−12n∑i=12(yi−Xiβ)σ2i(−Xi)=0,1≤i≤n.

Therefore, the maximum likelihood estimator (MLE) of satisfy the system of equations

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩^σ2i=(yi−Xi^β)2,^β=(n∑i=1X′iXi^σ2i)−1n∑i=1yiX′i^σ2i,1≤i≤n.

The corresponding maximized likelihood is

 L1=(2π)−n/2∏{(yi−Xi^β)2}−1/2exp(−n/2).

Notice that since the number of unknown parameters exceeds the sample size, this MLE cannot be a reliable estimator. Nevertheless, this likelihood concept will help us to define a meaningful test statistic for testing the homoscedasticity hypothesis as follows: we approximate the MLE in the maximized likelihood by the OLS to get an approximate value

 L∗1=(2π)−n/2∏{(yi−Xi^β0)2}−1/2exp(−n/2).

On the other hand under the homoscedasticity hypothesis, the OLS estimator and the estimator of the variance

 ^σ20=1nn∑i=1(yi−Xi^β0)2,

are in fact the MLEs. So the maximized likelihood under the null hypothesis is

 L0=(2π)−n/2(^σ20)−n/2exp(−n/2). (4)

Therefore, the approximate likelihood ratio, likelihood ratio is first derived by Mauchly (1940), is defined as

 L0L∗1=(^σ20)−n/2(∏ni=1(yi−Xi^β0)2)−1/2=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1n∑ni=1^ε2i(∏ni=1^ε2i)1/n⎫⎪ ⎪ ⎪⎬⎪ ⎪ ⎪⎭−n2,

where it is reminded that . This suggests to consider the approximate likelihood-ratio statistic

 T1=−2nlogL0L∗1=log1n∑ni=1^ε2i(∏ni=1^ε2i)1/n. (5)

Interestingly enough, the statistic depends on the ratio of the arithmetic mean of the squared residuals over their geometric mean: always and a large value of will indicate a significant deviation of the residuals from a constant, that is presence of heteroscedasticity. Meanwhile, this statistic has a scale-free property and is not affected by the magnitude of the variance under the null hypothesis. Therefore, without loss of generality for the study of , we assume that under the null. The asymptotic distribution of under the null is derived in the following theorem.

###### Theorem 1.

Assume that Assumptions (a)-(b)-(d) are satisfied for the regression model (1). Then under the null hypothesis of homoscedasticity, we have as

 √n(T1−[log2+γ])D⟶N(0,π22−2), (6)

where is the Euler constant.

The testing procedure using with the critical value from (6) is referred as the approximate likelihood-ratio test (ALRT). In addition to the scale-free property mentioned above, an attractive feature appears here is that the asymptotic distribution of is completely independent of , the relative magnitude of the dimension over the sample size . This prefigures a large applicability of the procedure to a wide range of combinations of in finite-sample situations. This robustness is indeed confirmed by the simulation study reported in Section 3.

The proof of Theorem 1 is based on the following lemma, which establishes the asymptotic limit of the joint distribution of and under the null.

###### Lemma 1.

Let be the sequence of the OLS residuals given in (3). Then, under and Assumptions (a)-(b)-(d), and as , we have

 Σ−1/21{(∑ni=1^ε2i∑ni=1log^ε2i)−μ1}D⟶N(0,I2), (7)

where

 μ1=(kn(−γ−log2+logcn)),

and

 Σ1=(2k2n2nn(π2/2+2/cn−2)).

The proofs of Lemma 1 and Theorem 1 are postponed to the appendix.

### 2.2 The coefficient-of-variation test

The departure of a sequence of numbers from a constant can also be efficiently assessed by its coefficient of variation. In multivariate analysis, this idea is closely related to optimal invariant tests, see John (1971). Applying this idea to the sequence of residuals leads to the following coefficient-of-variation statistic

 T2=1n∑ni=1(^ε2i−¯m)2¯m2,with ¯m=1nn∑i=1^ε2i. (8)

Obviously, the statistic becomes small and close to 0 under the null hypothesis of homoscedasticity, and larger under the alternative hypothesis of heteroscedasticity. Like the previous statistic , this statistic is also scale-free and again we can assume for under the null without loss of generality. The asymptotic distribution of under the null hypothesis is derived in the following theorem.

###### Theorem 2.

Assume that Assumptions (a)-(b)-(c) are satisfied for the regression model (1). Then under the null hypothesis of homoscedasticity, we have as

 √n(T2−2)D⟶N(0,24). (9)

The testing procedure using with the critical value from (9) is referred as the coefficient-of-variation test (CVT). Similar to the statistic , the asymptotic distribution of is also scale free and independent of , the relative magnitude of the dimension over the sample size .

The proof of Theorem 2 is based on the following lemma, which establishes the asymptotic limit of the joint distribution of and under the null.

###### Lemma 2.

Let be the sequence of the OLS residuals given in (3). Then, under and Assumptions (a)-(b)-(c), and as , we have

 Σ−1/22{(∑ni=1^ε4i∑ni=1^ε2i)−μ2}D⟶N(0,I2), (10)

where

 μ2=(3k(k+2)n+2k),

and

 Σ2=⎛⎜⎝24k4n3+72k3n212k2n+212k2n+22k⎞⎟⎠.

The proofs of Lemma 2 and Theorem 2 are postponed in the appendix.

## 3 Simulation experiments

We have undertaken an extensive simulation study to investigate the finite sample performance of the proposed tests, ALRT and CVT. Comparisons are also made with several existing popular methods: the BP test, proposed by Breusch and Pagan (1979) and modified by Koenker (1981); the White test (White, 1980); and the DM test (Dette and Munk, 1998).

Breusch and Pagan (1979) constructed a general test statistic, assuming that the conditional variance has a known functional form , where and . They proposed a Lagrange multiplier statistic to test the joint null hypothesis of while the intercept is unspecified. Koenker (1981) modified this test in order to improve its empirical size. This test has been widely used in the literature and is the representative one in the family of Lagrange multiplier or score tests, as it includes many other tests (e.g. Cook and Weisberg, 1983 and Eubank and Thomas, 1993) as special cases.

The White test fits an artificial regression of the squared OLS residuals on the elements of the lower triangle of the matrix , and the test statistic is the squared multiple correlation coefficient from this regression. The author proved that the statistic is asymptotically distributed as with degrees of freedom under the null hypothesis of homoscedasticity (as the sample size tends to infinity).

Dette and Munk (1998) proposed a nonparametric method, the DM test. It is constructed on estimation of empirical variance of expected squared residuals, and its asymptotic normality is given. This nonparametric test avoids the estimation of the regression curve directly, which makes it more robust and better than those tests based the estimated residuals.

### 3.1 Empirical sizes of the tests

We explore the performance of these tests using different combination of and . The sample sizes and ratios are considered. Each simulation is repeated 5000 times to test the stability of the method. Empirical size of a test is the percentage of rejected tested cases. According to the model (1), the design matrix are assumed to be multi-normal. The error is drawn from standard normal as the size and power of the proposed tests are invariant with respect to different scalings of variance function. The nominal test level is 5%.

Table 1 presents the empirical sizes of the ALRT, CVT, White and BP tests (values close to 5% are better). The proposed ALRT and CVT tests are consistently accurate in all tested combinations of (including the smallest ones); they largely outperform the White and BP tests. This good performance can be explained by a fast convergence in the limiting results of ALRT (Theorem 1) and CVT (Theorem 2). The ALRT test performs a little better than the CVT test for small value of the ratio , but the CVT test is preferred when is getting close to 1. The BP test loses its size from (approximately) 4% to 1% when the ratio increases from 0.05 to 0.5, while the White test has an empirical size of 0.16% when the ratio is and sample size is (Notice that this test is not applicable when due to its dimension-sample-size requirement ).

### 3.2 Empirical powers of the tests

To investigate the power of these tests, we follow Dette and Munk (1998) and consider the following three models with different error forms:

• Model 1: ;

• Model 2: ;

• Model 3: ;

where the vector is filled with elements and/or . The value corresponds to homoscedasticity, and we consider two levels of heteroscedasticity: with (1st component only) and (first 10% of components). Same setting with Section 3.1 is used and empirical powers of the tests are obtained using 5000 replications for each scenario.

Tables 2-4 present the empirical powers of the ALRT, CVT and BP tests for these three error models, respectively. Plots are also provided for the case of sample size for a easier comparison. The results of the White test are omitted here due to its worst performance in term of size in Table 1. As expected, for each model, the power becomes larger as the level of heteroscedasticity increases. In general, the empirical powers of all tests become smaller as the dimension goes up (ratio increases); the reason is that the BP test is not suitable for high-dimensional setting, and the ALRT and CVT tests are related to the degree of freedom of which becomes small when dimension increases. The CVT test is most powerful in all tested cases.

As for the three models considered, the results for Model 1 and Model 3 are similar with each other where the BP test show no power when while the ALRT and CVT tests have a reasonable power unless is close to 1. Recall that in such situation, the matrix is close to singularity, the OLS estimator is performing badly. However, our procedures still show a reasonable performance. The situation in Model 2 is radically different where the BP test has no power for all tested combinations of while the ALRT and CVT keep a reasonable power (unless is close to 1) as in Model 1 and 3. In conclusion, generally in all the tested situations, the proposed tests ALRT and CVT outperform the BP tests in a large extent. Table 2: Empirical powers of the ALRT, CVT and BP tests for Model 1 under two scenarios with sample size n=100,500,1000 and varying ratio p/n. Table 3: Empirical powers of the ALRT, CVT and BP tests for Model 2 under two scenarios with sample size n=100,500,1000 and varying ratio p/n. Table 4: Empirical powers of the ALRT, CVT and BP tests for Model 3 under two scenarios with sample size n=100,500,1000 and varying ratio p/n.

### 3.3 Non-Gaussian design

In this section, we investigate the performance of all tests applied to non-Gaussian design matrix. Here the entries of the design matrix are drawn from gamma distribution and uniform distribution , respectively. Except the design matrix , the same setting with Section 3.1 is used to obtain the empirical sizes of all tests, and the same setting with Section 3.2 is used to obtain the empirical powers of all tests. All results are obtained using 5000 replications for each scenario.

The empirical sizes and powers are presented in Tables 5 and 6, repectively. We find that there is no significant difference in terms of size and power between these two non-normal designs and the previously reported normal design. Similarly, the proposed ALRT and CVT perform well in all models and they are much better than the BP test. This suggests that the proposed tests are robust against the form or the distribution of the design matrix.

Simulation study is also conducted to explore the performance of these tests for fixed design. The design matrix is generated once and keep same for all replications. Even though our theoretic results are developed in the random design only, the inclusion of the fixed design simulation study is motivated by the believe that these asymptotic results of the ALRT and CVT tests remain useful in fixed design. As expected, the simulation results of empirical sizes and powers in fixed design are all similar to that in random design. These results are omitted here for brevity.

### 3.4 Small sample sizes

Simulation experiments are conducted to assess the performance of our tests for small sample size in a classical low-dimensional scenario. The DM test is compared here (notice that this test is not in Tables 1-4 since its implementation in a multivariate setting is unclear). Following the same set-up of Dette and Munk (1998), the design points are chosen as and the sample sizes are . The BP and White tests are not considered in this part due to the fact that the design matrix in the setting considered here is nearly singular, so that the OLS estimates used by these two tests are unreliable. The considered model is with three settings:

• S1: ,

• S2: ,

• S3: ,

with different values for (0, 0.5 and 1.0). is the mean function, so the linear model tested here is one dimension. And is the error term. The case corresponds to the null hypothesis of homoscedasticity and the choices and 1 correspond to two alternatives. We calculated the proportion of rejections of the tests using 5000 simulations for each scenario.

The empirical sizes and powers of the ALRT, CVT and DM tests are summarised in Table 7. The results of the DM test are from Tables 1 and 2 of Dette and Munk (1998). In term of empirical size, the ALRT test is conservative while the DM test is inclined to overestimate the size and both of them are close to the nominal level 0.05. But the ALRT test is more powerful than the DM test for settings S2 and S3. The ALRT test has similar performance with the DM test in setting S1 because it runs the OLS estimation for the sinusoidal mean function. The CVT only performs better than the DM test in term of power in several cases. Therefore, although the ALRT test is constructed under the high-dimensional framework, it is still a competitive procedure in classical low-dimensional regression even with a small sample size. This is also supported by the results for the cases in Tables 1-4.

## 4 Real data analyses

Though the newly proposed two tests seem to perform better than the classical ones in the simulation experiments, we now compare them on several real examples. According to the results of simulation, we use the BP test as the representation of classical tests.

### 4.1 Low-dimensional data sets

In order to check the performance of the proposed tests in low-dimensional situation, we analyse two data sets: the ‘bond yield’ data and the ‘currency substitution’ data111These two data sets are available in the R package ‘lmtest’.. The bond yield data set is a multivariate quarterly time series from 1961(1) to 1975(4) (sample size ) with seven variables, including RAARUS (difference of interest rate on government and corporate bonds), MOOD (measure of consumer sentiment), EPI (index of employment pressure), EXP (interest rate expectations), Y (joint proxies for the impact of callability) and K (artificial time series based on RAARUS). This data set is used to analyse the observed long-term bond yield differentials for different types of instruments. Two main works are Cook and Hendershott (1978) in which a linear regression of RAARUS on MOOD, EPI, EXP and RUS is fitted to find the factors contributed to the bond yield spreads, and Yawitz and Marshall (1981) in which another linear regression of RAARUS on MOOD, Y and K is fitted to see the effect of callability on bond yields. To investigate whether the homoscedasticity assumption in both models is justified, we applied the BP test, the ALRT test and the CVT test to each regression model. For the Cook-Hendershott model, we got three p-values of 0.5614 (BP), 0.3307 (ALRT) and 0.8333 (CVT). And the Yawitz-Marshall model yields three p-values of 0.3838 (BP), 0.7314 (ALRT) and 0.3885 (CVT). Hence, these tests show no evidence against the assumption of constant variability in both models.

The currency substitution data set is a multivariate quarterly times series from 1960(4) to 1975(4) (sample size ) with four variables, including logCUS (logarithm of the ratio of Canadian holdings of Canadian dollar balances and Canadian holdings of U.S. dollar balances), Iu (yield on U.S. treasury bills), Ic (yield on Canadian treasury bills) and logY (logarithm of Canadian real gross national product). This data set is used to analyse the effect of flexible exchange rates and studied by Bordo and Choudhri (1982) where a linear model is fitted for logCUS using the other three variables as covariates. Their results were obtained under the assumption that the error variances are constant, which is supported by our proposed test: the ALRT test reports a p-value of 0.5779 and the CVT test reports a p-value of 0.1309. However, the p-value obtained by the BP test is 0.01324 which is inconsistent with the results in Bordo and Choudhri (1982).

### 4.2 High dimensional data sets

In this part, we evaluate the performance of our proposed tests on two data sets with medium and high dimensions: the ‘international economic growth’ data222Available on the web-site: https://stuff.mit.edu/ vchern/NBER/ Belloni et al. (2011) and the ‘eminent-domain’ data333Available on the web-site: http://faculty.chicagobooth.edu/christian.hansen/research/#Code. (Belloni et al., 2012).

The international economic growth data set concerns the national growth rates in GDP per capita with covariates including education, science policies, strength of market institutions, trade openness, saving rates and others. The sample size is . There is no unmeasurable underlying variable in this example so the regression model with all variables has constant disturbance. The CV and BP tests provide same conclusion by reporting p-values 0.5822 and 0.9436, respectively. Belloni et al. (2011) used covariate selection procedure to select significant variables among 62 covariates and the variable black market premium” is selected. Actually, this variable has important economic meaning as it characterises trade openness. Hence, the regression model without this variable will have heteroscedastic errors, and this conjecture is supported by our proposed CV test with a p-value of 0.0686 (compared with the value 0.5822 in full model). However, the p-value obtained by the BP test is 0.9808 which is inconsistent with the result in Belloni et al. (2011).

Belloni et al. (2012) studied the effects of federal appellate court decisions regarding eminent domain on a variety of economic outcomes. To explore the effect of the characteristics of three-judge panels on judicial decisions, the data set ‘eminent-domain’ containing explanatory variables (gender, race, religion, political affiliation, etc.) is used with sample size . The ratio of dimension and sample size is larger than 0.8. Belloni et al. (2012) argued that much heteroscedasticity exists in this data set and used heteroscedasticity consistent standard error estimator in their analysis. Applying ALRT and CVT tests on this data set, we found a p-value of and , respectively, strongly supporting these authors’ approval. On the other hand, the BP test cannot detect the existence of heteroscedasticity by reporting a p-value of 0.3331.

These results of real data sets analysis are consistent with the conclusion drawn from the simulation part that our newly proposed tests can provide accurate detection of heteroscedasticity under the medium or high dimensional situations, while the BP test, constructed under the low-dimensional scheme, not only cannot possess a correct size, but also loses power when heteroscedasticity exists.

## 5 Conclusion and discussion

For high-dimensional linear regression model, we propose two simple and efficient tests to detect the existence of heteroscedasticity. The asymptotic normalities of test statistics with simple form are constructed under the assumption that the degree of freedom is large compared to the sample size with as and are thus appropriate for analyzing high-dimensional data sets. Extensive Monte-Carlo experiments demonstrates the superiority of our proposed tests over some popular existing methods in terms of size and power. The good performance of our tests is also illustrated by several real data analyses. Surprisingly enough, these high-dimensional tests when used in the tested low-dimensional situations also show a performance comparable to that of the existing classical tests which are designed specifically under low-dimensional scheme.

There are still several avenues for future research. For example, the asymptotic results of the tests proposed here are based on the normality assumption for both the error and the random design. It is highly valuable to investigate the non-Gaussian setting. Although we have shown some robustness of the proposed procedures against non-Gaussian design in simulation experiments, a thorough investigation is missing. It is however clear that new theoretical tools will be needed to tackle with such non-Gaussian setting.

Lastly, our procedures rely on the OLS residuals, therefore have some limitations. First, it is required that even though both of them can be large. How to address the case where remains an open question. Second, it is well-known that the OLS estimates lack robustness against outliers. It is very likely that our tests possess same weakness.

## References

• Azzalini and Bowman (1993) A. Azzalini and A. Bowman. On the use of nonparametric regression for checking linear relationships. Journal of the Royal Statistical Society. Series B (Methodological), 55(2):549–557, 1993.
• Bean et al. (2013) D. Bean, P. J. Bickel, N. El Karoui, and B. Yu. Optimal m-estimation in high-dimensional regression. Proceedings of the National Academy of Sciences, 110(36):14563–14568, 2013.
• Belloni et al. (2011) A. Belloni, V. Chernozhukov, and C. Hansen. Inference for high-dimensional sparse econometric models. Advances in Economics and Econometrics, 10th World Congress of Econometric Society, 2011.
• Belloni et al. (2012) A. Belloni, D. Chen, V. Chernozhukov, and C. Hansen. Sparse models and methods for optimal instruments with an application to eminent domain. Econometrica, 80(6):2369–2429, 2012.
• Belloni et al. (2014a) A. Belloni, V. Chernozhukov, and C. Hansen. High-dimensional methods and inference on structural and treatment effects. The Journal of Economic Perspectives, 28(2):29–50, 2014a.
• Belloni et al. (2014b) A. Belloni, V. Chernozhukov, and L. Wang. Pivotal estimation via square-root lasso in nonparametric regression. The Annals of Statistics, 42(2):757–788, 2014b.
• Bordo and Choudhri (1982) M. D. Bordo and E. U. Choudhri. Currency substitution and the demand for money: some evidence for canada. Journal of Money, Credit and Banking, 14(1):48–57, 1982.
• Breusch and Pagan (1979) T. S. Breusch and A. R. Pagan. A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5):1287–1294, 1979.
• Cook and Weisberg (1983) R. D. Cook and S. Weisberg. Diagnostics for heteroscedasticity in regression. Biometrika, 70(1):1–10, 1983.
• Cook and Hendershott (1978) T. Q. Cook and P. H. Hendershott. The impact of taxes, risk and relative security supplies on interest rate differentials. The Journal of Finance, 33(4):1173–1186, 1978.
• Daye et al. (2012) Z. J. Daye, J. Chen, and H. Li. High-dimensional heteroscedastic regression with an application to eqtl data analysis. Biometrics, 68(1):316–326, 2012.
• Dette and Munk (1998) H. Dette and A. Munk. Testing heteroscedasticity in nonparametric regression. Journal of the Royal Statistical Society: Series B (Methodology), 60(4):693–708, 1998.
• Diblasi and Bowman (1997) A. Diblasi and A. Bowman. Testing for constant variance in a linear model. Statistics & Probability letters, 33(1):95–103, 1997.
• El Karoui et al. (2013) N. El Karoui, D. Bean, P. J. Bickel, C. H. Lim, and B. Yu. On robust regression with high-dimensional predictors. Proceedings of the National Academy of Sciences, 110(36):14557–14562, 2013.
• Eubank and Thomas (1993) R. L. Eubank and Will Thomas. Detecting heteroscedasticity in nonparametric regression. Journal of the Royal Statistical Society: Series B (Methodological), 55(1):145–155, 1993.
• Ferrari and Cribari-Neto (2002) S. L. P. Ferrari and F. Cribari-Neto. Corrected modified profile likelihood heteroskedasticity tests. Statistics & Probability letters, 57(4):353–361, 2002.
• Godfrey (1996) L. G. Godfrey. Some results on the glejser and koenker tests for heteroskedasticity. Journal of Econometrics, 72(1):275–299, 1996.
• Godfrey and Orme (1999) L. G. Godfrey and C. D. Orme. The robustness, reliabiligy and power of heteroskedasticity tests. Econometric Reviews, 18(2):169–194, 1999.
• John (1971) S. John. Some optimal multivariate tests. Biometrika, 58(1):123–127, 1971.
• Koenker (1981) R. Koenker. A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1):107–112, 1981.
• Koenker and Bassett (1982) R. Koenker and G. Bassett. Robust tests for heteroscedasticity based on regression quantiles. Econometrica, 50(1):43–61, 1982.
• Lee (1992) B. J. Lee. A heteroskedasticity test robust to conditional mean misspecification. Econometrica, 60(1):159–171, 1992.
• Mauchly (1940) J. W. Mauchly. Significance test for sphericity of a normal n-variate distribution. The Annals of Mathematical Statistics, 11(2):204–209, 1940.
• Muirhead (1982) R. J. Muirhead. Aspects of Multivariate Statistical Theory. Wiley, 1982.
• Newey and Powell (1987) W. K. Newey and J. L. Powell. Asymmetric least squares estimation and testing. Econometrica, 55(4):819–847, 1987.
• Paul and Aue (2014) D. Paul and A. Aue. Random matrix theory in statistics: A review. Journal of Statistical Planning and Inference, 150:1–29, 2014.
• Réffy (2005) J. Réffy. Asymptotics of Random Unitaries. PhD thesis, Budapest University of Technology and Economics. http://www.math.bme.hu/~reffyj/disszer.pdf., 2005.
• Song and Gupta (1997) D. Song and A. K. Gupta. Lp-norm uniform distribution. Proceedings of the American Mathematical Society, 125(2):595–601, 1997.
• Su and Ullah (2013) L. Su and A. Ullah. A nonparametric goodness-of-fit-based test for conditional heteroskedasticity. Econometric Theory, 29(1):187–212, 2013.
• Taddy (2013) M. Taddy. Multinomial inverse regression for text analysis. Journal of the American Statistical Association, 108(503):755–770, 2013.
• White (1980) H. White. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4):817–838, 1980.
• Yao et al. (2015) J. Yao, S. Zheng, and Z. Bai. Large Sample Covariance Matrices and High-dimensional Data Analysis. Cambridge University Press, 2015.
• Yawitz and Marshall (1981) J. B. Yawitz and W. J. Marshall. Measuring the effect of callability on bond yields. Journal of Money, Credit and Banking, 13(1):60–71, 1981.

## Appendix A Technical proofs

According to (3), the OLS residuals are normal distributed , where is a projection matrix of rank . Let . Since has i.i.d. zero-mean normal variables, it is easily seen that has the same distribution as for any orthogonal matrix . Therefore is a -frame, that is, it is distributed as columns of a Haar matrix (Muirhead, 1982, Chapter 2). Furthermore, since , if we complement to an orthogonal matrix , we have then and becomes a -frame () distributed as columns of a Haar matrix. Therefore, we have

 ^ε=UU′ε=UZ, (11)

where under the null hypothesis. Notice that despite the multiplication by , is indepedent of (since its conditional distribution given is independent of ). Rewrite as

 (12)

Then the components (residuals) of can be expressed as

 ^εi=v′iZ=k∑j=1vijzj,for i=1,…,n. (13)

The proofs below rely on precise properties of the -frame ( columns of a Haar matrix). These useful properties are recalled in the next section, followed by the proofs of the main results of the paper.

### a.1 Haar matrix and related results

Here we present some important results of Haar matrix that will be used afterwards. First, the elements of in (12) have the same marginal distribution by symmetry and the square of each element has a beta distribution with parameter , see for example Réffy (2005). Their (marginal) moments are thus easily known. For example, we have

 E(v211)=1n;E(v411)=3n(n+2); E(v611)=15n(n+2)(n+4);E(v811)=105n(n+2)(n+4)(n+6). (14)

In addition, these elements are not independent, but weakly correlated, the moments of their products can be obtained using the following facts of an orthogonal matrix:

 .n∑i=1v2ij=1 1≤j≤k; .n∑i=1vijvij′=0,1≤j≠j′≤k.

Meanwhile, by Lemma 3.4 of Réffy (2005), for positive integers ,

if is odd for some , or is odd for some . This leads to the following list of cross-moment identities that will be used in upcoming proofs. The cross-moments of two elements in a same row (or same column) are as follows

 E(v211v212)=1n(n+1); E(v411v212)=3n(n+2)(n+4); E(v611v212)=15n(n+2)(n+4)(n+6); E(v411v412)=9n−6n(n+2)2(n+4)(n+6). (15)

The cross-moments of two elements in different rows and different columns are

 E(v411v222)=3(n+3)n(n−1)(n+2)(n+4); E(v411v422)=9n2+81n+222n(n−1)(n+2)2(n+4)(n+6). (16)

The cross-moments of three elements in a same row (or same column) are

 E(v211v212v213)=1n(n+2)(n+4); E(v411v212v213)=3(n2+4)n(n−2)(n+2)2(n+4)(n+6). (17)

The cross-moments of three elements in different rows or different columns are

 E(v211v212v222)=n+1n(n−1)(n+2)(n+4); E(v211v212v223)=n+3n(n−1)(n+2)(n+4); E(v411v221v222)=3n2+15n+42n(n−1)(n+2)2(n+4)(n+6); E(v411v222v223)=3n3+21n2+12n−156n(n−1)(n−2)(n+2)2(n+4)(n+6). (18)

The cross-moments of four elements in the same row (or same column) is

 E(v211v212v213v214)=n3−3n2−4n−60n(n−2)(n−3)(n+2)2(n+4)(n+6). (19)

The cross-moments of four elements in different rows or different columns are

 E(v211v212v221v222)=n3+3n2−4n−36n(n−1)(n−2)(n+2)2(n+4)(n+6); E(v211v212v221v223)=n4+3n3−10n2−36n+96n(n−1)(n2−4)2(n+4)(n+6); E(v311v12v21v22)=−3n(n−1)(n+2)(n+4); E(v211v212v223v224)=n4+5n3−10n2−44n+120n(n−1)(n2−4)2(n+4)(n+6); E(v311v12v321v22)=−9n−6n(n−1)(n+2)2(n+4)(n+6); E(v311v12v21v322)≈E(v311v12v321v22). (20)

The last approximate expression is due to the symmetry between the elements. Finally, some useful cross-moments of more than four elements in different rows or different columns are as follows:

 E(v211v12v13v22v23)=−1n(n−1)(n+2)(n+4); E(v311v12v21v22v223)=−3n2−6n−48n(n−1)(n−2)(n+2)2(n+4)(n+6); E(v211v12v13v221v22v23)=−n3−6n2+20n−48n(n−1)(n2−4)2(n+4)(n+6); E(v11v12v13v14v21v22v23v24)=3(n3−6n2+20n−48)n(n−1)(n−3)(n2−4)2(n+4)(n+6); E(v211v12v13v22v23v224)≈E(v211v12v13v221v22v23). (21)

Next, by Theorem 2.1 of Song and Gupta (1997), the joint distribution of all the squared elements in in (12) is known to be

 (v2i1,v2i2,…,v2ik)∼Dk(12,⋯,12;n−k2), (22)

where is the Dirichlet distribution with positive parameters . Therefore, has beta distribution with parameters . It follows that

 E(||vi||2)=cn,var(||vi||2)=2cn(1−cn)n+2, (23)
 cov(||vi||2,||vj||2)=2cn(cn−1)(n−1)(n+2),for i≠j, (24)
 E(log||vi||2)=logcn+1n−1k+O(1n2)+O(1k2), (25)
 E(||vi||2log||vi||2)=cn(logcn+1k−1n)+O(1n2)+O(1k2), (26)
 E(log||vi||2)2 = (logcn)2+2(1n−1k)logcn+2k−2n (27) +O(1n2)+O(1k2)+O(1nk),

Next, we derive the asymptotic limits for some joint distributions of .

###### Lemma 3.

Based on the above results on , as , we have

 (28)
###### Proof.

For , the multivariate central limit theorem states that

 Σ−1/20⋅√n(||v1||2−cn||v2||2−cn)D⟶N(0,I2),

where

By the previous results (23) and (24), we obtain that

 n⋅var(||v1||2)=n⋅var(||v2||2)=2cn(1−cn), n⋅cov(||v1||2,||v2||2)=2cn(cn−1)n→0as n→∞.

Then, Lemma 3 follows. ∎

There are two corollaries (easy consequences) of (28) by delta method:

 √n⎛⎜ ⎜⎝(√2cn(1−cn))−1(||v1||2−cn)(√2(1−cn)/cn)−1(log(||v2||2)−logcn)⎞⎟ ⎟⎠D⟶N(0,I2),

and

 √ncn/2(1−cn)⎛⎜⎝log(||v1||2)−logcnlog(||v2||2)−