1 INTRODUCTION
In many modern economic applications, datasets often contain a large number of covariates p relative to the sample size n. When performing high-dimensional regression on such data, parameter estimation typically assumes homoscedastic errors for simplicity. However, heteroscedasticity is common in practice, and ignoring it can lead to inefficient estimates and invalid inference. Despite its practical importance, testing for heteroscedasticity remains both underexplored and challenging in high-dimensional settings.
Classical heteroscedasticity tests, such as those proposed in Breusch and Pagan (Reference Breusch and Pagan1979), White (Reference White1980), and Koenker and Bassett (Reference Koenker and Bassett1982), were developed under the assumption of fixed p, and are no longer valid when
$p> n$
, which is the primary focus of this article. To address this gap, we propose a high-dimensional heteroscedasticity test in the proportional regime where
$n/p \to \delta \in (0,1]$
. Our approach extends the classical test of Newey and Powell (Reference Newey and Powell1987), which considers the following heteroscedastic linear model:
$$ \begin{align} Y = X\beta_0 + \sigma(X \gamma_0)\varepsilon, \end{align} $$
where
$X = (X_1, \ldots , X_p) \in \mathbb {R}^p$
,
$\beta _0, \gamma _0 \in \mathbb {R}^p$
are unknown parameter vectors, and
$Y, \varepsilon \in \mathbb {R}$
, with
$\varepsilon $
independent of X. The function
$\sigma (\cdot )$
captures heteroscedasticity, allowing the error variance to vary with the linear combination
$X\gamma _0$
.
To detect nonconstant variance, we use the expectile as a distributional statistic, evaluated at various levels
$\tau \in (0,1)$
, defined as
$$ \begin{align} u_\tau(Y) - \mathbb{E}(Y) = \frac{2\tau - 1}{1 - \tau} \int_{[u_\tau(Y), \infty)} (y - u_\tau(Y)) \, dF_Y(y). \end{align} $$
In model (1.1), the
$\tau $
th expectile of Y given X is
$u_\tau (Y \mid X) = X\beta _0 + u_\tau (\sigma (X\gamma _0)\varepsilon )$
. A key insight is that under heteroscedasticity, the population parameter
$\beta _0$
is affected by an additional contribution from
$\sigma (X\gamma _0)$
, scaled by the expectile of the error term
$\varepsilon $
. In contrast, under homoscedasticity,
$\beta _0$
remains invariant across different values of
$\tau $
. This observation becomes even clearer when a local linear approximation of
$\sigma (\cdot )$
in (2.3) is considered, and it motivates a testing framework based on distributional statistics (Koenker and Bassett, Reference Koenker and Bassett1982; Newey and Powell, Reference Newey and Powell1987), where the hypothesis concerns the influence of
$\gamma _0$
.
Compared to previous work (Newey and Powell, Reference Newey and Powell1987), this article addresses two additional challenges commonly encountered in contemporary datasets: (i) the number of parameters p exceeds the sample size n and diverges at a linear rate:
$n/p \rightarrow \delta \in (0,1]$
and (ii) the p-vectors
$\beta _0$
and
$\gamma _0$
are sparse. This linear rate regime was first discussed in the early approximate message passing (AMP) literature (Donoho, Maleki, and Montanari, Reference Donoho, Maleki and Montanari2009, Reference Donoho, Maleki and Montanari2010) and then formally stated in Bayati and Montanari (Reference Bayati and Montanari2011a, Thm. 1). The AMP algorithm has attracted significant attention for the theoretical analysis of regularized M-estimators and low-rank matrix estimation. Subsequent work often separately discusses
$\delta \in (0,1]$
for sparse regression with regularization (e.g., Bradic, Reference Bradic2016; Su, Bogdan, and Candès, Reference Su, Bogdan and Candès2017; Weinstein, Barber, and Candès, Reference Weinstein, Barber and Candès2017; Weinstein et al., Reference Weinstein, Su, Bogdan, Barber and Candès2023; Zhou, Claeskens, and Bradic, Reference Zhou, Claeskens and Bradic2020) and
$\delta \in (1, \infty )$
for unstructured regression coefficient vectors (Bradic, Reference Bradic2016; Donoho and Montanari, Reference Donoho and Montanari2016; Sur, Chen, and Candès, Reference Sur, Chen and Candès2019; Sur and Candès, Reference Sur and Candès2019; Candès and Sur, Reference Candès and Sur2020). We refer to Feng et al. (Reference Feng, Venkataramanan, Rush and Samworth2022) for a thorough literature review. In addition to the AMP framework, El Karoui et al. (Reference El Karoui, Bean, Bickel, Lim and Yu2013) and El Karoui (Reference El Karoui2018) introduced an alternative way of analyzing the estimators in this
$n, p$
proportional regime, with conclusions agreeing with the AMP analysis. Specifically, the analysis addresses unregularized and ridge regression, where a sparse structure is not assumed for the parameter vector. In this work, we focus on sparse regression and address the
$n/p \to \delta \in (0,1]$
setting.
In the past decade, fundamental research for high-dimensional linear regression has focused on sparse regularization (Tibshirani, Reference Tibshirani1996; Candès and Tao, Reference Candès and Tao2007; Zou, Reference Zou2006; Fan and Lv, Reference Fan and Lv2008; Zhang, Reference Zhang2010). We refer to Bühlmann and van De Geer (Reference Bühlmann and van De Geer2011) and Fan et al. (Reference Fan, Li, Zhang and Zou2020) for a systematic review of this topic. Most existing literature in this field assumes homoscedastic regression errors, that is,
$\sigma (\cdot )$
is a constant, as opposed to the fact that heteroscedasticity is commonly observed in real data, see, for example, Wang, Wu, and Li (Reference Wang, Wu and Li2012) and Daye, Chen, and Li (Reference Daye, Chen and Li2012) on genomic data and Belloni, Chernozhukov, and Wang (Reference Belloni, Chernozhukov and Wang2014) on economic data.
There are natural solutions to address heteroscedasticity in model estimation. Cattaneo, Jansson, and Newey (Reference Cattaneo, Jansson and Newey2018) and Jochmans (Reference Jochmans2022) addressed this by splitting the p covariate vector X into a fixed low-dimensional subvector
$X_s$
of main interest with unknown parameter
$\beta _s$
and a nuisance subvector
$X_{p-s}$
with nonnegligible size
$q =p-s$
. Their high dimensionality refers to a setting where
$n> p$
but p is considerably larger than in the classic setting, where
$p/n \to 0$
. The asymptotic normality is further derived for
$\hat \beta _s$
, and the main focus is on heteroscedastic-robust asymptotic covariance estimators. Addressed by Jochmans (Reference Jochmans2022, Sect. 3.2, 3.3), the consistency of the heteroscedasticity-robust covariance estimator requires the length of the nuisance parameter to satisfy
$\lim _n \sup q/n < 1$
, and handling the case of
$1/2 <\lim _n \sup q/n < 1$
is already challenging. Another approach considers replacing the least-squares loss function with asymmetric loss functions, such as the quantile loss functions (Koenker and Bassett, Reference Koenker and Bassett1978; Koenker, Reference Koenker2005) and the expectile loss functions (Newey and Powell, Reference Newey and Powell1987). A thorough review of quantile regression can be found in Koenker (Reference Koenker2017). High-dimensional quantile regression, addressing
$p> n$
setting, has been studied by several authors (Wu and Liu, Reference Wu and Liu2009; Belloni and Chernozhukov, Reference Belloni and Chernozhukov2011; Wang et al., Reference Wang, Wu and Li2012; Zheng, Gallagher, and Kulasekera, Reference Zheng, Gallagher and Kulasekera2013; Gu et al., Reference Gu, Fan, Kong, Ma and Zou2018; He et al., Reference He, Pan, Tan and Zhou2023). Gu and Zou (Reference Gu and Zou2016) studied high-dimensional expectile regression and proposed COSALES, which couples two expectile loss functions at expectile level
$0.5$
and another level of interest. By this construction, COSALES can simultaneously estimate parameter vectors
$\beta _0$
and
$\gamma _0$
in (2.3). Additionally, Zhou and Zou (Reference Zhou and Zou2023) proposed a cross-fitted residual regression for estimating heteroscedasticity in (2.3).
Despite recent efforts in estimating heteroscedasticity through methods, such as quantile regression, expectile regression, or cross-refitted regression, testing for the presence of heteroscedasticity remains largely unexplored in the high-dimensional literature. Since quantile and expectile loss functions are robust to heteroscedastic regression errors, using these two asymmetric loss functions to construct heteroscedastic tests is natural and thus taken in this article. Moreover, as previously mentioned, the Newey–Powell test is shown to be superior to the quantile-based test in low-dimensional settings (Newey and Powell, Reference Newey and Powell1987). Therefore, we focus on the high-dimensional Newey–Powell test in this article. In fact, the contrast between expectile regression and quantile regression is even more significant in high dimensions, owing to the smoother loss function in expectile regression. Gu and Zou (Reference Gu and Zou2016) show that fitting a high-dimensional expectile regression model is as simple as fitting Lasso regression. In contrast, serious efforts have been dedicated to handling the nonsmooth loss function in high-dimensional quantile regression, including ADMM (Gu et al., Reference Gu, Fan, Kong, Ma and Zou2018) and smoothing approximation (He et al., Reference He, Pan, Tan and Zhou2023).
Given the above considerations, we extend the Newey–Powell test (Newey and Powell, Reference Newey and Powell1987) to sparse high-dimensional settings where
$n/ p \to \delta \in (0, 1]$
,
$n, p\to \infty $
, and
$\ell _1$
-regularization or other sparse regularization is necessary to obtain good estimators. The heteroscedasticity of regression error is tested via
$H_{0,j}: \gamma _{0, j} = 0$
versus
$H_{1, j}: \gamma _{0, j} \neq 0$
,
$j = 1, \ldots , p$
. By using the AMP theory, we establish the joint asymptotic normality of multiple expectile estimators, which we then use to construct the test statistic. We analytically examine the size and power of the proposed test.
The rest of the article is organized as follows. In Section 2, we first introduce the notation and assumptions for the high-dimensional heteroscedastic linear model in (2.3), following the expectile regression in Section 2.1. The AMP algorithm is introduced in Section 3, and the main technical results are included in Section 3.3. Section 4 presents the proposed heteroscedasticity test. The analytical expression of the power function is presented in Theorem 2. We discuss a practically important decorrelation step for using AMP in applications in Section 5. Section 6 contains simulation studies and real data examples from Belloni et al. (Reference Belloni, Chernozhukov and Hansen2013) and Lan et al. (Reference Lan, Zhong, Li, Wang and Tsai2016). Section 7 concludes the article. The assumptions, proofs of lemmas and theorems, and auxiliary lemmas and definitions are provided in Appendixes A–C, respectively. Additional supporting results are included in the Supplementary Material (Zhou and Zou, Reference Zhou and Zou2025).
2 HETEROSCEDASTIC LINEAR REGRESSION
We restate the well-studied form of the heteroscedastic model from Koenker and Bassett (Reference Koenker and Bassett1982) and Newey and Powell (Reference Newey and Powell1987):
$$ \begin{align} Y = X\beta_0 + \big(1 + X\gamma_0\big) \varepsilon, \end{align} $$
where 
are unknown parameter vectors. In addition, we have n independent and identically distributed (i.i.d.) copies of the pair
$(Y, X)$
denoted by
$(Y_i, X_i), i= 1, \ldots , n$
. Further, we denote n-vectors
${\boldsymbol Y} = (Y_1, \ldots , Y_n)^\top $
,
${\boldsymbol \varepsilon } = (\varepsilon _{1}, \ldots , \varepsilon _{n})^\top $
; the design matrix
${\boldsymbol X} = (X_1, \ldots \, X_n)^\top $
. The components of the vector
${\boldsymbol \varepsilon }$
satisfy Assumption (A1).
-
(A1)
$\varepsilon _1, \ldots , \varepsilon _n$
are i.i.d. random variables with mean zero, a finite
$(2\kappa - 2)$
th moment for
$\kappa \ge 2$
.
The constant
$\kappa $
depends on the loss function used to estimate the unknown parameters
$\beta _0,\gamma _0$
; in our case, the expectile loss function in (2.6). This statement will become clear in Section 3.2 when we discuss the AMP in detail. Further, following Bayati and Montanari (Reference Bayati and Montanari2011a), the empirical measure of a p-vector is defined as putting
$1/p$
weight on each component of the vector. Then, the unknown parameter vector
$\gamma _0$
complies with Assumption (A2).
-
(A2) The empirical measure of the components of the p-vector
$\gamma _0$
, when
$p\to \infty $
, converges weakly to the probability measure of a random variable
$\Gamma _0$
with bounded
$(2\kappa -2)$
th moment for
$\kappa \ge 2$
.
This assumption precedes a standard Assumption (A3) on the parameter vector
$\beta _0$
in the AMP literature. In practice, the proposed test is valid for
$\beta $
with finite p components (see Section 4), and the hypothesis in (4.33) tests the components
$\gamma _j$
’s for finite p. We also demonstrate applying the proposed test to two real data examples with finite p in Section 6.3.
-
(A3) The empirical measure of the components of the p-vector
$\beta _0$
, when
$p\to \infty $
, converges weakly to the probability measure of a random variable
$B_0$
with bounded
$(2\kappa -2)$
th moment for
$\kappa \ge 2$
.
Notice
$X\gamma _0 = \sqrt {n} X \frac {\gamma _0}{\sqrt {n}},$
where the scaled design
$\sqrt {n} X $
and parameter vector
$\frac {\gamma _0}{\sqrt {n}}$
agree with the assumptions in Koenker and Bassett (Reference Koenker and Bassett1982). Further, the sequence
$\gamma _n = \frac {\gamma _0}{\sqrt {n}}$
restricts the calculation to local alternatives and guarantees well approximation of sequences of multiplicative heteroscedastic models.
Following Newey and Powell (Reference Newey and Powell1987, Thm. 1(iii)), the
$\tau $
th expectile of Y given X can be expressed as
$$ \begin{align} u_\tau(Y \mid X) = X\beta_0 + (1 + X\gamma_0)u_\tau(\varepsilon) = X \beta_\tau + u_\tau(\varepsilon), \end{align} $$
where
$\beta _\tau = \beta _0 + \gamma _0 u_\tau (\varepsilon )$
is the population parameter vector corresponding to the
$\tau $
th expectile of Y, approximated by
$X\beta _\tau $
. Here,
$u_\tau (\varepsilon )$
is the
$\tau $
th expectile of the error term, as defined in (1.2), and will be denoted simply as
$u_\tau $
henceforth. The parameter
$\beta _\tau $
is fixed and depends on the scalar
$u_\tau $
. In this spirit, we rewrite the model as follows:
$$ \begin{align} Y = X \beta + \varepsilon, \end{align} $$
where
$\beta = \beta _0 + \varepsilon \gamma _0$
and the error
$\varepsilon $
satisfies Assumption (A1). In addition, we assume the ratio
$n / p \to (0, 1]$
when
$n, p \to \infty $
. Let
$s = \|\beta \|_0$
, where
$\|\cdot \|_0$
denotes the
$\ell _0$
-norm of a vector. Obviously,
$s \le \|\beta _0\|_0+\|\gamma _0\|_0$
. The way of parameterization in (2.5) has two reasons: (1) we want to express the model to fit it in the AMP framework for asymptotic analysis, where an independent error term stands alone and (2) the error term
$\varepsilon $
determines the location change of
$\beta $
from
$\beta _0$
. Identifying
$\gamma _0$
relies on making use of different spots of the error distribution and requires a pair of regularized estimators. This is technically challenging due not only to the nontrivial asymptotic behavior of regularized estimators when
$p>n$
, but also to the correlation structure among pairs of estimators whose population parameters do not have the same center.
2.1 Expectile Regression under Heteroscedasticity
The expectile regression (Newey and Powell, Reference Newey and Powell1987) estimates the
$\tau $
th conditional expectile of Y by minimizing an asymmetric quadratic loss function defined as
$$ \begin{align} \rho_{\tau}(x; u_\tau) =|\tau - I\{x \leq u_\tau\}| (x - u_\tau)^2 , \end{align} $$
where
$\tau \in (0, 1)$
is the expectile level and
$u_\tau $
denotes the
$\tau $
th expectile of error
$\varepsilon $
. which is the solution to the nonlinear equation (Newey and Powell, Reference Newey and Powell1987, Eq. (2.7)). A special case is
$\tau = 1/2$
, where the expectile loss function reduces to the least-squares loss. Under our model, the conditional expectile of Y at expectile level
$\tau $
follows
$X(\beta _0 + \gamma _0 u_\tau ) + u_\tau $
(Newey and Powell, Reference Newey and Powell1987) and the coefficient vector
$\beta _\tau = \beta _0 + \gamma _0 u_\tau $
, where
$u_\tau $
is the expectile of error
$\varepsilon $
at
$\tau $
in (1.2). Since
$\varepsilon $
has mean zero (Assumption (A1)), expectile regression at
$\tau = 1/2$
does not identify the combination of
$\gamma _0$
and
$\beta _0$
. Instead,
$\beta _{1/2}$
only estimates
$\beta _0$
.
The proposed heteroscedasticity test is constructed based on the estimator of the parameter vector
$\beta $
via minimizing the empirical expectile loss function. Specifically, we consider the
$\ell _1$
-regularized expectile estimator defined as follows:
where
$\|\cdot \|_1$
denotes the
$\ell _1$
-norm and
$\lambda $
denotes the tuning parameter. Due to technical restrictions of the AMP algorithm (see Section 3), we do not jointly optimize the intercept term
$u_\tau $
and the slope vector
$\beta $
. Instead, we treat the term
$u_\tau $
as a parameter of the expectile loss function. Recentering Y by the intercept
$u_\tau $
is important to guarantee the following AMP analysis when considering robust loss functions. In our simplest homoscedastic error case, plugging in
$u_\tau $
reduces the conditional
$\tau $
-mean in (2.4) to only
$X\beta _\tau $
. This intercept term is only ignorable when we address the conditional mean of Y given X by Assumption (A1), the zero mean of
$\varepsilon $
. In homoscedastic settings, the error expectile
$u_\tau $
can be estimated by solving
$$ \begin{align} \hat u_\tau - \frac{1}{n}\sum_{i = 1}^n r_i = [(2\tau - 1)/(1 - \tau)] \frac{1}{n}\sum_{i = 1}^n (r_i - \hat u_\tau) I\{r _i \ge \hat u_\tau \}, \end{align} $$
where
$r_i = Y_i - \widehat \beta _\tau X_i$
’s are obtained using an initial estimator
$\widehat \beta ,$
such as the regularized least-squares estimators, the regularized Huber estimator, etc. However, estimating
$u_\tau $
becomes challenging under heteroscedasticity since solving (2.8) estimates the
$\tau $
’s expectile of
$\sigma (X\gamma _0) \varepsilon $
instead of
$\varepsilon $
. We observe that
$u_\tau $
is reflected in the intercept term of
$\gamma _0u_\tau $
in the composition
$\beta _\tau = \beta _0 + \gamma _0 u_\tau $
. Since the coupled sparse asymmetric least squares (COSALES) (Gu and Zou, Reference Gu and Zou2016) estimates
$\beta _0$
and
$\gamma _0 u_\tau $
simultaneously for any
$\tau \in (0, 1), \tau \neq 1/2$
, we use the intercept estimator of
$\gamma _0 u_\tau $
from the COSALES as a surrogate for the true
$u_\tau $
in (1.2). Further, to evaluate the estimation error of
$u_\tau $
on the test, we compared the numerical results using the true and estimated
$u_\tau $
. We observe that using the estimated
$u_\tau $
from the COSALES does not affect the test results. We also investigated the performance of the test when plugging in (2.8) under the null in Table S.1 in the Supplementary Material (Zhou and Zou, Reference Zhou and Zou2025). The numerical performance appears similar.
Importantly, when the regression errors are i.i.d., that is,
$\sigma (\cdot ) = 1$
, the slope
$\beta $
does not vary across different expectile levels in the classical setting with fixed p and
$n\to \infty $
. We propose to test the violation of homoscedasticity under diverging
$n, p$
based on the ground of common slope under homoscedasticity at different expectile levels
$\tau _k, k = 1, \ldots , K$
.
3 APPROXIMATE MESSAGE PASSING FOR EXPECTILE REGRESSION
The AMP analysis starts by regularizing the loss functions
$\rho $
by using the proximal mapping operator defined as
where the parameter
$b>0$
is a positive scalar that controls the step size moving toward the minimum of the loss function
$\rho $
. In our case, the proximal operator for expectile loss (2.6) can be obtained following the derivation in Section B.1 and can be expressed as follows:
$$ \begin{align} \text{ Prox }_{\tau}(z;b) &= \left\{\begin{array}{@{}ll} \big(z+ 2b (1-\tau)u_\tau \big)\big({2b(1-\tau) + 1}\big)^{-1} , & z\leq u_\tau \nonumber\\[4pt] \big({z+ 2b \tau u_\tau}\big)\big({2b\tau +1 }\big), & z> u_\tau \end{array}\right. \\ &= \frac{z - 2bu_\tau (\tau - I\{z \leq u_\tau\})}{1- 2b(\tau - I\{z \leq u_\tau\})}. \end{align} $$
With the analytical expressions of the proximal operator for the expectile loss function (3.10), we further work out the expression of the effective score function (Donoho and Montanari, Reference Donoho and Montanari2016) defined as
$$ \begin{align} \widetilde G(z;b) = b \partial\rho(x)|_{x = \mathrm{Prox}(z;b)}, \text{ with } b>0, \end{align} $$
where
$\partial \rho $
is the subgradient of
$\rho $
defined as
$\partial \rho (x) = \{y: \rho (u)\ge \rho (x) + y(u-x), \text {for all } u \}$
. This leads to the subgradient of the expectile loss function in (3.12) as
$$ \begin{align} \partial\rho_{\tau}(z) = 2 |I\{z\leq u_\tau\} - \tau| (z- u_\tau). \end{align} $$
One could interpret the effective score function by associating it with the score function. By plugging (3.10) in (3.12) and using the definition of the effective score function in (3.11), we obtain the effective score function as follows:
$$ \begin{align} \widetilde G_{\tau}(z;b) = \left\{\begin{array}{@{}ll} \frac{2b(1-\tau)}{2b(1-\tau) + 1 } (z - u_\tau) , & z \leq u_\tau \\[4pt] \frac{2b\tau}{2b\tau + 1} (z- u_\tau) , & z>u_\tau. \end{array}\right. \end{align} $$
To analyze regularized estimators under a sparse structure, specifically, when
$ s/p \to \omega $
, where s denotes the number of nonzero components of
$\beta $
, Bradic (Reference Bradic2016) and Zhou et al. (Reference Zhou, Claeskens and Bradic2020) incorporated sparsity by introducing the rescaled effective score function
$$ \begin{align} G(z; b) = \delta \omega^{-1} \widetilde{G}(z; b). \end{align} $$
3.1 The Approximate Message Passing Algorithm
We consider the AMP algorithm to approximate the estimators in (2.7) in the sense that the averaged difference between the estimators obtained from the AMP algorithm and the estimators by minimizing the loss functions converges to zero almost surely, see Zhou et al. (Reference Zhou, Claeskens and Bradic2020, p. 2562, Eq. (23) and below) for a summary of discussions in the relevant references (Bayati and Montanari, Reference Bayati and Montanari2011b; Bradic, Reference Bradic2016; Donoho and Montanari, Reference Donoho and Montanari2016; Sur et al., Reference Sur, Chen and Candès2019; Huang, Reference Huang2020).
The AMP algorithm is an iterative algorithm indexed by t consisting of three steps that belong to the general recursion in Bayati and Montanari (Reference Bayati and Montanari2011a). We will describe the general recursion when stating the theory in Section 3.3. Let the matrix
${\boldsymbol X}$
satisfy Assumption (A4).
-
(A4) A standard Gaussian design:
$X_i, i = 1, \ldots , n,$
are i.i.d. copies of X, where with
$X_{j}\sim N(0,1/n), j=1,\ldots ,p,$
i.i.d. components.
with 
Since
$\beta _\tau $
is a sparse vector with most components being zero, the algorithm starts with an initialization
$\widehat \beta _{(0)} = 0, t = 0$
and proceed at
$t \ge 1$
following steps:
-
Step 1. Adjust residuals The adjusted residuals
$z_{(t)}$
are obtained by (3.15)where G is defined in (3.14) and
$\eta (\cdot; \theta _{(t-1)})$
is the soft-thresholding operator with parameter
$\theta _{(t-1)}$
defined as
$\eta (\beta; \theta ) = \text { sgn }(\beta )\max (|\beta | - \theta ,0)$
.
-
Step 2. Update the parameter for effective score function The scalar parameter
$b_{(t)}$
is chosen such that the average of the effective score function has slope 1, that is,
$ \frac {1}{n} \sum _{i = 1}^n \partial _1 G (z_{i, (t)}; b_{(t)}) = 1 $
, where
$\partial _1 G(z; b)$
denotes the derivative of G w.r.t. the first argument z in case of differentiable function G. Hence, for the expectile regression with calculations in Section B.2, the parameter
$b_{(t)}$
is updated by solving (3.16)where
$$ \begin{align} \frac{\omega}{\delta} = \frac{2b(1-\tau)}{2b(1-\tau) + 1 } \frac{\#\{i: z_{(t), i} \leq u_\tau\}}{n} + \frac{2b\tau}{2b\tau + 1 } \frac{\#\{i: z_{(t), i}> u_\tau\}}{n}, \end{align} $$
$\omega $
is estimated by the fraction of the nonnull components in an initial estimator
$\widehat \beta $
, and
$u_\tau $
is replaced by plugging in an estimator from the COSALES algorithm (see discussions under (2.8)). The update
$b_{(t)}$
at step t is taken to be the first value on a fine grid that solves (3.16).
-
Step 3. Update the estimators for
$\beta $
Update the coefficient estimator for
$\beta $
as follows: (3.17)
$$ \begin{align} \widehat\beta_{(t+1)} = \eta (\widetilde\beta_{(t)}; \theta_{(t)}), \text{ where } \widetilde\beta_{(t)} = \widehat\beta_{(t)} + {{\boldsymbol X}}^\top G(z_{(t)}; b_{(t)}). \end{align} $$
The estimator
$\widehat \beta $
is an approximation for regularized estimators in (2.7), where the soft-thresholding operator
$\eta (\cdot; \theta )$
is the proximal operator for the regularizer
$\|\beta \|_1$
. The tuning parameter
$\theta $
is chosen to minimize the mean squared error (MSE) of
$\widehat \beta $
. As suggested in Donoho and Johnstone (Reference Donoho and Johnstone1994, Reference Donoho and Johnstone1998), Donoho et al. (Reference Donoho, Maleki and Montanari2009), and Bayati and Montanari (Reference Bayati and Montanari2011b), we take
$\theta _{(t)} = \alpha \bar \zeta _{\mathrm {emp}(t)}$
for scalar
$\alpha $
and (3.18)
$$ \begin{align} \bar\zeta^2_{\mathrm{emp}, (t)} = \frac{1}{n} \sum_{i = 1}^n G(z_{i, (t)}; b_{(t)})^2. \end{align} $$
In practice, one chooses a grid of values for
$\alpha $
, then obtains the corresponding estimator
$\widehat \beta $
and the estimated MSE for each value of
$\alpha $
on the grid. The optimal
$\alpha $
is chosen to be the one with the minimum estimated MSE. The details can be found in Zhou et al. (Reference Zhou, Claeskens and Bradic2020, Algorithm 1). We also implement practical stopping rules: the maximum number of iterations is set to T = 100, and convergence is declared when
$\text { MSE }(\widehat {\beta }_{(t+1)}, \widehat {\beta }_{(t)}) = \frac {1}{p} \sum _{j=1}^p (\widehat {\beta }_{j,(t+1)} - \widehat {\beta }_{j,(t)})^2 < \text { tol }$
. In the simulations, the tolerance level
$\text { tol }$
is set to
$10^{-8}$
for high sparsity settings and
$10^{-7}$
for medium sparsity settings.Details of the population version of the parameter
$\bar \zeta _{\mathrm {emp}, (t)}^2$
can be found in Section 3.2. The estimator
$\widetilde \beta $
is of major interest in this project since it is connected to the debiased estimator (Javanmard and Montanari, Reference Javanmard and Montanari2014; van de Geer et al., Reference van de Geer, Bühlmann, Ritov and Dezeure2014) and follows a Gaussian distribution asymptotically (Bayati and Montanari, Reference Bayati and Montanari2011b; Bayati, Erdogdu, and Montanari, Reference Bayati, Erdogdu and Montanari2013). We will discuss this in Section 3.3.
3.2 State Evolution
When
$p, n \to \infty $
, the averaged-over-components performance of estimators
$\widetilde \beta _{(t)}$
and
$\widehat \beta _{(t)}$
from the AMP algorithm can be described using two important parameters
$\bar \zeta _{(t)}^2$
,
$\bar \sigma _{(t)}^2$
. Assume that the following limit exists
$\bar \sigma _{(0)}^2 = \lim _{p\to \infty }(p \delta )^{-1} \|q_{(0)}\|^2 < \infty $
. With such initialization, the state evolution recursion for 
is defined as
$$ \begin{align} \bar\zeta_{(t)}^2 = \mathbb E [G(\varepsilon + \bar\sigma_{(t)}Z; b_{(t)})^2], \end{align} $$
$$ \begin{align} \bar\sigma_{(t)}^2 = \delta^{-1} \mathbb E[ (\eta(B + \bar\zeta_{(t-1)}Z; \theta_{(t - 1)} ) - B )^2], \end{align} $$
where B denotes the limiting form of the components of
$\beta $
in (2.5) and
$Z \sim N(0, 1)$
is a standard Gaussian random variable that is independent of B and
$\varepsilon $
. The parameter
$\bar \zeta _{(t)}$
is the limiting version of (3.18) when
$n, p \to \infty $
(Zhou et al., Reference Zhou, Claeskens and Bradic2020, Sect. B.2.3).
The readers might wonder if (3.19) and (3.20) still hold under heteroscedasticity since the relevant literature has not discussed this issue and often assumes homoscedasticity. It is worth mentioning that the state evolution recursion in (3.19) and (3.20) is valid as long as the empirical distribution of
$\varepsilon $
and the components of
$\beta $
converge weakly (Bayati and Montanari, Reference Bayati and Montanari2011a, Thm. 2). We will discuss this in detail in Section 3.3.
Using the two parameters, we can analyze the limiting performance of the estimator
$\widehat \beta $
from the AMP algorithm. A typical example is analyzing the asymptotic MSE (AMSE), defined as
$\mathrm {AMSE}(\widehat \beta _{(t+1)}, \beta ) = \lim _{p\to \infty } \frac {1}{p}\sum _{j = 1}^p(\widehat \beta _{(t+1),j} - \beta _j)^2$
. It is easy to show that AMSE complies with
$$ \begin{align*} \mathrm{AMSE}(\widehat\beta_{(t+1)}, \beta) = \lim_{p \to \infty} \frac{1}{p} \sum_{j = 1}^p \Big(\eta(\widetilde\beta_{(t),j}; \theta_{(t)}) - \beta_{j}\Big)^2 \stackrel{a.s.}{=} \mathbb E[\{\eta(B + \bar\zeta_{(t)}Z; \theta_{(t}) ) - B \}^2], \end{align*} $$
which corresponds to Bradic (Reference Bradic2016, Eq. (3.4)) with
$\widetilde \beta _{(t),j}$
in (3.17). An important observation from (3.20) is that, at each iteration step t, the estimator
$\widetilde \beta _{(t),j}$
is obtained by adding a Gaussian noise
$N(0, \bar \zeta ^2_{(t)})$
to the jth component of the true parameter vector
$\beta _{ j}$
. This suggests that
$\widetilde \beta _{(t)} \sim N(\beta , \bar \zeta ^2_{(t)}I_{p\times p})$
for finite marginals. Similar results are also obtained for the Lasso estimator in Bayati and Montanari (Reference Bayati and Montanari2011b), which is a special case of the
$\ell _1$
-regularized expectile regression when
$\tau = 1/2$
. We will come back to the asymptotic normality of
$\widetilde \beta _{(t)}$
in Section 3.3 after introducing the main technical lemma for multiple estimators
$\widetilde \beta _{k, (t)}, k= 1,1 \ldots , K$
.
3.3 Asymptotic Analysis for Multiple Estimators
The heteroscedasticity test for expectile regression is based on the difference of the estimators at different expectile levels
$\tau _k, k = 1, \ldots , K$
. The test statistic is based on asymptotic analysis for multiple estimators
$\widehat \beta _{k}$
’s among which the correlation is nontrivial when estimated using the same dataset (Zhou et al., Reference Zhou, Claeskens and Bradic2020, Lemma 1). Since estimating the correlation requires information on the unknown parameter vector
$\beta $
, we further extend the discussion on the correlations to the adjusted residuals
$z_{k, (t)}$
’s in the main technical result—Lemma 1 and construct an adjusted-residual-based estimator for the covariance between
$\widetilde \beta _{k}$
’s. Further, we want to address that the technical results, including the general recursion and Lemma 1 in this section, hold for similar loss functions, such as the quantile loss function.
We describe the general recursions in Bayati and Montanari (Reference Bayati and Montanari2011a) for the completeness of the article and a better connection of Lemma 1 and the AMP algorithm in Section 3.1. Under the heteroscedastic model (2.5), the coefficient vector
$\beta _k = \beta _0 + u_k \gamma _0$
, where
$u_k$
is a constant featuring the strength of heteroscedasticity which depends on
$\varepsilon $
. For expectile regression,
$u_k$
is the error expectile. By Assumptions (A2) and (A3), the empirical distribution of the parameter vector
$\beta _k$
converges weakly to
$B_k = B_0 + u_k \Gamma _0$
. In cases where
$u_k = 0$
or homoscedasticity such that
$\Gamma _0 = 0$
degenerates,
$B_k = B_0$
.
In addition, let a noise vector 
and the coefficient vector 
satisfy Assumptions (A1) and (A3). By fixing the initial condition
$q_{k, (0)}$
and recall the convention
$m_{k,(-1)} = 0$
, the general recursion for
$t\ge 0$
is defined as
$$ \begin{align} &&h_{k, (t+1)} = {{\boldsymbol X}}^\top m_{k, (t)} - \xi_{1,k, (t)}q_{k, (t)},\ \ \ m_{k, (t)} = g_{1,k, (t)}(d_{k, (t)}, {\boldsymbol \varepsilon}) , \nonumber\\ &&d_{k, (t)} = {{\boldsymbol X}} q_{k, (t)} - \xi_{2,k,(t)} m_{k, (t-1)},\ \ \ q_{k,(t)} = g_{2,k,(t)}(h_{k,(t)}, \beta_k) , \end{align} $$
where
$$ \begin{align} \xi_{1,k, (t)} = n^{-1} \sum_{i = 1}^n \partial_1 g_{1,k,(t)}(d_{k,(t), i}, \varepsilon_i), \end{align} $$
$$ \begin{align} \xi_{2, (t)} = (\delta p)^{-1}\sum_{j = 1}^p\partial_1 g_{2,k,(t)}(h_{k,(t), j}, \beta_{k,j}). \end{align} $$
As a special case, the AMP algorithm takes
$$ \begin{align*} &h_{k, (t+1)} = \beta_{k} - {{\boldsymbol X}}^\top G_k(z_{k, (t)}; b_{k, (t)}) - \widehat\beta_{k, (t)}, \\ &q_{k, (t)} = \widehat\beta_{k, (t)} - \beta_k, \qquad z_{k,(t)} = {\boldsymbol \varepsilon} - d_{k, (t)}, \end{align*} $$
which indicates
$$ \begin{align} d_{k,(t)} = {\boldsymbol \varepsilon} - z_{k,(t)}, \qquad m_{k,(t)} = - {G_k}(z_{k, (t)}; b_{k, (t)}), \end{align} $$
with the functions
$g_{1,k,(t)}(x_1,x_2) = - G_k(x_2-x_1; b_{k,(t)})$
and
$g_{2, k, (t)}(x_1, \beta ) = \eta (\beta - x_1; \theta _{k, (t)}) - \beta $
. Notice
$\widehat \beta _{k, (t)}$
the regularized estimator in (3.17), and
$$ \begin{align} h_{k, (t+1)} = \beta_k - \widetilde\beta_{k, (t)}. \end{align} $$
This detail is mathematically trivial but allows us to apply Lemma 1 directly to obtain the asymptotic normality of any pairs
$(\widetilde \beta _{k_1, (t)}, \widetilde \beta _{k_2,(t)})$
, as well as later obtain the power function for the proposed hypothesis test in Theorem 2.
Before presenting the main technical lemma, we assume a moment condition (Assumption (A5)), which also complements Assumptions (A1)–(A3). Recall in Assumptions (A1)–(A3), we state that the empirical measures of the components of the vectors
$\beta _0$
,
$\gamma _0$
, and
${\boldsymbol \varepsilon }$
converge weakly to the probability measures of random variables
$B_0$
,
$\Gamma _0$
, and
$\varepsilon $
, respectively. Let
$\widehat {f}_{\beta _0(p)}$
,
$\widehat {f}_{\gamma _0(p)}$
, and
$\widehat {f}_{{\boldsymbol \varepsilon }(p)}$
denote the empirical measures of
$\beta _0$
,
$\gamma _0$
, and
${\boldsymbol \varepsilon }$
, that is, putting equal weights on the components of the vectors. Specifically, for the sequence of vectors
${\boldsymbol \varepsilon }(p)$
, we index the sequence by p to unify the notation following Bayati and Montanari (Reference Bayati and Montanari2011a). This sequence index is legit since our asymptotic theory is derived under
$n/p \to \delta \in (0, 1]$
. Then, Assumption (A5) can be stated as follows:
-
(A5) For some
$\kappa>1 $
, let the empirical measures
$\widehat {f}_{\beta _0(p)}$
,
$\widehat {f}_{\gamma _0(p)}$
, and
$\widehat {f}_{{\boldsymbol \varepsilon }(p)}$
converge weakly to the probability measures
$f_{B_0}$
,
$f_{\Gamma _0}$
, and
$f_\varepsilon $
. In addition, assume the following
$(2\kappa - 2)$
th moment conditions hold:-
(a)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{\beta _0(p)}}(B_0^{2\kappa - 2}) = \mathbb E_{f_{B_0}}(B_0^{2\kappa - 2}) < \infty $
; -
(b)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{\gamma _0(p)}}(\Gamma _0^{2\kappa - 2}) = \mathbb E_{f_{\Gamma _0}}(\Gamma _0^{2\kappa - 2}) < \infty $
; -
(c)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{{\boldsymbol \varepsilon } (p)}}(\varepsilon ^{2\kappa - 2}) = \mathbb E_{f_{\varepsilon }}(\varepsilon ^{2\kappa - 2}) < \infty $
; -
(d)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{q_{k}(p)}}(B_k^{2\kappa - 2}) < \infty $
.
-
The constant
$\kappa $
depends on the pseudo-Lipschitz order of function
$g_{1, k, (t)}$
, see Definition 1 for the definition of pseudo-Lipschitz functions and Lemma 1 for associating the order
$\kappa $
with the function
$g_1 = \tilde \psi _{c'}, \tilde \psi _{c"}$
. As described above, in the special case where the general recursion (3.21) is the AMP algorithm, function
$g_1$
takes a more specific form, that is, the negative effective score function
$-G_k(\cdot; \cdot )$
.
Lemma 1. Let the sequences of design matrices
$\{{\boldsymbol X}(p)\}$
, location coefficient vectors
$\{\beta _0(p)\}$
, heteroscedasticity parameter vectors
$\{\gamma _0(p)\}$
, error vectors
$\{{\boldsymbol \varepsilon }(p)\}$
, and initial condition vectors
$\{q_k(p)\}$
be weakly convergent sequences for K recursions satisfying Assumptions (A1)–(A5). Let
$\{\bar \sigma _{k, (t)}^2, \bar \zeta _{k, (t)}^2\}$
be defined uniquely by the recursions in (3.19) and (3.20). These are the state evolution parameters for the kth estimation with initialization
$\bar \sigma ^2_{k, (0)} = \lim _{n\to \infty }(p\delta )^{-1} \|q_{k,(0)}\|^2$
. Further, we assume a constant
$u_k$
for each of the K recursions specifying the signal strength of the heteroscedasticity parameter, such that the parameter vector (2.5) has the composition
$\beta _k = \beta _0 + u_k \gamma _0$
. Assume Lemma 1 in Bayati and Montanari (Reference Bayati and Montanari2011a) holds individually for each of the K recursions. Additionally, for all pseudo-Lipschitz functions 
of order
$\kappa _{\mathrm {c}}$
for some
$1 \leq \kappa _{\mathrm {c}} = \kappa /2$
with
$\kappa $
as in (A5) and the iteration index 
,
$$ \begin{align} &{\lim_{p\to \infty} \frac{1}{p}\sum_{j = 1}^p \tilde\psi_{\mathrm{c}}(h_{k_1, (1), j},\ldots, h_{k_1, (t+1), j}, \beta_{k_1, j}) \tilde\psi_{\mathrm{c}}(h_{k_2, (1), j}, \ldots, h_{k_2, (t+1), j}, \beta_{k_2, j})} \stackrel{a.s.}{=} \nonumber\\ &\mathbb E[\tilde\psi_{\mathrm{c}}(\bar\zeta_{k_1, (0)}Z_{k_1, (0)},\ldots, \bar\zeta_{k_1, (t)}Z_{k_1, (t)}, B_{k_1}) \tilde\psi_{\mathrm{c}}(\bar\zeta_{k_2, (0)}Z_{k_2, (0)},\ldots, \bar\zeta_{k_2, (t)}Z_{k_2, (t)}, B_{k_2})]. \nonumber\\ \end{align} $$
And, for all pseudo-Lipschitz functions 
of order
$\kappa _{\mathrm {c'}}, \kappa _{\mathrm {c"}}$
for some
$1 \leq \kappa _{\mathrm {c'}}, \kappa _{\mathrm {c"}} \leq \kappa /2$
, where
$\kappa _c' + \kappa _c" = \kappa $
with
$\kappa $
as in (A5),
$$ \begin{align} &{\lim_{n\to \infty} \frac{1}{n}\sum_{i = 1}^n \tilde\psi_{\mathrm{c'}}(d_{k_1, (0), i},\ldots, d_{k_1,(t), i}, \varepsilon_i) \tilde\psi_{\mathrm{c"}}(d_{k_2, (0), i}, \ldots, d_{k_2, (t), i}, \varepsilon_i)} \stackrel{a.s.}{=} \nonumber\\ &\mathbb E[\tilde\psi_{\mathrm{c'}}(\bar\sigma_{k_1, (0)}\hat Z_{k_1, (0)},\ldots, \bar\sigma_{k_1, (t)}\hat Z_{k_1, (t)}, \varepsilon) \tilde\psi_{\mathrm{c"}}(\bar\sigma_{k_2, (0)}\hat Z_{k_2, (0)},\ldots, \bar\sigma_{k_2, (t)}\hat Z_{k_2, (t)}, \varepsilon)], \nonumber\\ \end{align} $$
where
$(Z_{k, (0)}, \ldots , Z_{k, (t)}), (\hat Z_{k, (0)}, \ldots , \hat Z_{k, (t)})$
are multivariate Gaussian with,
$k = k_1, k_2$
, is a
$(t+1)$
-dimensional zero-mean multivariate standard normal vector independent of
$B_k$
,
$\varepsilon $
; at iteration t,
$(Z_{k_1, (t)}, Z_{k_2, (t)})$
and
$(\hat Z_{k_1, (t)}, \hat Z_{k_2, (t)})$
are standard normal 2-vectors with covariance not necessarily equal to zero.
Further, for
$0<t_1, t_2<t$
, with
$\stackrel {P}{\to }$
denoting convergence in probability, it holds that
$$ \begin{align} \frac{1}{n} \sum_{i = 1}^n m_{k_1, (t_1), i} m_{k_2, (t_2), i} \stackrel{P}{\to}\lim_{p\to \infty} \frac{1}{p}\sum_{j = 1}^p h_{k_1, (t_1+1), j} h_{k_2, (t_2+1), j}. \end{align} $$
Lemma 1 generalizes Bayati and Montanari (Reference Bayati and Montanari2011a, Lemma 1(b)) by considering the correlation structure between multiple recursions. Further, Lemma 1 extends the analysis to account for error heteroscedasticity, which has not been addressed in the relevant literature. When
$k_1 = k_2$
and the heteroscedastic parameter vector
$\gamma _0 = {\boldsymbol 0}$
is a null vector, Lemma 1 reduces to Bayati and Montanari (Reference Bayati and Montanari2011a, Lemma 1(b)). Further, the covariance matrices of the standard normal 2-vectors
$(Z_{k_1, (t)}, Z_{k_2, (t)})$
and
$(\hat Z_{k_1, (t)}, \hat Z_{k_2, (t)})$
will be involved in constructing the test statistic in Section 4. The proof of Lemma 1 is provided in the Supplementary Material (Zhou and Zou, Reference Zhou and Zou2025).
Remark 1. We thank a reviewer for drawing our attention to Theorem 1 and Lemma 2 in Javanmard and Montanari (Reference Javanmard and Montanari2013). Their result agrees with our conclusion in (3.26), which shows the asymptotic approximation of the covariance of the structure of multiple estimators. Our additional contributions lie in establishing (3.27) and (3.28), which are of independent interest to statistics and econometrics, as they are important for establishing Theorem 1, (3.30), and demonstrate that the correlation between estimators is linked to the correlation of residuals.
To guarantee an estimable test statistic, we propose an adjusted-residual-based estimator in Theorem 1, which is directly obtainable from the AMP algorithm. The construction of the covariance matrix estimators is included in Section B.4 by directly applying (3.28). It is worth noting that (3.28) is a contribution to the existing AMP literature. The significance of (3.28) is that it shows the correlation between estimators (noticing
$h_{k,(t+1)} = \beta - \widetilde \beta _{k, (t)}$
in (3.25)) from two sequences is asymptotically equivalent to the product of two effective score functions. This avoids the problem of lacking information about the true parameter vector
$\beta $
. Additionally, Theorem 1 validates that
$Z_{k_1}$
and
$Z_{k_2}$
are correlated by (3.29).
Theorem 1. Under the assumptions in Lemma 1 and
$\kappa = 2$
for any two recursions
$k_1, k_2 = 1, \ldots , K$
, with probability 1, it holds that
$$ \begin{align} \lim_{p\to \infty} \frac{1}{p}\sum_{j = 1}^p \big(\widetilde\beta_{k_1, (t),j} - \beta_{k_1, j}\big)\big(\widetilde\beta_{k_2, (t), j} - \beta_{k_2, j}\big) = \bar\zeta_{k_1, (t)} \bar\zeta_{k_2, (t)}\mathrm{Cov}(Z_{k_1} ,Z_{k_2}), \end{align} $$
The estimator constructed as follows:
$$ \begin{align} \widehat\zeta_{k_1, k_2, (t)} = \frac{1}{n}\sum_{i=1}^n G_{k_1}(z_{k_1,(t), i}; b_{k_1, (t)}) G_{k_2}(z_{k_2, (t), i}; b_{k_2, (t)}) \end{align} $$
is a consistent estimator of the covariance
$\bar \zeta _{k_1, (t)}\bar \zeta _{k_2, (t)} \mathrm {Cov}(Z_{k_1, (t)}, Z_{k_2, (t)})$
.
Theorem 1 first validates our claim in Lemma 1 that
$Z_{k_1, (t)}, Z_{k_2, (t)}$
are correlated. The correlation between
$Z_{k_1, (t)}, Z_{k_2, (t)}$
describes the averaged correlation between the components of
$\widetilde \beta _{k_1, (t)}$
and
$\widetilde \beta _{k_2, (t)}$
. A similar conclusion holds for
$\hat Z_{k_1, (t)}, \hat Z_{k_2, (t)}$
. Further, the covariance between
$Z_{k_1, (t)}, Z_{k_2, (t)}$
can be approximated by using the adjusted residuals without any information about
$\beta $
.
The last piece of information before moving to the test statistic is the joint asymptotic distribution of the pair
$(\widetilde \beta _{k_1}, \widetilde \beta _{k_2})$
.
Lemma 2 (Joint normality)
Under the assumptions in Lemma 1 for any two recursions
$k_1, k_2 = 1, \ldots , K$
and for any pseudo-Lipschitz function 
of order
$\kappa _c$
in (3.26), the following convergence for estimators from the AMP algorithm holds with probability 1:
$$ \begin{align} \lim_{p\to \infty}\frac{1}{p}\sum_{j = 1}^p \psi(\widetilde\beta_{k_1, (t), j})\psi(\widetilde\beta_{k_2, (t), j}) = \mathbb E\big[\psi(B_{k_1}+ \bar\zeta_{k_1, (t)}Z_{k_1}) \psi(B_{k_2} + \bar\zeta_{k_2, (t)} Z_{k_2})\big], \end{align} $$
where
$(Z_{k_1}, Z_{k_2})$
is a standard normal random 2-vector.
The proof of Lemma 2 is a direct application of Lemma 1 (see Section B.3). When writing in the matrix form, the conditional joint asymptotic normality
${\boldsymbol {\tilde {\xi }}} - {\boldsymbol \xi } \sim N\Big ((u_{k_1}, u_{k_2}) \otimes \gamma _0 , \Sigma _{(t)} \otimes {\boldsymbol I}_{p\times p}\Big )$
holds for
$k_1, k_2 = 1, \ldots , K$
, where
${\boldsymbol {\tilde {\xi }}} = (\widetilde \beta _{k_1, (t)},\widetilde {\beta }_{k_2, (t)})$
,
${\boldsymbol \xi } = (1,1)\otimes \beta _0$
, and
${\boldsymbol I}_{p\times p}$
is a
$p\times p$
identity matrix. The covariance matrix
$\Sigma _{(t)}$
has
$(k_1, k_2)$
th component
$\Sigma _{k_1, k_2, (t)} = \bar \zeta _{k_1, k_2, (t)}$
, which can be estimated by (3.31).
Under homoscedasticity where
$\gamma _0 = {\boldsymbol 0}$
,
$\widetilde \beta _{k_1}$
and
$\widetilde \beta _{k_2}$
have a common center
$\beta _0$
. In addition, variance components
$\bar \zeta _{k, (t)}^2, k = 1, \ldots , K,$
are inflated under heteroscedasticity. This phenomenon can be observed by the inflation of the standard deviation
$\bar \sigma _{k, (t)}^2$
in (3.19), which is caused by the extra heteroscedastic parameter in (3.20).
4 HETEROSCEDASTICITY TEST
Although the asymptotic normality of
${\boldsymbol {\widetilde {\xi }}}$
holds for all iterations
$t = 0, 1, \ldots $
of the AMP algorithm, we use the estimators
${\boldsymbol {\widetilde {\xi }}}$
at
$t \to \infty $
to construct test statistics. When the iteration
$t \to \infty $
, the
$l_2$
-norm of the difference of the AMP output and the estimator in (2.7) converges to zero almost surely (see Bradic Reference Bradic2016, Thm. 2). In practice, we consider a stopping rule and define the algorithm’s convergence in Section 3.1 under (3.18). At the stop, we obtain an estimator approximating the estimator defined in (2.7). Hence, the iteration index t is dropped in the following sections. We consider testing the components of the vector
$\gamma _0$
$$ \begin{align} H_{0, j}: \gamma_{0, j} = 0 \text{ versus } H_{1, j}: \gamma_{0, j} \neq 0, j = 1, \ldots, p. \end{align} $$
Lemma 2 states that the K expectile estimators follow a multivariate Gaussian distribution asymptotically with a common center
$\beta _0$
under homoscedasticity. And the rows of
${\boldsymbol {\widetilde {\xi }}}$
are i.i.d. following a joint Gaussian distribution with covariance matrix
$\Sigma $
.
For
$\gamma _{0, j}$
, the test statistic is defined as
$$ \begin{align} T_j = \frac{\Delta {\boldsymbol {\tilde{\xi}}}_j^\top }{\sqrt{\Delta \Sigma \Delta^\top}}, j = 1, \ldots, p. \end{align} $$
Following Newey and Powell (Reference Newey and Powell1987, Sect. 4.1), we consider the contrast between two expectile levels for the testing problem, that is,
${\boldsymbol {\tilde {\xi }}} = (\widetilde \beta _{\tau _1}, \widetilde \beta _{\tau _2})$
,
$\Delta = (1, -1)$
, and the components of the covariance matrix
$\Sigma $
are replaced by their estimators in (3.31). Then, by Lemma 2, the test statistic
$T_j$
follows a standard normal distribution
$N(0, 1)$
. For the two-sided alternative in (4.33), the p-value can be calculated by
$$ \begin{align} P_{j} = 2\left\{1 - \Phi \left(|\frac{\Delta {\boldsymbol {\tilde{\xi}}}_j^\top }{\sqrt{\Delta \Sigma \Delta^\top}} |\right) \right\}, \end{align} $$
with rejection region
$R = (-\infty , \Phi ^{-1}(\alpha /2)) \cup (\Phi ^{-1}(1 - \alpha /2), \infty )$
, where
$\Phi $
denotes the cumulative distribution function of a standard normal distribution
$N(0,1)$
. The null hypothesis is rejected if
$P_j \le \alpha $
at nominal level
$\alpha $
. To analyze the power function on the population scale when
$n, p \to \infty $
, we consider for any jth component
$\vartheta (\gamma _{0, j}) = I\big ( T_j \in R, \Gamma _0 = \gamma _{0, j}\big )$
. Then, the limiting expression of the function
$\vartheta $
is obtained in Theorem 2 with proof in Appendix B.5.
Theorem 2. The function denoted by
$\vartheta (\gamma _{0, j}), j = 1, \ldots , p$
complies with
$$ \begin{align} &\lim_{p\to\infty}\frac{1}{p} \sum_{j = 1}^p \vartheta(\gamma_{0, j}) \stackrel{a.s.}{=} \\ &\mathbb E_{\Gamma_0}\Big[1 - \Phi\big(\Phi^{-1}(1 - \alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\big) + \Phi\big(\Phi^{-1}(\alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\big) \Big].\nonumber \end{align} $$
Theorem 2 can be used to confirm the size and power of the proposed test. Assume the nominal significance level of the test is
$\alpha $
. The size of the test
$ P(T_j \in R | \Gamma _0 = 0)$
is approximated by
$$ \begin{align*} &\frac{\lim_{p\to \infty}\frac{1}{p}\sum_{j = 1}^p I(T_j \in R ,\Gamma_0 = 0)}{P(\Gamma_0 = 0)}\\ &\stackrel{a.s.}{=} \frac{\mathbb E_{\Gamma_0}\Big[1 - \Phi\big(\Phi^{-1}(1 - \alpha/2) - 0\big) + \Phi\big(\Phi^{-1}(\alpha/2) - 0\big) \mid \Gamma_0 = 0\Big] P(\Gamma_0 = 0)}{P(\Gamma_0 = 0)}\\ &= \alpha, \end{align*} $$
where the conditional expectation is obtained by (4.36) at
$\Gamma _0 = 0$
. This confirms that the size of a nominal
$\alpha $
test is
$\alpha $
asymptotically.
Further, the empirical size of the test is defined in (6.37) and used to evaluate the performance of the proposed test. As presented in Tables 1 and 2, the empirical size, averaged over the components of
$\gamma _0$
, is approximately equal to the nominal level
$\alpha $
. Similarly, the power of the test is obtained by considering the nonnull components
$\gamma _{0,j}$
’s. However, from the right-hand side of (4.36), we notice that the power of the test depends on the averaged signal strength of
$\Gamma _0$
, the expectiles
$(u_{\tau _1}, u_{\tau _2})$
, and the covariance matrix
$\Sigma $
. The empirical power of the test defined in (6.38) is also used to evaluate the numerical performance of the test in Section 6.
Test (
$\alpha = 0.05$
) results for different distributions of
$\varepsilon $
under homoscedastic and heteroscedastic variance
Note: The pair of expectile levels is
$(\tau _1, \tau _2)= (0.2,0.8)$
. The top half presents results for
$\varepsilon $
following
$N(0,1)$
,
$t_3$
,
$\text {Lognormal}(0,1 )$
; the bottom half for
$0.9N(-0.2, 0.25) + 0.1N(1.8, 0.01)$
(mixNormal 1),
$0.9N(0.2, 0.25) + 0.1N(-1.8, 0.01)$
(mixNormal 2),
$0.95N(0, 0.25) + 0.05 N(0, 4)$
(mixNormal 3). Each simulation setting is replicated
$R = 400$
times. The average FP proportions are calculated for homoscedastic variance, where the heteroscedastic parameter vector
$\gamma _0 = {\boldsymbol 0}$
. For heteroscedastic variance, where
$\gamma _0$
has five nonnull components
$(3, 1, -5, -5, -3)$
, the average FP and TP proportions are calculated.
Test (
$\alpha = 0.05$
) results for different distributions of
$\varepsilon $
under homoscedastic and heteroscedastic variance
Note: The pair of expectile levels is
$(\tau _1, \tau _2)= (0.6,0.8)$
. The top half presents results for
$\varepsilon $
following
$N(0,1)$
,
$t_3$
,
$\text {Lognormal}(0,1 )$
; the bottom half for
$0.9N(-0.2, 0.25) + 0.1N(1.8, 0.01)$
(mixNormal 1),
$0.9N(0.2, 0.25) + 0.1N(-1.8, 0.01)$
(mixNormal 2),
$0.95N(0, 0.25) + 0.05 N(0, 4)$
(mixNormal 3). Each simulation setting is replicated
$R = 400$
times. The average FP proportions are calculated for homoscedastic variance, where the heteroscedastic parameter vector
$\gamma _0 = {\boldsymbol 0}$
. For heteroscedastic variance, where
$\gamma _0$
has five nonnull components
$(3, 1, -5, -5, -3)$
, the average FP and TP proportions are calculated.
5 DECORRELATION IN AMP
It is well known that the standard Gaussian design assumption, which is stated in Assumption (A4), is important for the AMP framework. Researchers have made several attempts to relax this assumption. One direction is to relax the Gaussian distribution assumption. For example, Bayati, Lelarge, and Montanari (Reference Bayati, Lelarge and Montanari2015) and Chen and Lam (Reference Chen and Lam2021) showed that the AMP framework could be extended to other random matrices; Rangan, Schniter, and Fletcher (Reference Rangan, Schniter and Fletcher2019) and Fan (Reference Fan2022) discussed the extension of AMP algorithms to right rotationally invariant design. Li and Sur (Reference Li and Sur2026) showed the statistical applications of these algorithms. However, the independence assumption is still essential in these works. A more relevant direction for this work is to relax the independence assumption in order to handle general design matrices with arbitrary covariance matrices. Consider a general design matrix in the sparse high-dimensional heteroscedastic linear model:
$ \tilde {{\boldsymbol Y}} = \tilde {{\boldsymbol X}} \beta + \tilde {{\boldsymbol \varepsilon }}, $
where 
, 
, 
, and the rows of
$\tilde {{\boldsymbol X}}$
are i.i.d. from a multivariate Gaussian distribution
$N_p(0, {\tilde \Sigma })$
. Several papers have discussed incorporating a general Gaussian design with strictly positive-definite covariance matrices. For example, Donoho and Montanari (Reference Donoho and Montanari2016, Def. 4.5, Cor. 4.6), Candès and Sur (Reference Candès and Sur2020, Sect. 3.1.1), Zhao et al. (Reference Zhao, Sur and Candès2022, Sect. 2), and Huang (Reference Huang2022) have proposed a method to decorrelate the design matrix using the nonsingular covariance matrix
${\boldsymbol {\tilde {\Sigma }}}$
by applying the transformation
$({\boldsymbol {\tilde X}\tilde {\Sigma }}^{-1/2})/\sqrt {n}$
. This decorrelation allows for satisfying Assumption (A4) and expands the applicability of AMP-based methods and theories. However, this decorrelation approach has limitations when dealing with sparse parameter vectors
$\beta $
because the null components of
$\beta $
become nonnull under the linear transformation
$\sqrt {n}{\boldsymbol {\tilde {\Sigma }}}^{1/2}\beta $
and the corresponding back transformation. Additionally, in practical applications, the covariance matrix
${\boldsymbol {\tilde {\Sigma }}}$
is often unknown, and estimating it becomes challenging when
$p>n$
unless certain structure assumptions are made. These challenges restrict the use of this approach in real applications.
In applications, we consider a different method for decorrelation in the AMP. The validity of this type of decorrelation has been tested numerically in Zhou and Claeskens (Reference Zhou and Claeskens2023, Table 3). We utilize a “generalized puffer transformation” on the data, resulting in transformed variables denoted as
${\boldsymbol {Y}} = F \tilde {{\boldsymbol Y}}$
and
${\boldsymbol {X}} = F \tilde {{\boldsymbol X}}$
, where the transformed data satisfy the model
$ {\boldsymbol Y}= {\boldsymbol X}\beta + F { \tilde {{\boldsymbol \varepsilon }}}. $
It is important to note that in our setting,
$p>n$
. To perform the transformation, we utilize matrices 
, 
, and 
obtained from the singular value decomposition
${\boldsymbol {\tilde {X}} = UDV^\top }$
. The “generalized puffer transformation” was originally proposed in Jia and Rohe (Reference Jia and Rohe2015), where the authors aimed to decorrelate the design matrix to bypass the so-called “irrepresentable condition” for the Lasso regression. The “generalized puffer transformation” uses
$F = \sqrt {p/n} {\boldsymbol U {\hat {D}} U}^\top $
, where
${\boldsymbol {\hat {D}}}$
is a diagonal matrix with diagonal elements
${\boldsymbol {\hat {D}}}_{ii} = I({\boldsymbol {D}}_{ii} \le 1/\sqrt {n})\sqrt {n}+I({{\boldsymbol D}}_{ii}> 1/\sqrt {n})/{{\boldsymbol D}}_{ii}, i=1,\ldots ,n$
. Subsequently, we apply the AMP algorithm and conduct the hypothesis test using the transformed data
$({\boldsymbol {Y}}, {\boldsymbol {X}})$
.
Unlike the decorrelation step in Donoho and Montanari (Reference Donoho and Montanari2016), Candès and Sur (Reference Candès and Sur2020), and Zhao et al. (Reference Zhao, Sur and Candès2022), our approach does not involve a transformation of
$\beta $
. Thus, it can preserve sparse structures. When
$p> n$
, we cannot guarantee that
${\boldsymbol X}$
precisely satisfies Assumption (A4), but the correlation is largely reduced, as observed in the numerical study. Moreover, on the transformed data, the error term becomes
$F \tilde {{\boldsymbol \varepsilon }}$
. Under the null hypothesis that the model (2.3) is homoscedastic, we have verified that
$F\tilde {{\boldsymbol \varepsilon }}$
is independent of the transformed design matrix
${\boldsymbol X} = \sqrt {p/n}{\boldsymbol U V^\top }$
. This observation suggests that the “generalized puffer transformation” does not alter the null hypothesis. This explains why the proposed test retains its validity size, as demonstrated in our simulation study. For further details, please refer to Section 6.2.
6 NUMERICAL PERFORMANCE
6.1 Simulation
We consider simulation settings where
$p = 500$
and
$n= 250$
. The samples
$X_i$
are generated from a standard Gaussian design, as specified in Assumption (A4). The nonnull components of the true parameter vector
$\beta _0$
are generated using a standard normal
$N(0, 1)$
distribution and then multiplied by 5. Each setting is replicated
$R = 400$
times. The parameter vector
$\beta _0$
and design matrix
${\boldsymbol X}$
are fixed across R replications using the same seed number, while the regression errors
${\boldsymbol \varepsilon }$
vary in each replication. We evaluate the test under six different error distributions:
$N(0,1 )$
,
$t_3$
,
$\text { Lognormal }(0,1)$
,
$0.9N(-0.2, 0.25) + 0.1N(1.8, 0.01)$
,
$0.9N(0.2, 0.25) + 0.1N(-1.8, 0.01)$
, and
$0.95 N(0, 0.25) + 0.05 N(0, 4)$
. The first three error distributions are randomly generated, centered, and scaled to have a standard deviation of 0.5. We aim to assess the test’s performance using symmetric errors (represented by
$N(0,1)$
), heavy-tailed errors (represented by
$t_3$
), and fat-tailed errors (represented by
$\text { Lognormal }(0,1)$
). Additionally, we consider three mixture normal distributed errors. The first two,
$0.9N(-2, 0.25) + 0.1N(18, 0.01)$
and
$0.9 N(2, 0.25) + 0.1N(-18, 0.01)$
, model clustered outliers on the right and left sides of the main error distribution, respectively. The last distribution—
$0.95 N(0, 0.25) + 0.05 N(0, 4)$
—adds additional variability to
$N(0, 0.25)$
. For heteroscedasticity, we assume five components of
$\gamma _0$
are nonnulls
$(3, 1, -5, -5, -3)^\top $
, and the rest components are all zeros. Similar to
$\beta _0$
, the components of
$\gamma _0$
are fixed throughout R replications using the same seed number.
The individual test is evaluated by the averaged size (or false positive) proportion and power (or true positive proportion). Both proportions are obtained by averaging over simulation replications and the components of
$\gamma _0$
, which agrees with the construction of the power function in Theorem 2. At significance level
$\alpha $
, the size can be estimated via
$$ \begin{align} \text{ FP }(\alpha) =(p-s)^{-1}\sum_{j \in S^{\mathrm{c}}}{\sum_{r=1}^R I\{P_{r,j} \leq \alpha \}}/{R} \qquad\text{ (empirical size),} \end{align} $$
and the power can be estimated by
$$ \begin{align} \text{ TP }(\alpha) =s^{-1}{\sum_{j \in S}{\sum_{r=1}^R I\{P_{r,j} \leq \alpha \}}/{R}} \qquad \text{ (empirical power)}. \end{align} $$
For homoscedastic cases, the components of
$\gamma _0$
are all zeros and
$\|\gamma _0\|_0 = 0$
; hence only FP values are calculated. For heteroscedasticity cases, five components of
$\gamma _0$
are nonnull with magnitude
$(3, 1, -5, -5, -3)^\top $
. The signal strength of heteroscedasticity
$u_\tau \gamma _0$
depends on the error expectile
$u_\tau $
. Since the errors are scaled to a standard deviation of 0.5, the heteroscedasticity signal strength is largely weakened by
$u_\tau $
.
We consider two sparsity levels for fixed ratio
$\delta = 0.5$
,
$p = 500$
, and
$n = 250$
.
-
1. (High-sparsity)
$\beta _0$
has
$s= 5$
nonnull components. -
2. (Medium-sparsity)
$\beta _0$
has
$s = 40$
nonnull components.
The above two settings intend to investigate the impact of sparsity. Further, we consider homoscedastic and heteroscedastic regression errors for each sparsity level, where the errors are generated from six different distributions. Simulation results are presented in Tables 1 and 2. The test statistic in (4.34) is constructed using pairs of expectile estimators
${\boldsymbol \xi } = (\widetilde \beta _{\tau _1}, \widetilde \beta _{\tau _2})$
, where seven pairs—
$\tau _2 = 0.8$
,
$\tau _1 = 0.1, 0.2, 0.6, $
and
$\tau _2= 0.9$
,
$\tau _1 = 0.1, 0.2, 0.6, 0.8$
—are tested. Since the test results for four pairs are similar, we only report two pairs—
$(\tau _1, \tau _2) = (0.2, 0.8)$
and
$(0.6, 0.8)$
. Table 1 uses expectile levels at
$(\tau _1, \tau _2)= (0.2, 0.8)$
, and Table 2 uses expectile levels at
$(\tau _1, \tau _2) = (0.6, 0.8)$
.
In addition, under homoscedasticity of the regression error, the test statistic in (4.34) can be seen as samples from a standard normal
$N(0,1)$
distribution. To confirm the theory, we also inspected the normality of the test statistic constructed using different pairs of expectile estimators. Example QQ-plots are included in Figure S.1 in the Supplementary Material (Zhou and Zou, Reference Zhou and Zou2025); the examples are randomly chosen among simulation replications to avoid “cherry-picking.”
We make several observations from Tables 1 and 2. When comparing Tables 1 and 2, we observe that using the expectile level pair
$(\tau _1, \tau _2) = (0.2, 0.8)$
leads to higher test power, which agrees with Theorem 2 that shows the power of the test depends on the difference
$(u_{\tau _2 } - u_{\tau _1})$
. Placing two expectile levels far enough would preserve the signal strength of
$\gamma _0$
to a certain extent. However, we do not suggest using two very extreme expectiles, since estimating them is more challenging than estimating those closer to the center.
When comparing the results for different error distributions, we note that the average
$\text {FP}$
proportions are close to the nominal
$\alpha $
, which complies with Theorem 2. The average
$\text {FP}$
proportions are slightly greater than the nominal level for the mixture normal distribution
$0.9N(-0.2, 0.25) + 0.1 N(1.8, 0.01)$
. At the same time, this error distribution tends to have the highest empirical power among all six error distributions. In medium-sparsity settings, where the number of nonnull components of
$\beta _0$
is 40, estimating
$\beta _{\tau }$
becomes more challenging, and it is not surprising that the power of the test in this setting dropped compared to Setting 2, where the number of nonnull components of
$\beta _0$
is 5. The estimation variances diminish the signal strength when
$s=40$
. This setting is very challenging for all tests, not just for ours.
6.2 Check the Decorrelation Proposal
We use simulations to check the validity of the decorrelation step proposed in Section 5. We assume an autoregressive (
$AR(1)$
) covariance matrix for
${\boldsymbol {\tilde {X}}} \sim N(0, {\boldsymbol {\tilde \Sigma }})$
, where the components
${\boldsymbol {\tilde \Sigma }}_{ij} = \rho ^{|i - j|}$
for
$i, j = 1, \ldots , p$
. The error distributions considered are
$N(0,1)$
and
$t_3$
. To assess the impact of correlation on test performance, we consider
$\rho = 0.3$
and
$0.7$
.
We see in Tables 3 and 4 that the test is still valid numerically using the transformed data by the “generalized puffer transformation.” We have checked that on the transformed data, the error is independent of the covariates, which suggests that the null hypothesis remains unchanged after the transformation. Simulation results demonstrate that the size of the test is approximately equal to the nominal
$\alpha $
. The power of the test decreases as the “generalized puffer transformation” reduces the signal-to-noise ratio. Furthermore, the power of the test decreases when the correlation parameter
$\rho $
increases from 0.3 to 0.7.
Test (
$\alpha = 0.05$
) results are presented for the decorrelation procedure under both homoscedastic and heteroscedastic variance
Note: The pair of expectile levels is
$(\tau _1, \tau _2) = (0.2, 0.8)$
. The covariance matrix for
$\boldsymbol {\tilde X} \sim N(0,\boldsymbol {\tilde \Sigma })$
has components
$\boldsymbol {\tilde \Sigma }_{ij} = \rho ^{|i - j|}$
for
$i, j = 1, \ldots , p$
. The error distributions considered are
$N(0,1)$
and
$t_3$
. Each simulation setting is replicated
$R = 400$
times. The average FP proportions are calculated for homoscedastic variance, where the heteroscedastic parameter vector
$\gamma _0 = \mathbf {0}$
. For heteroscedastic variance, where
$\gamma _0$
has five nonnull components
$(3, 1, -5, -5, -3)$
, the average FP and TP proportions are calculated.
Test (
$\alpha = 0.05$
) results are presented for the decorrelation procedure under both homoscedastic and heteroscedastic variance
Note: The pair of expectile levels is
$(\tau _1, \tau _2) = (0.6, 0.8)$
. The covariance matrix for
$\boldsymbol {\tilde X} \sim N(0, \boldsymbol {\tilde \Sigma })$
has components
$\boldsymbol {\tilde \Sigma }_{ij} = \rho ^{|i - j|}$
for
$i, j = 1, \ldots , p$
. The error distributions considered are
$N(0,1)$
and
$t_3$
. Each simulation setting is replicated
$R = 400$
times. The average FP proportions are calculated for homoscedastic variance, where the heteroscedastic parameter vector
$\gamma _0 = {\boldsymbol 0}$
. For heteroscedastic variance, where
$\gamma _0$
has five nonnull components
$(3, 1, -5, -5, -3)$
, the average FP and TP proportions are calculated.
6.3 Real Data Examples
6.3.1 International Economic Growth Data
We apply the proposed test on the “international economic growth” data used in Belloni et al. (Reference Belloni, Chernozhukov and Hansen2013), where the authors mainly focused on variable selection and estimation. We now want to test whether heteroscedasticity is present in this data. The dataset is available on the website https://stuff.mit.edu/~vchern/NBER/ and consists of
$n = 90$
samples. Further, we exclude the intercept term in the unprocessed dataset. The original dataset consists of 61 covariates measuring aspects, such as education, trading openness, and science policies. We further expand the covariate set by including the polynomial terms up to order 4 of all 61 covariates. The expanded dataset consists of 244 covariates with a ratio
$\delta $
approximately 0.37. The primary interest is the response variable—the national growth rate of GDP per capita. We apply the decorrelation procedure outlined in Section 5 before constructing the test. Similar to Section 6.1, the test statistic is calculated using two
$\ell _1$
-regularized expectile estimators. We consider expectile estimators at level
$\tau _2 = 0.8$
and at
$\tau _1 = 0.2$
. According to Belloni et al. (Reference Belloni, Chernozhukov and Hansen2013), the dataset is homoscedastic since no unmeasured underlying variables exist. Our test results show that none of the covariates cause heteroscedasticity across the different pairs. To assess the stability of the method when using pairs at different expectile levels, we also test
$\tau _1= 0.4, 0.5, 0.6$
. The test results agree with
$\tau _1 = 0.2$
for this dataset.
6.3.2 Supermarket Data
We also apply the proposed heteroscedasticity test to the “supermarket data” used in Lan et al. (Reference Lan, Zhong, Li, Wang and Tsai2016). The authors extended the t-test in their study and focused on testing
$\beta _{0, j} \neq 0$
in a high-dimensional setting (
$p \gg n$
). The dataset consists of
$n = 464$
daily records, where the response variable is the number of customers, and the covariates represent the sales volume of 6,398 different products. Although the dataset is potentially heteroscedastic, the test proposed in Lan et al. (Reference Lan, Zhong, Li, Wang and Tsai2016) is robust to heteroscedastic errors.
Similarly, we utilize the fused Kolmogorov filter (Mai and Zou, Reference Mai and Zou2015) to reduce the number of covariates to
$p = 928$
, achieving
$\delta = n/p=0.5$
. The ratio
$\delta =n/p= 0.5$
. Note that this reduced dimension is well beyond the theoretical limit
$O(n/\log (n))$
, which ensures the sure screening property of the fused Kolmogorov filter. Thus, we have high confidence that the screening step only removes noisy features and retains all important features. After applying the decorrelation procedure, we employ a pair of
$\ell _1$
-regularized expectile estimators at
$(\tau _1, \tau _2) = (0.2, 0.8)$
to construct the test. The test results indicate that 61 covariates potentially contribute to heteroscedasticity. We also tested
$\tau _1 = 0.4, 0.5, 0.6$
and
$\tau _2= 0.8$
. When using the pair
$(\tau _1, \tau _2) = (0.6, 0.8)$
, potentially 43 covariates contribute to heteroscedasticity. In addition, there are 34 common covariates detected by both
$\tau _1 = 0.2$
and
$\tau _1 = 0.6$
. The other two pairs using
$\tau _1 = 0.4, 0.6$
have more stable performance, overlapping 43 and 47 covariates with the pair using
$\tau _1 = 0.2$
.
7 CONCLUDING REMARKS
In this article, we propose a high-dimensional Newey–Powell heteroscedasticity test based on contrasting expectile regression at two different expectile levels. The testing procedure is based on a bias-corrected estimator directly obtained from the AMP algorithm in Section 3, which is tracked by the state evolution recursion asymptotically when
$n, p\to \infty $
. We clarify that we want to offer an alternative debiasing framework in addition to the existing one, which requires adding a bias-correction term to the regularized estimators. Such a well-explored debiasing framework requires estimating sparse sample precision matrices and a thorough theoretical investigation for the bias-corrected estimator, which requires a considerable amount of work. The proposed test has stable performance when
$\gamma _0$
is highly sparse, meaning that only a few covariates
$X_j$
contribute to the heteroscedasticity. The numerical performance suggests that the error distributions also have an impact on the size and power of the test, which is worth further investigation. Although the AMP theory is developed for multiple levels, we did not explore contrasting multiple expectiles for testing heteroscedasticity in this work. This could be interesting for future research.
Upon examining Theorem 2, we observe that the strength of the signal in
$\gamma _0$
and the covariance of the two expectile estimators also impact the power of the test. An open question remains regarding the optimal choice of expectile levels.
While the current literature on the AMP framework assumes an independent Gaussian design matrix, it is important to acknowledge that this assumption can be relaxed. However, a significant challenge in AMP remains unresolved: the lack of known theories to relax the assumption of independent covariates, accounting for both numerical and theoretical aspects. It is crucial to emphasize that this challenge is not unique to our work but is a broader issue in the AMP literature. To address this limitation in practical applications, we propose a “generalized puffer transformation” to enhance the applicability of the proposed test to real-world data. Through a simulation study, we demonstrate that the resulting test maintains its validity by preserving its nominal size. However, a further theoretical investigation is necessary to understand and address this challenge, and it remains an open question for future research.
APPENDIX
A ASSUMPTIONS
For the sake of completeness, we summarize the assumptions used in our theory.
-
(A1)
$\varepsilon _1, \ldots , \varepsilon _n$
are i.i.d. random variables with mean zero, a finite
$(2\kappa - 2)$
th moment for
$\kappa \ge 2$
. -
(A2) The empirical measure of the components of the p-vector
$\gamma _0$
, when
$p\to \infty $
, converges weakly to the probability measure of a random variable
$\Gamma _0$
with bounded
$(2\kappa -2)$
th moment for
$\kappa \ge 2$
. -
(A3) The empirical measure of the components of the p-vector
$\beta _0$
, when
$p\to \infty $
, converges weakly to the probability measure of a random variable
$B_0$
with bounded
$(2\kappa -2)$
th moment for
$\kappa \ge 2$
. -
(A4) A standard Gaussian design:
$X_i, i = 1, \ldots , n,$
are i.i.d. copies of with
$X_{j}\sim N(0,1/n), j=1,\ldots ,p,$
i.i.d. components. -
(A5) For some
$\kappa>1$
,-
(a)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{\beta _0(p)}}(B_0^{2\kappa - 2}) = \mathbb E_{f_{B_0}}(B_0^{2\kappa - 2}) < \infty $
; -
(b)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{\gamma _0(p)}}(\Gamma _0^{2\kappa - 2}) = \mathbb E_{f_{\Gamma _0}}(\Gamma _0^{2\kappa - 2}) < \infty $
; -
(c)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{{\boldsymbol \varepsilon } (p)}}(\varepsilon ^{2\kappa - 2}) = \mathbb E_{f_{\varepsilon }}(\varepsilon ^{2\kappa - 2}) < \infty $
; -
(d)
$\lim _{p \to \infty } \mathbb E_{\widehat {f}_{q_{k}(p)}}(B_k^{2\kappa - 2}) < \infty $
.
-
with 
B PROOFS
B.1 Proof of (3.10)
Proof. The proximal operator is defined as the minimizer of the objective function
$b\rho (x) + \frac {1}{2}(x - z)^2$
, which, by elementary calculus, is the root of the subgradient of the objective function. The proximal operator is then obtained by solving an equation of x as follows:
$$\begin{align*}b\partial\rho_{\tau}(x;u_\tau) + x - z = 0, \end{align*}$$
where
$\partial \rho _{\tau }$
is the gradient of the differentiable expectile loss function
$\rho _{\tau }$
in (3.12). Calculations follow similar steps and arguments for the two cases, that is,
$x \leq u_\tau $
and
$x> u_\tau $
. We present the calculation steps for
$x \leq u_\tau $
in detail. When
$x \leq u_\tau $
,
$\partial \rho _{\tau }(x;u_\tau ) = 2(1 - \tau ) (x - u_\tau )$
. Then, by plugging in
$\partial \rho _{\tau }$
and rearranging
$2b(1-\tau ) (x- u_\tau ) +x - z = 0$
, we obtain
$$\begin{align*}x = \frac{z + 2b (1 - \tau)u_\tau}{2b(1-\tau) + 1}.\end{align*}$$
The domain of z is obtained by
$\frac {z + 2b (1 - \tau )u_\tau }{2b(1-\tau ) + 1} \leq u_\tau $
which is
$z\leq u_\tau $
.
B.2 Proof of (3.16)
Proof. By (3.14), the subgradient of the effective score function
$\tilde G_{\tau }(z; b)$
in (3.13) w.r.t. z follows
$$\begin{align*}\partial_1 G_{\tau}(z; b) = \frac{\delta}{\omega}\Big\{\frac{2b(1-\tau)}{2b(1-\tau) + 1} I\{z \leq u_\tau\} + \frac{2b\tau}{2b\tau + 1} I\{z> u_\tau\}\Big\}. \end{align*}$$
The parameter b is updated by finding solutions
$b>0$
, s.t.
$$\begin{align*}1 = \frac{1}{n} \sum_{i = 1}^n \partial_1 G_{\tau}(z_i; b) = \frac{\delta}{\omega} \Big\{\frac{1}{n}\sum_{i = 1}^n \frac{2b(1-\tau)}{2b(1-\tau) + 1} I\{z_i \leq u_\tau\} + \frac{2b\tau}{2b\tau + 1} I\{z_i> u_\tau\}\Big\}. \end{align*}$$
B.3 Proof of Lemma 2
B.4 Proof of Theorem 1
We first prove (3.29) and (3.30). By choosing
$\widetilde \psi _c(x_1,\ldots , x_{t+1}, y) = - x$
in (3.26), we obtain
The last equality holds since
$(Z_{k_1, (t)}, Z_{k_2, (t)})$
is a standard normal random 2-vector. And we drop the iteration index t for convenience.
Similarly, we choose
$\tilde \psi _{c'}(x_{1}, \ldots , x_{t}, y)$
and
$\tilde \psi _{c"}(x_{1}, \ldots , x_{t}, y)$
in (3.27) to be the rescaled effective score functions
$G_{k_1}(y - x_{t}; \cdot )$
and
$G_{k_2}(y - x_{t}; \cdot )$
. By (3.24), we obtain
The last equality holds by
$Z \stackrel {d}{= } - Z$
for
$Z\sim N(0,1)$
.
To obtain (3.31), we first combine (3.24), (3.25), and (3.28) which leads to
Then, when
$n, p$
are sufficiently large, we use (3.31) to approximate the covariance of
$\widetilde \beta _{k_1, (t+1)}$
and
$\widetilde \beta _{k_2, (t+1)}$
.
B.5 Proof of Theorem 2
Proof.
$$ \begin{align*} &\vartheta(\gamma_{0, j}) = I\{T_j \in R_j, \Gamma_0 = \gamma_{0, j}\} \\ &= I \Big\{\frac{\Delta {\boldsymbol {\tilde{\xi}}}_j^\top }{\sqrt{\Delta \Sigma \Delta^\top}}> \Phi^{-1}(1 - \alpha/2)\Big\} + I\Big\{\frac{\Delta {\boldsymbol {\tilde{\xi}}}_j^\top }{\sqrt{\Delta \Sigma \Delta^\top}} < \Phi^{-1}(\alpha/2)\Big\}\\ &= I\Big\{\frac{\Delta ({\boldsymbol {\tilde{\xi}}}_j - (u_{\tau_1}, u_{\tau_2})\gamma_{0, j})^\top }{\sqrt{\Delta \Sigma \Delta^\top}}> \Phi^{-1}(1 - \alpha/2) - \gamma_{0, j}\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\Big\} \\ &\quad + I\Big\{\frac{\Delta({\boldsymbol {\tilde{\xi}}}_j- (u_{\tau_1}, u_{\tau_2})\gamma_{0, j})^\top }{\sqrt{\Delta \Sigma \Delta^\top}} < \Phi^{-1}(\alpha/2) - \gamma_{0, j}\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\Big\} .\\ \end{align*} $$
By Lemma 1 and define
$\widetilde \Xi = (B_0 + u_{\tau _1}\Gamma _0 + \bar \zeta _{1}Z_1, B_0 + u_{\tau _2}\Gamma _0 + \bar \zeta _{2}Z_2)$
, where
$(\bar \zeta _1 Z_1, \bar \zeta _2 Z_2) \sim N(0, \Sigma )$
, we obtain
$$ \begin{align*} &\lim_{p\to\infty}\frac{1}{p} \sum_{j = 1}^p \vartheta(\gamma_{0, j}) \\ &\stackrel{a.s.}{=} \mathbb E I\Big\{\frac{\Delta(\widetilde\Xi - (u_{\tau_1}, u_{\tau_2})\Gamma_0)^\top}{\sqrt{\Delta \Sigma \Delta^\top}}> \Phi^{-1}(1 - \alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}} \Big\}\\ &\quad + \mathbb E I\Big\{\frac{\Delta(\widetilde\Xi - (u_{\tau_1}, u_{\tau_2})\Gamma_0)^\top}{\sqrt{\Delta \Sigma \Delta^\top}} < \Phi^{-1}(\alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}} \Big\} \\ &= \mathbb E I \{ Z_3> \Phi^{-1}(1 - \alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}} \Big\} + \mathbb E I\{Z_3 < \Phi^{-1}(\alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}} \Big\} \\ &= \mathbb E_{\Gamma_0}\Big[1 - \Phi\big(\Phi^{-1}(1 - \alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\big) +\Phi\big(\Phi^{-1}(\alpha/2) - \Gamma_0\frac{\Delta (u_{\tau_1}, u_{\tau_2})^\top}{\sqrt{\Delta \Sigma \Delta^\top}}\big) \Big]. \end{align*} $$
C AUXILIARY DEFINITIONS AND LEMMAS
Definition 1 (Pseudo Lipschitz function)
A function 
is pseudo-Lipschitz of order
$\kappa \geq 1$
, if there exists a constant
$L> 0$
, such that for all 
,
$$ \begin{align*} |\phi(x) - \phi(y)| \leq L(1 + \| x\|^{\kappa - 1} + \|y\|^{\kappa -1})\|x - y\|. \end{align*} $$
It follows that if
$\phi $
is a pseudo Lipschitz function of order
$\kappa $
, then there exists a constant
$L'$
such that for all 
.
Lemma 3. For any pseudo Lipschitz functions 
of order
$\kappa _{c'}$
and
$\kappa _{c"}$
, their product function
$\phi _{c'}\phi _{c"}$
is a pseudo Lipschitz function of order
$\kappa _{c'} + \kappa _{c"}$
.
Lemma 4 (Bayati and Montanari (Reference Bayati and Montanari2011a, Lemma 2))
For any deterministic unit vectors 
and 
, and a random matrix
$\tilde {{\boldsymbol X}}$
distributed as 
, the following conclusions hold:
-
1.
$v^\top \tilde {{\boldsymbol X}} u \stackrel {d}{=} Z/\sqrt {n}$
, where
$Z \sim N(0,1)$
. -
2.
$\lim _{n\to \infty } \|\tilde {{\boldsymbol X}} u\|^2 = 1$
almost surely. -
3. Consider a d-dimensional subspace
with orthonormal basis
$s_1, \ldots , s_d$
and
$\|s_i\|^2 = n, i = 1, \ldots , d$
. Denote the matrix
${\boldsymbol D} = (s_1, \ldots , s_d)$
and the orthogonal projection onto as . Then, we have with satisfying
$\lim _{n\to \infty }\|x\| \stackrel {a.s.}{=}0$
and
$x = o_d(1)$
.
with orthonormal basis 
as 
. Then, we have 
with 
satisfying 
Lemma 5 (Jameson (Reference Jameson2014, Thm. 1))
If
$x_{i} \ge 0$
where
$i = 1, \ldots , n$
and
$p \ge 1$
, then
$ \sum _{i = 1}^n x_{i}^p \leq (\sum _{i = 1}^n x_{i})^p \leq n^{p-1} \sum _{i = 1}^n x_{i}^p$
. The reversed inequality holds for
$p \in (0, 1)$
.
Lemma 6 (Bayati and Montanari (Reference Bayati and Montanari2011a, Lemma 4))
Let
$\kappa \ge 2$
and a sequence of vectors
$\{v(N)\}_{N\ge 0}$
whose empirical distribution converges weakly to probability measure
$f_{V}$
on 
with bounded
$\kappa $
th moment; additionally, assume that
$\lim _{p \to \infty } \mathbb E_{\widehat f_{v}}(V^\kappa ) = \mathbb E_{f_{V}}(V^\kappa )$
. Then for any pseudo Lipschitz function 
of order
$\kappa $
,
$$ \begin{align*} \lim_{N\to \infty} \frac{1}{N} \sum_{j = 1}^{N} \psi(v_{j}) \stackrel {\mathrm{a.s.}}{=} \mathbb {E}[\psi(V)]. \end{align*} $$
Lemma 7. For any fixed unit vectors 
, 
and random matrices
$\tilde {{\boldsymbol X}}$
distributed as 
, it holds that
$\tilde {{\boldsymbol X}} u = \Big (\sum _{j = 1}^p \tilde {{\boldsymbol X}}_{1, j} u_j, \dots , \sum _{j = 1}^p \tilde {{\boldsymbol X}}_{n, j} u_j \Big )^\top \stackrel {d}{=} \frac {1}{\sqrt {n}}( Z_{1}, \ldots , Z_{n})^\top \sim N(0, \frac {1}{n}I_n)$
, where
$Z_{i}$
’s are i.i.d. samples of a standard normal
$N(0,1)$
distributed random variable Z. Similarly,
$\tilde {{\boldsymbol X}}^\top v = \Big (\sum _{i = 1}^n \tilde {{\boldsymbol X}}_{i, 1} u_i, \dots , \sum _{i = 1}^n \tilde {{\boldsymbol X}}_{i, p} u_i \Big )^\top \stackrel {d}{=} \frac {1}{\sqrt {n}}( Z^{\prime }_{1}, \ldots , Z^{\prime }_{p})^\top \sim N(0, \frac {1}{n}I_p)$
, where
$Z^{\prime }_{j}$
’s are i.i.d. samples of a standard normal
$N(0,1)$
distributed random variable
$Z'$
.
Lemma 8 (Shi et al. (Reference Shi, Shen, Hill and Van Den Hengel2012, Thm. 5))
Assume any 
and any random Gaussian matrix
$\tilde {{\boldsymbol X}}$
distributed as 
. For any
$\epsilon \in (0,1)$
, if the inner product
$\langle m_1, m_2 \rangle> 0$
, then with probability at least
$1 - 6 \exp (-\frac {n}{2}(\frac {\epsilon ^2}{2} - \frac {\epsilon ^3}{3})) $
, the following holds:
COMPETING INTEREST STATEMENT
The authors declare that no competing interests exist.
FUNDING STATEMENT
H.Z. is supported in part by NIH 1R01GM163244-01.
SUPPLEMENTARY MATERIAL
Zhou, J., & Zou, H. (2025). Supplement to “High-dimensional Newey–Powell test via approximate message passing,” Econometric Theory Supplementary Material. To view, please visit: https://doi.org/10.1017/S0266466626100450.
Open access
