Abstract
Testing heteroscedasticity of the errors is a major challenge in highdimensional regressions where the number of covariates is large compared to the sample size. Traditional procedures such as the White and the BreuschPagan 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 lowdimensional setting where the number of variables is fixed while the sample size tends to infinity, and the proportional highdimensional 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 highdimensional 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, highdimensional regression, hypothesis testing, Haar matrix.
Testing for Heteroscedasticity in Highdimensional 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
(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
(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 highdimensional 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 . Highdimensional 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 highdimensional 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 squareroot Lasso methods, respectively. However, if the errors are homoscedastic, these heteroscedasticityconsistent methods may lose efficiency as suggested by the phenomenon arising in lowdimensional 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 Lassotype of methods result in biased estimates of the coefficients, and the least squares method is preferable to other Mestimators in highdimensional 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 highdimensional estimation methodologies highlights the importance of conducting heteroscedasticity detection as a preliminary step in practice in order to select a suitable estimation method for highdimensional regressions.
Heteroscedastic testing has been extensively studied for classical lowdimensional 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 highdimensional regressions. The large sample theory of all the existing tests discussed above is developed under the lowdimensional framework where the dimension should be fixed while the sample size tends to infinity. By referring to recent advances in highdimensional 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 lowdimensional 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 highdimensional effect is an interesting phenomenon shown in Ferrari and CribariNeto (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 highdimensional regression.
In this paper, we propose two new procedures for testing heteroscedasticity, which are dimensionproof in the sense that they are valid for a wide range of dimension (covering both low and highdimensional settings). More precisely, our procedures are theoretically valid once the degree of freedom is large enough (precisely when ). This includes for instance the lowdimensional setting where and the highdimensional 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 highdimensional regressions. More surprisingly, even in lowdimensional 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 setup 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 lowdimensional situation where is a constant and . Therefore, the procedure derived under this setting will be applicable to both the high and lowdimensional 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 likelihoodratio 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
Without assuming the homoscedasticity, the likelihood is maximised by solving the system of equations
Therefore, the maximum likelihood estimator (MLE) of satisfy the system of equations
The corresponding maximized likelihood is
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
On the other hand under the homoscedasticity hypothesis, the OLS estimator and the estimator of the variance
are in fact the MLEs. So the maximized likelihood under the null hypothesis is
(4) 
Therefore, the approximate likelihood ratio, likelihood ratio is first derived by Mauchly (1940), is defined as
where it is reminded that . This suggests to consider the approximate likelihoodratio statistic
(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 scalefree 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
(6) 
where is the Euler constant.
The testing procedure using with the critical value from (6) is referred as the approximate likelihoodratio test (ALRT). In addition to the scalefree 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 finitesample 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
(7) 
where
and
2.2 The coefficientofvariation 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 coefficientofvariation statistic
(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 scalefree 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
(9) 
The testing procedure using with the critical value from (9) is referred as the coefficientofvariation 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
(10) 
where
and
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 multinormal. 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%.
n=500  n=1000  

ALRT  CVT  White  BP  ALRT  CVT  White  BP  ALRT  CVT  White  BP  
0.05  4.62  4.24  5.92  3.74  4.88  5.66  0.16  4.64  5.16  5.48  NA  4.16 
0.1  5.00  4.64  0.28  3.80  4.86  5.78  NA  3.64  5.44  5.20  NA  3.60 
0.3  4.98  4.70  NA  1.66  4.60  6.06  NA  1.88  5.38  5.06  NA  2.32 
0.5  5.30  4.72  NA  0.52  4.84  4.80  NA  0.70  5.04  5.14  NA  0.72 
0.7  4.66  4.58  NA  0  5.60  5.50  NA  0.02  5.38  5.70  NA  0.02 
0.9  5.06  4.28  NA  0  5.48  5.24  NA  0  4.44  5.04  NA  0 
NA denotes Not Applicable”
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 dimensionsamplesize 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 24 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 highdimensional 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.
3.3 NonGaussian design
In this section, we investigate the performance of all tests applied to nonGaussian 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.
Gamma design  Uniform design  
ALRT  CVT  White  BP  ALRT  CVT  White  BP  
0.05  4.68  5.80  0.22  4.12  4.48  4.84  0.14  4.18 
0.1  4.94  4.80  NA  4.28  5.08  4.84  NA  3.96 
0.3  5.14  5.62  NA  2.42  5.02  4.72  NA  2.10 
0.5  5.60  5.76  NA  0.60  5.26  5.20  NA  0.68 
0.7  5.60  6.00  NA  0.08  4.86  4.86  NA  0.02 
0.9  5.72  6.20  NA  0  4.70  4.26  NA  0 
NA denotes Not Applicable”
Gamma design  Uniform design  

Setting  ALRT  CVT  BP  ALRT  CVT  BP  
S1  0.05  1  1  0.9448  1  1  0.9998 
0.1  1  1  0.8046  1  1  0.9754  
0.3  1  1  0.4534  1  1  0.2964  
0.5  1  1  0.2278  1  1  0.0412  
0.7  0.9818  1  0.0304  0.9576  1  0.0032  
0.9  0.2760  0.9978  0  0.2356  0.9408  0  
S2  0.05  1  1  0.0366  1  1  0.0374 
0.1  1  1  0.0404  1  1  0.0374  
0.3  0.9222  0.9952  0.0230  0.9238  0.9942  0.0244  
0.5  0.4550  0.8270  0.0056  0.4412  0.8046  0.0052  
0.7  0.1410  0.3236  0.0004  0.1272  0.2914  0.0002  
0.9  0.0582  0.0812  0  0.0542  0.0630  0  
S3  0.05  1  1  1  0.6688  0.9138  0.9998 
0.1  0.9992  1  0.9998  0.3910  0.7128  0.9202  
0.3  0.3648  0.7522  0.5480  0.1006  0.1902  0.0804  
0.5  0.1144  0.2454  0.0352  0.0712  0.0866  0.0116  
0.7  0.0624  0.0946  0.0014  0.0542  0.0588  0.0002  
0.9  0.0472  0.0664  0  0.0484  0.0456  0 
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 nonnormal 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 lowdimensional scenario. The DM test is compared here (notice that this test is not in Tables 14 since its implementation in a multivariate setting is unclear). Following the same setup 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 highdimensional framework, it is still a competitive procedure in classical lowdimensional regression even with a small sample size. This is also supported by the results for the cases in Tables 14.
Settings  

Setting  ALRT  CVT  DM  ALRT  CVT  DM  
S1  0  0.047  0.032  0.056  0.046  0.020  0.053 
0.5  0.057  0.072  0.084  0.051  0.036  0.072  
1.0  0.131  0.213  0.148  0.096  0.093  0.089  
S2  0  0.047  0.032  0.053  0.046  0.020  0.052 
0.5  0.601  0.485  0.276  0.292  0.206  0.101  
1.0  0.884  0.792  0.365  0.570  0.406  0.094  
S3  0  0.047  0.032  0.054  0.046  0.020  0.053 
0.5  0.094  0.135  0.113  0.077  0.061  0.076  
1.0  0.250  0.331  0.198  0.152  0.145  0.114 
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 Lowdimensional data sets
In order to check the performance of the proposed tests in lowdimensional situation, we analyse two data sets: the ‘bond yield’ data and the ‘currency substitution’ data^{1}^{1}1These 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 longterm 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 CookHendershott model, we got three pvalues of 0.5614 (BP), 0.3307 (ALRT) and 0.8333 (CVT). And the YawitzMarshall model yields three pvalues 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 pvalue of 0.5779 and the CVT test reports a pvalue of 0.1309. However, the pvalue 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’ data^{2}^{2}2Available on the website: https://stuff.mit.edu/ vchern/NBER/ Belloni et al. (2011) and the ‘eminentdomain’ data^{3}^{3}3Available on the website: 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 pvalues 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 pvalue of 0.0686 (compared with the value 0.5822 in full model). However, the pvalue 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 threejudge panels on judicial decisions, the data set ‘eminentdomain’ 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 pvalue of and , respectively, strongly supporting these authors’ approval. On the other hand, the BP test cannot detect the existence of heteroscedasticity by reporting a pvalue 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 lowdimensional scheme, not only cannot possess a correct size, but also loses power when heteroscedasticity exists.
5 Conclusion and discussion
For highdimensional 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 highdimensional data sets. Extensive MonteCarlo 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 highdimensional tests when used in the tested lowdimensional situations also show a performance comparable to that of the existing classical tests which are designed specifically under lowdimensional 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 nonGaussian setting. Although we have shown some robustness of the proposed procedures against nonGaussian design in simulation experiments, a thorough investigation is missing. It is however clear that new theoretical tools will be needed to tackle with such nonGaussian 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 wellknown 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 mestimation in highdimensional 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 highdimensional 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. Highdimensional 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 squareroot 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. Highdimensional 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 highdimensional 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 CribariNeto (2002) S. L. P. Ferrari and F. CribariNeto. 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 nvariate 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. Lpnorm uniform distribution. Proceedings of the American Mathematical Society, 125(2):595–601, 1997.
 Su and Ullah (2013) L. Su and A. Ullah. A nonparametric goodnessoffitbased 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 heteroskedasticityconsistent 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 Highdimensional 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. zeromean 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
(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
(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
(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:
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 crossmoment identities that will be used in upcoming proofs. The crossmoments of two elements in a same row (or same column) are as follows
(15) 
The crossmoments of two elements in different rows and different columns are
(16) 
The crossmoments of three elements in a same row (or same column) are
(17) 
The crossmoments of three elements in different rows or different columns are
(18) 
The crossmoments of four elements in the same row (or same column) is
(19) 
The crossmoments of four elements in different rows or different columns are
(20) 
The last approximate expression is due to the symmetry between the elements. Finally, some useful crossmoments of more than four elements in different rows or different columns are as follows:
(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
(22) 
where is the Dirichlet distribution with positive parameters . Therefore, has beta distribution with parameters . It follows that
(23) 
(24) 
(25) 
(26) 
(27)  
Next, we derive the asymptotic limits for some joint distributions of .
Lemma 3.
Based on the above results on , as , we have
(28) 