Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,470 characters · 7 sections · 92 citation commands
Testing Heteroskedasticity Under Measurement Error
\doublespacing
As white1980heteroskedasticity famously observed, \textquotedblleft the presence of heteroskedasticity in the disturbances of an otherwise properly specified linear model leads to consistent but inefficient parameter estimates and inconsistent covariance matrix estimates \textquotedblright, testing whether the error terms in regression models satisfy the homoskedasticity assumption is of fundamental importance, in order to avoid invalid hypothesis testing and a substantial deterioration in the predictive accuracy of the model.
An extensive literature has emerged on testing for heteroskedasticity across a wide range of regression frameworks over the past half-century. Work includes bickel1978using, breusch1979simple, white1980heteroskedasticity, cook1983diagnostics, dette1998testing, zhu2001heteroscedasticity, dette2002consistent, zheng2009testing, su2013nonparametric, guo2020pairwise, tan2021testing and xu2021distance. Despite the maturity of this line of research, almost all existing procedures rest on a strong maintained assumption that all random variables are observed without error. This no–measurement-error assumption is rarely credible in empirical work, and the literature’s longstanding silence on it is largely driven by the severe technical difficulties.
A growing literature has examined inference and testing problems in regression models with measurement error. Early work includes carroll2006measurement and hu2012estimation, which primarily focus on identification and estimation. Building on these developments, more recent work further extends these ideas: dong2021average study the density-weighted average derivative estimator in measurement error models; otsu2021specification develop smoothing-based specification tests, and dong2022estimation develop semiparametric estimation for varying coefficient models with mismeasured regressors; dong2022nonparametric propose nonparametric significance tests based on deconvolution estimators. More fundamentally, for estimating characteristic functions under unknown measurement error distributions, delaigle2008deconvolution and kurisu2022uniform provide nonparametric identification and uniform convergence results, respectively. However, in this setting, heteroskedasticity testing becomes more complex, as the researcher must consistently estimate both the underlying regression function and the disturbance variance structure, potentially without knowledge of the measurement error distribution.
Unfortunately, if researchers ignore the noise in the data and mechanically apply conventional procedures to test for heteroskedasticity, they can easily be led to misleading conclusions. From a theoretical perspective, the seminal analysis of heteroskedasticity tests in linear errors-in-variables models in wooldridge1996solutions shows that, once the regressors are measured with error, the asymptotic behavior of the usual statistics based on auxiliary regressions of squared residuals no longer follows from the classical framework and has to be re-derived under additional regularity conditions. On the applied side, alica2025comparison, by comparing the simulation performance of traditional tests (such as those discussed in goldfeld1965some, park1966estimation, glejser1969new, harvey1976estimating, breusch1979simple and white1980heteroskedasticity) when the data are contaminated by measurement error, demonstrates that ignoring measurement error can cause the empirical size of these tests to deviate substantially from the nominal level and their power to deteriorate dramatically.
We next review the limited literature addressing tests for heteroskedasticity in the presence of measurement error. In an early contribution, carroll1992diagnostics proposed corrected residual plots that adjust for measurement error, together with associated formal tests. wooldridge1996solutions provided an asymptotic analysis of familiar Lagrange multiplier–type tests for homoskedasticity when the regressors are subject to classical measurement error. wallentin2002test showed, by simulation, that least-squares residuals may still be informative for detecting heteroskedasticity even when parameters are estimated by instrumental variables, whereas naive implementations of White’s test based on instrumental-variables residuals and observed regressors can exhibit severely inflated rejection rates. Most recently, romeo2024detecting developed a measurement-error–adapted version of White’s test, implemented via a model-based bootstrap. However, these contributions are all derived in simple linear regression settings that are difficult to satisfy in many empirical applications involving multiple covariates and more flexible parametric specifications.
The primary focus of our analysis is to extend heteroskedasticity tests for measurement error to parametric regression models that are more representative of empirical practice. Specifically, we develop integrated conditional moment (ICM)–type tests based on a residual-marked empirical process in the spirit of bierens1982consistent and bierens1990consistent. The presence of measurement error in the regressors is handled by employing deconvolution kernel estimators, as motivated by carroll1988optimal and stefanski1990deconvolving, to construct residuals. Critical values are then obtained through a computationally effective multiplier bootstrap. In implementing this bootstrap, we project the empirical process orthogonally onto the tangent space of the nuisance parameter to remove the \textquotedblleft parametric estimation effect\textquotedblright, which is initially discussed in durbin1973distribution and has long been regarded as an unavoidable and notoriously difficult obstacle to the practical implementation of ICM-type tests, thereby ensuring the validity of the bootstrap approximation.
It is important to highlight our contribution to addressing the long-standing difficulty of jointly recovering the regression relationship and the disturbance variance structure in the presence of measurement error. In particular, the resulting test statistics enjoy a parametric convergence rate, the associated bootstrap implementation is fast and conceptually transparent, and the procedure is robust with respect to tuning parameters, making it a promising tool for other testing problems beyond the specific setting considered here. In addition, we develop a companion testing procedure for the empirically realistic case in which the distribution of the measurement error is unknown, and we provide both asymptotic theory and empirical evidence for its performance.
The rest of the paper is organized as follows. In Section (ref), we outline the testing framework and the construction of our statistics. Then the asymptotic properties of the test statistics with some reasonable assumption are discussed under the null, local alternatives, and global alternative in Section (ref). The cases of unknown measurement error distribution are addressed in Section (ref) and the analysis about asymptotic behaviors of the proposed test statistics in Section (ref) parallels that in Section (ref) and differ only in the strength of assumptions. In Section (ref), we detail the implementation of the projection-based multiplier bootstrap procedure, emphasizing that the customized, explicitly constructed projection is particularly appealing due to its ease of interpretation and computational convenience. Results of Monte Carlo simulations and empirical studies are presented in Section (ref) and (ref), respectively. Additional simulation results and proofs are provided in the online supplementary appendix.
Let $Y$ be a response variable and $X$ be the scalar unobservable explanatory variable. Suppose that we consider a potentially heteroskedastic parametric regression model
where $g(\cdot;\theta_0)$ is the conditional mean function of a known parametric form characterized by the unknown parameter $\theta_0\in\Theta\subset\mathbb{R}^p$ and $U$ is the error term, which is the part of $Y$ not explained by $X$. In addition, the condition $\mathbb{E}[U\vert X]=0$ almost surely (a.s.) is imposed to ensure the correct specification of the conditional mean function $g(X;\theta_0)$. Moreover, since the true regressor $X$ is not directly observed, we instead observe a scalar variable $W$ through an additive measurement error model
where $\epsilon$ is assumed to follow the classical measurement error assumption in the sense that $\epsilon$ is independent of $X$. The density of $\epsilon$, $f_\epsilon$, is assumed to be known here and in Section (ref), while the case of an unknown density $f_\epsilon$ is addressed in Section (ref). We are interested in testing whether the conditional variance function of $U$ given $X$ is equal to an unknown positive constant. That is, our null hypothesis of interest is
and the alternative hypothesis is
which is the negation of $H_0$. Note that $\sigma^2_0$ is also the unconditional variance of $U$ when $H_0$ holds.
Following bierens1982consistent, bierens1990consistent, and bierens1997asymptotic, the conditional moment restriction in (ref) that characterizes $H_0$ can be equivalently expressed as a continuum number of unconditional moment restrictions as follows:
where ${\rm i}=\sqrt{-1}$ denotes the imaginary unit, $\Pi$ is a properly chosen compact set with nonempty interior, and the weighting function ${\rm e}^{{\rm i} X\xi}$ is selected from a parametric family indexed by $\xi$. In the absence of measurement error, the parameters $\theta_0$ and $\sigma^2_0$ are estimable via parametric methods, enabling the construction of test statistics based on a sample analog of equation (ref). However, in the presence of measurement error, the explanatory variable $X$ is unobservable and only the contaminated version $W$ is available. Motivated by dong2022nonparametric, equation (ref) can be reformulated in terms of the joint density of $(Y,X)$, which is denoted by $f_{Y,X}(y,x)$, as detailed below
Given a random sample $\{(Y_i,W_i)'\}_{i=1}^n$ of size $n\geq 1$ and motivated by carroll1988optimal and stefanski1990deconvolving, the unknown density $f_{Y,X}(y,x)$ can then be estimated using deconvolution kernels associated with the measurement error $\epsilon$,
where $K$ is a kernel, $b$ is bandwidth shrinking to zero at suitable rates and
For notational simplicity, throughout the paper, let $\mathcal{K}_\epsilon(x) = b\mathcal{K}_b(x)$ and let $f_\eta^{\text{ft}}(t)=\int {\rm e}^{{\rm i} tx}f_\eta(x)\,dx$ denote the Fourier transform of the generic random variable $\eta$. Let $\theta_n$ denote the estimator for $\theta_0$ based on the sample $\{(Y_i,W_i)'\}_{i=1}^n$ when the measurement error density $f_\epsilon$ is assumed to be known. The estimator for the variance $\sigma_0^2$ can be obtained by expressing $\mathbb{E}[(Y-g(X;\theta_0))^2]$ as an integral by replacing the joint density $f_{Y,X}(y,x)$ with its estimator $\hat{f}_{Y,X}(y,x)$ in (ref), which, after rearrangement, is given by
Based on the above estimators, we can express the sample version of $S(\xi,\theta_0,\sigma^2_0)$ in (ref) as the following empirical process:
which is expected to be close to zero under the null and deviate from zero under the alternative.
The test statistics are constructed based on appropriate distances from $S_n(\cdot,\theta_n,\sigma^2_n)$ to zero. We adopt the commonly used Kolmogorov--Smirnov (KS)-type statistic, which is based on the sup norm, and the Cram\'{e}r--von Mises (CvM)-type statistic, which relies on the squared norm, respectively:
where the uniform integrating measure on $\Pi$ is employed as recommended in dong2022nonparametric. The null hypothesis $H_0$ is rejected when the test statistics $KS_{n}$ and $CvM_{n}$ exceed their critical values, which are obtained via a computationally attractive multiplier bootstrap procedure detailed in Section (ref).
In this section, we study the asymptotic properties of the $KS_n$ and $CvM_n$ statistics introduced in Section (ref), under the assumption that the distribution of the measurement error is known. The case with unknown measurement error distribution will be discussed in detail in Section (ref). Based on the smoothness of the measurement error $\epsilon$, we distinguish between the ordinary smooth and supersmooth cases. For each case, we provide the required assumptions and establish the limiting distribution of the empirical process $S_n(\cdot,\theta_n,\sigma^2_n)$ in (ref) under the null. A sequence of local alternatives converging to the null at the parametric rate $n^{-1/2}$ is subsequently investigated, and the corresponding asymptotic local power is derived. Finally, the global power of the proposed test is discussed by analyzing the asymptotic behavior of $S_n(\cdot,\theta_n,\sigma^2_n)$ under the alternative.
In the following, let $f_X(x)$ denote the density function of $X$, $g^{(p)}(x;\theta)$ denote $p$-th derivative of function $g(x;\theta)$ with respect to $x$ and $\binom{a}{b}=a!/[(a-b)!b!]$ denote the binomial coefficient for nonnegative integers $a\geq b$. To facilitate the theoretical analysis, we begin by imposing a set of regularity conditions that hold for both the ordinary smooth and supersmooth cases.
Assumption (ref)({\romannumeral1}) imposes random sampling and the existence of the fourth moment of $Y$, which are necessary to ensure that the statistic constructed from $Y^2$ has a finite second-order moment. Assumption (ref)({\romannumeral2}) relaxes the convergence rate condition on the parameter compared with the case without measurement error. As pointed out in the literature, estimation in parametric models involving measurement error often fails to attain the parametric rate unless additional assumptions are imposed, see taupin1998estimation. To account for this, we relax the convergence rate assumption to allow rates between $n^{-1/4}$ and $n^{-1/2}$. Finally, Assumption (ref)({\romannumeral3}) adopts the classical measurement error assumption, which ensures the asymptotic zero mean of the test statistic and the existence of its variance. Following delaigle2008deconvolution, which highlights that ordinary smooth and supersmooth cases can be distinguished by the decay rate of the characteristic function of the measurement error, we first consider the setting where the characteristic function decays polynomially, and state the corresponding assumptions.
To facilitate a detailed analysis of the structural functions at the observed sample points $W_i$ using Taylor expansion, Assumption (ref)({\romannumeral1}) places smoothness restrictions on the structural functions. This follows the assumptions of the ordinary smooth case in existing work on measurement error (e.g., dong2022nonparametric), which impose Lipschitz continuity and integrability of the corresponding Lipschitz coefficients. Moreover, to ensure that the variance of the empirical process $S_n(\cdot,\theta_n,\sigma^2_n)$ in (ref) is not distorted by \textquotedblleft parameter estimation uncertainty\textquotedblright, initially discussed in durbin1973distribution, we strengthen the smoothness requirements to hold uniformly over a neighborhood of the true value $\theta_0$. Assumption (ref)({\romannumeral1}) strengthens the conventional ordinary smooth assumption by characterizing the exact limiting behavior of $f^{\mathrm{ft}}_\epsilon(\cdot)$, which is generalized to the form $f^{\mathrm{ft}}_\epsilon(t) = \exp({\rm i} t\zeta)/(c_0^{os}+c_{1}^{os}t+\cdots+c_{\alpha}^{os}t^{\alpha})$ for some real number $\zeta$ in fan1995average. This refinement is essential for deriving the precise asymptotic form of the test statistic; otherwise, stronger assumptions would be required to impose on a more complex expression involving the Fourier transform to ensure the existence of the variance and the asymptotic negligibility of the expectation. In addition, commonly used distributions such as the Laplace and Gamma distributions are included under this assumption. Accordingly, we introduce higher-order kernel functions, following the construction method provided in alexander2009deconvolution for kernels of arbitrary order in Assumption (ref)({\romannumeral3}). Kernel functions are used to characterize the integration involving the deconvolution kernel, thereby deriving the precise asymptotic form of the test statistic. In addition, higher-order kernels are also employed to establish the asymptotic negligibility of the expectation. Assumption (ref)({\romannumeral4}) requires that the bandwidth converges to zero sufficiently fast as the sample size increases. This so-called undersmoothing assumption, commonly used in the literature, is essential to guarantee the asymptotic negligibility of the expectation of the proposed statistics. Finally, Assumption (ref)({\romannumeral5}) is a necessary condition to ensure the boundedness of the variance of $S_n(\cdot,\theta_n,\sigma^2_n)$, as noted in dong2022nonparametric, and has been widely adopted in the literature on measurement error models.
Based on the assumptions stated above, the following theorem characterizes the asymptotic behavior of the empirical process $S_n(\cdot,\theta_n,\sigma^2_n)$ in (ref) under the null, in the case where the measurement error distribution is known and is of the ordinary smooth type. Let \textquotedblleft $\Longrightarrow$\textquotedblright denote weak convergence on $(l^{\infty}(\Pi),\mathcal{B}_\infty)$ in the sense of Hoffmann--J\orgensen, where $\mathcal{B}_\infty$ denotes the corresponding Borel $\sigma$-algebra, see, e.g., Definition $1.3.3$ in van1996weak.
Theorem (ref) shows that the proposed empirical process $S_n(\cdot,\theta_n,\sigma^2_n)$ converges at the parametric rate (i.e., $\sqrt{n}$) to a centered Gaussian process. As a consequence, the tests $KS_n$ and $CvM_n$ based on $S_n(\cdot,\theta_n,\sigma^2_n)$ also achieve the parametric rate.\footnote{Given the conditions in Theorem (ref), the continuous mapping theorem (see van1996weak) implies that the proposed $KS_n$ and $CvM_n$ statistics converge to the sup norm and the squared $L_2$-norm of the limiting Gaussian process $S_\infty^{os}(\cdot,\theta_0,\sigma^2_0)$, respectively. As similar techniques are used to derive the limiting distributions of these statistics under the null and the alternative hypotheses, for both ordinary smooth and supersmooth measurement errors, and in the multiplier bootstrap procedure, we do not repeat the details in the subsequent analysis.} Notably, this parametric rate is not affected by the use of deconvolution kernels, which are nonparametric estimators, highlighting a key advantage of the proposed procedure. In fact, such an advantage has also been observed in various methods that address measurement error problems using nonparametric estimators, including semiparametric estimators for regression functions, nonlinear regression estimation, and significance testing; see, for example, the discussions in hall2007semiparametric and dong2022nonparametric. It is also worth noting that although we use a deconvolution kernel to estimate $f_{Y,X}(y,x)$, the proposed test statistics evaluate the discrepancy over the entire parameter space for $\xi$. As a result, both the size accuracy and power of the tests exhibit robustness with respect to the bandwidth choice, which will be further discussed in Section (ref).
Next, we consider the scenario where the characteristic function of the measurement error decays at an exponential rate, a setting commonly referred to as the supersmooth case, as described in Assumption (ref)({\romannumeral2}). The required assumptions are presented below.
It is worth noting that for the supersmooth case, Assumptions (ref)({\romannumeral1}), ({\romannumeral2}), and ({\romannumeral3}) impose stronger conditions than their counterparts for the ordinary smooth case in Assumption (ref). Specifically, we require that the characteristic function of the measurement error be of exponential type, the structural functions be infinitely differentiable, and the kernel functions be of infinite order. While these assumptions are more restrictive, they are satisfied by many practical examples—for instance, using polynomials, exponentials, or circular functions as structural functions (with additional examples provided in alexander2009deconvolution); normally distributed measurement error; and the infinite-order kernel constructed in mcmurry2004nonparametric, which is also employed in our simulation study in Section (ref). Moreover, such strengthened assumptions are necessary in the supersmooth setting. Without these strengthened assumptions, two main difficulties arise. First, it becomes challenging to ensure the boundedness of the variance of the test statistic, as required in Assumption (ref)({\romannumeral5}). Second, one would need to impose additional conditions to guarantee the asymptotic negligibility of the bias, such as the undersmoothing bandwidths and the Lipschitz continuity of the structural functions stated in Assumption (ref), which would significantly complicate the set of required assumptions. In addition, the conditions can be relaxed to permit the measurement error to be distributed as the convolution of a Gaussian density and any ordinary smooth density that satisfies Assumption (ref)({\romannumeral2}), thereby broadening the applicability of the supersmooth case.
Theorem (ref) shows that our empirical process $S_n(\cdot,\theta_n,\sigma^2_n)$ also achieves the parametric rate in the supersmooth case, and that the limiting Gaussian process exhibits a variance structure similar to that of the ordinary smooth case. However, existing literature suggests that, in other measurement error-related problems, the convergence rate and the variance structure of the limiting null process may differ substantially between the ordinary smooth and the supersmooth cases. In particular, otsu2021specification demonstrates that, in specification testing based on nonparametric estimators, the supersmooth case often yields a slower convergence rate and a more complex asymptotic structure. A detailed comparison between the ordinary smooth and the supersmooth cases is beyond the scope of this paper. Instead, in Section (ref), we focus on constructing alternative test statistics when the distribution of the measurement error is unknown.
We next introduce a sequence of local alternatives and evaluate the asymptotic behavior of the proposed tests under these alternatives, demonstrating that the tests possess nontrivial local power to detect nonconstant variance. Specifically, we assume
where $\Delta:\mathbb{R}\to\mathbb{R}$ is a nonzero function satisfying $\mathbb E[\Delta(X)]=0$ and $\mathbb{E}|\Delta(X)|<\infty$, so that the sequence of local alternatives converges to the null hypothesis at the parametric rate. We then investigate the asymptotic convergence results of the test statistic separately under the ordinary smooth and the supersmooth cases.
Theorem (ref) implies that under the sequence of local alternatives $H_{1n}$ in (ref), the empirical process $S_{n}(\cdot,\theta_n,\sigma^2_n)$ still converges at the parametric rate to its limiting distribution. Notably, for both the ordinary smooth and the supersmooth cases, the limiting distribution under $H_{1n}$ differs from that under $H_0$ by a deterministic shift function $\mu_\Delta(\xi)$. This shift corresponds to the covariance between $\Delta(X)$ and ${\rm e}^{{\rm i} X\xi}$, ensuring that the limiting distributions under $H_{1n}$ and $H_0$ are different and thus yielding nontrivial asymptotic local power for the proposed $S_{n}(\cdot,\theta_n,\sigma^2_n)$.
At the end of this section, we establish the global consistency of $S_{n}(\cdot,\theta_n,\sigma^2_n)$ by showing that it possesses nontrivial power under the alternative hypothesis $H_1$ in (ref), as formalized in the following theorem. Let $\sigma_\ast^2$ denote the probability limit of $\sigma_n^2$ under $H_1$, i.e., the pseudo-true value and the unconditional variance of $U$. Note that $\sigma_\ast^2=\sigma_0^2$ under $H_0$.
Under $H_1$ where the conditional variance $\mathbb E[U^2|X]$ is a nonconstant function of the unobserved regressor $X$, Theorem (ref) implies that $S_{n}(\xi,\theta_n,\sigma^2_n)$ converges in probability uniformly to the nonzero constant function $S(\xi,\theta_0,\sigma^2_\ast)$. Consequently, the test statistics $KS_n$ and $CvM_n$ based on $\sqrt nS_{n}(\xi,\theta_n,\sigma^2_n)$ diverge to positive infinity in probability, thereby guaranteeing the consistency of $KS_n$ and $CvM_n$ against $H_1$.
In this section, we focus on constructing tests for heteroskedasticity when information about the measurement error distribution is unavailable. The necessity arises because a misspecified measurement error model can lead to biased estimators, as noted in hall2007semiparametric. The assumption of replicated measurements has been widely adopted in numerous theoretical studies, see, e.g., carroll1992diagnostics and delaigle2008deconvolution. In practice, repeated measurements can be obtained in certain settings, for example, through high-frequency observations over a short time interval or by rapidly generating data from simulated systems. As noted in allen1999checking, climate scientists, when applying optimal fingerprinting methods to detect and attribute anthropogenic influences on climate, can readily obtain repeated measurements from a climate model's control run. As commonly stated in the literature carroll1992diagnostics, repeated measurements on unobservable regressors are usually required in multiple replicates to obtain accurate information. One of the notable advantages of our method is that a single set of repeated measurements suffices to guarantee favorable asymptotic behavior of the test statistics, as detailed in the setting below:
where $W$ and $W^r$ denote a pair of repeated measurements on the latent variable $X$ and $(\epsilon,\epsilon^r)$ are assumed to be i.i.d. measurement errors. Additionally, we make the following assumptions,
Based on the repeated measurements described in Assumption (ref)({\romannumeral1}), Assumption (ref)({\romannumeral2}) first imposes a symmetry condition on the distribution of the measurement error, which is equivalent to requiring that its characteristic function be real-valued. It then places relatively strong moment conditions on the measurement error $\epsilon$. Commonly used measurement error distributions, such as the Laplace and normal distributions, satisfy these conditions. While some studies have considered relaxations of this assumption, for instance, li1998nonparametric relaxes the symmetry condition, and delaigle2016methodology weakens the repeated measurement requirement by imposing stronger restrictions on the true regressor, we do not pursue these extensions here for simplicity. Assumption (ref)({\romannumeral2}), motivated by the later analysis of the asymptotic behavior of the derivative of the unknown characteristic function of the measurement error, is discussed in detail in kurisu2022uniform. It characterizes the integration of the estimated deconvolution kernel in our work.
Given the availability of replicate measurements for each sample $X_i$, we can obtain an estimator for $\theta_0$, which we denote as $\hat\theta_n$. Note that $\hat\theta_n$ is different from $\theta_n$, which is constructed under the known measurement error distribution assumption. The repeated measurements are also used to estimate the characteristic function of the measurement error, which is then plugged into the previously constructed deconvolution kernel to obtain an estimator, as given by
The estimated kernel $\hat{\mathcal{K}}_b(x)$ is subsequently used in place of the original kernel in the proposed estimator of the variance (ref) and the empirical process (ref) to obtain
and
The corresponding test statistics with repeated measurements, denoted by $\widehat{KS}_n$ and $\widehat{CvM}_n$, are then constructed in the same manner as the original $KS_n$ and $CvM_n$ statistics, respectively. The bootstrap procedure for obtaining critical values will be discussed in detail in Section (ref). In addition, we introduce the notation $\Pi_\epsilon$ to evaluate the accuracy of $\hat{f}^{\mathrm{ft}}_\epsilon(\cdot)$ as an estimator of $f^{\mathrm{ft}}_\epsilon(\cdot)$,
In fact, due to the repeated measurements and the symmetry of the error distribution as assumed in Assumption (ref), $\Pi_\epsilon$ is unbiased.
We still distinguish between the ordinary smooth and the supersmooth cases based on the decay rate of the characteristic function of the measurement error in the case of an unknown measurement error distribution due to the differing conditions required to derive the asymptotic properties of the reconstructed test statistics. In the first case, the ordinary smooth scenario, we add the following assumption to Assumption (ref).
Assumption (ref) is essential, as the estimation of $\hat{f}^{\mathrm{ft}}_\epsilon(\cdot)$ introduces additional higher-order terms and, unfortunately, distorts the variance structure of the main term. By imposing Assumption (ref)({\romannumeral1}), we establish a lower bound on the bandwidth, which ensures the negligibility of the unavoidable higher-order terms and provides a theoretically justified bandwidth range combined with the upper bound specified in Assumption (ref)({\romannumeral4}). In addition, estimating the error characteristic functions modifies the form of the main term of the original empirical process. Specifically, it adds an additional component to the leading term $r^{os}_{\infty}(Y,W,W^r;\cdot,\theta_0,\sigma^2_0)$ that involves the stochastic term $r^{\epsilon,os}_{\infty}(Y,W,W^r;\cdot,\theta_0,\sigma^2_0)$ mentioned in Assumption (ref)({\romannumeral2}), thereby necessitating a strengthening of Assumption (ref)({\romannumeral5}) through the requirement of bounded second moment of $r^{\epsilon,os}_{\infty}(Y,W,W^r;\cdot,\theta_0,\sigma^2_0)$ to ensure the finiteness of the variance. We then establish the asymptotic theory for the reconstructed test statistics for the ordinary smooth case with unknown measurement error distribution, as stated in the following theorem.
Theorem (ref) establishes that the reconstructed empirical process still converges at the parametric rate\footnote{The asymptotic null distributions of statistics $\widehat{KS}_n$ and $\widehat{CvM}_n$ are established in the same manner as in Theorem (ref). By applying the continuous mapping theorem, it follows directly that $\widehat{KS}_n$ and $\widehat{CvM}_n$ converge to the sup norm and the squared $L_2$-norm of the limiting Gaussian process $\hat{S}^{os}_\infty(\cdot,\theta_0,\sigma^2_0)$, respectively.}, unaffected by the error brought by the estimation of the error characteristic functions. In contrast to the case with a known measurement error distribution, the covariance structure of the limiting process differs and requires more restrictive assumptions on the variance.
For the supersmooth case, we assume
Analogously, for the supersmooth case with unknown measurement error, the estimation brings a higher-order term, rendered asymptotically negligible through a stronger bandwidth condition in Assumption (ref)({\romannumeral1}), as well as additional uncertainty in the main term, which necessitates a more stringent condition to ensure variance boundedness. Compared with the ordinary smooth case under an unknown measurement error distribution, Assumption (ref)({\romannumeral2}) is imposed to ensure the finiteness of the variance. A similar comparison can be observed in Assumptions (ref) and (ref) when the error distribution is known.
In general, the lack of measurement error information and the supersmooth case are typically considered to lead to slower convergence rates. Nevertheless, the reconstructed empirical process proposed in this paper retains the parametric rate of convergence, a desirable property not well established in the testing literature when measurement error is present. Overall, the findings from the two cases suggest that estimating the error characteristic functions does not alter the requirements on the smoothness of the structural functions or on the order of the kernel functions. Moreover, it does not affect the convergence rate of the test statistic toward its asymptotic null distribution. However, it leads to different admissible range for the bandwidth and changes the asymptotically equivalent form of $\hat{S}_{n}(\cdot,\hat\theta_n,\hat{\sigma}^2_n)$, denoted respectively by $\hat{r}^{os}_{\infty}(Y,W,W^r;\cdot,\theta_0,\sigma^2_0)$ and $\hat{r}^{ss}_{\infty}(Y,W,W^r;\cdot,\theta_0,\sigma^2_0)$ for ordinary smooth and supersmooth, resulting in different variance structures of the limiting process.
Similar to the case where the measurement error distribution is known, we claim that the proposed test exhibits nontrivial local power even in the absence of the information about the measurement error distribution, as demonstrated by studying the asymptotic properties of the reconstructed statistics under the sequence of local alternative hypotheses described in (ref). We present the following theorem to establish this result.
We observe that when the alternative hypothesis deviates from the null at $\sqrt{n}$-rate, the limiting distribution of the test statistic, whether for the ordinary smooth or supersmooth case, adds a deterministic shift term to the corresponding null distribution, aligning with the case where the measurement error distribution is known. Furthermore, we proceed to analyze the limiting distribution of the test statistics under the alternative hypothesis (ref) to investigate the global power of our test.
Theorem (ref) shows that $\sqrt{n}\hat{S}_{n}(\xi,\hat{\theta}_n,\hat{\sigma}^2_n)$ diverges to infinity as $n\to\infty$ under $H_1$, thereby confirming that the consistency of the tests is not affected by the absence of information about the measurement error distribution. To summarize, when the measurement error distribution is unknown, the convergence rates of the test statistics under the null, as well as their local and global powers, remain unchanged, thereby facilitating the bootstrap procedure developed in Section (ref).
Theorems in Sections (ref) and (ref) demonstrate that the limiting processes in both cases are case-dependent; specifically, the variance functions depend on the underlying data distribution, which motivates the use of a bootstrap method. In this section, we derive the corresponding test statistics for implementing a computationally attractive multiplier bootstrap and obtain the critical values.
In classical bootstrap methods (see, for example, koul1994bootstrapping), researchers typically compute residuals from observed data and estimate their joint population distribution. Bootstrap samples are then generated by resampling from the estimated distribution, and the bootstrap test statistics are constructed in the same manner as the original statistics. This so-called residual-based bootstrap approach, however, faces substantial challenges in the presence of measurement error. As noted in dong2022nonparametric, the unavailability of the true regressors prevents direct estimation of the residual distribution. Consequently, researchers have resorted to deconvolution kernel techniques, as suggested in hall2007testing. Nonetheless, such methods suffer from slow convergence rates, difficulties in selecting tuning parameters, and high computational costs. These limitations motivate the development of alternative bootstrap procedures. Inspired by van1996weak, we compute the corresponding term from each observation and multiply each of these terms by an independent mean-zero, unit-variance random variable before summing them to obtain the bootstrap version of the statistics. By performing $B$ rounds of multiplier resampling, we obtain $B$ bootstrap versions of the test statistics.
Unfortunately, the application of the multiplier bootstrap is complicated by the \textquotedblleft parameter estimation effect\textquotedblright, as noted in durbin1973distribution. Specifically, the previously proposed empirical process for the presence of measurement error can be rewritten as
and a similar formulation for the case where the measurement error distribution is unknown. The first term, which corresponds to the empirical process constructed using the true parameter values, is asymptotically equivalent to the leading term in each case. As for the second term, we claim that it converges to zero at $\sqrt{n}$-rate in probability, under reasonable assumptions on the structural functions, kernels and bandwidth (as discussed in Sections (ref) and (ref)), and by employing the parameter estimation method satisfying the conditions in Assumption (ref)({\romannumeral2}). That is, although estimation uncertainty for $\theta_0$ exists, it does not affect the asymptotic behavior of the test statistics. The third term, commonly referred to as \textquotedblleft parameter estimation effect\textquotedblright, converges to the limiting process at rate $O_p(\sigma_n^2-\sigma_0^2)$ in the original statistic but at a slower rate during the bootstrap procedure. Henceforth, the notation $O_p(\sigma_n^2-\sigma_0^2)$ denotes the convergence rate in probability of the stochastic process to its limiting distribution, which is the same as the rate at which $\sigma^2_n$ converges in probability to $\sigma_0^2$. We address this issue by modifying the weight function ${\rm e}^{{\rm i} X\xi}$ in (ref) to its mean-centered version in the construction of the bootstrap version of statistics, that is, subtracting the sample analog of its expectation,
where
and $\{V_i\}_{i=1}^n$ is a sequence of i.i.d. random variables with mean zero, variance one, bounded support, and independent of the sample $\{(Y_i,W_i)'\}_{i=1}^n$, e.g., mammen1993bootstrap's two-point distribution: $$ V_i=\left\{
\right. $$
The above construction ensures that $(\sigma_n^2-\sigma_0^2)\mathbb{E}({\rm e}^{{\rm i} X\xi}-\mathbb{E}{\rm e}^{{\rm i} X\xi}) = 0$, which in turn asymptotically eliminates the influence of the \textquotedblleft parameter estimation effect\textquotedblright arising from $\sigma_n^2$. As a result, the bootstrap versions of empirical processes $S_n^{\ast}(\cdot,\theta_n,\sigma^2_n)$ converge to the same limiting processes as the original counterparts $S_n(\cdot,\theta_n,\sigma^2_n)$ for both cases under the null hypothesis, as shown in Theorem (ref). In the bootstrap procedure, the test statistics $KS^{\ast}_n$ and $CvM_n^{\ast}$ can be constructed in a completely analogous way to those in Section (ref), based on the proposed $S_{n}^{\ast}(\cdot,\theta_n,\sigma^2_n)$. Specifically, we replace the original process $S_{n}(\cdot,\theta_n,\sigma^2_n)$ in $KS_n$ and $CvM_n$ with its bootstrap counterpart $S_{n}^{\ast}(\cdot,\theta_n,\sigma^2_n)$ and compute the sup norm and the squared $L_2$-norm accordingly.
Defining \textquotedblleft $\overset{\ast}{\Longrightarrow}$\textquotedblright as weak convergence and $\mathbb{P}_n^\ast$ as the bootstrap probability under the bootstrap law, i.e., conditional on the original sample, see, e.g., Section 2.9 of van1996weak, the validity of the proposed multiplier bootstrap is formally established by the following theorem.
The centering adjustment of the weighting function does not affect the parametric-rate convergence of the bootstrap version of the empirical process to the limiting process. Moreover, the bootstrap version $S_{n}^{\ast}(\cdot,\theta_n,\sigma^2_n)$ converges to the same limiting process as that of the original test process $S_{n}(\cdot,\theta_n,\sigma^2_n)$ under the null hypothesis, for both the ordinary smooth and supersmooth cases, as established in Section (ref). Notably, the bootstrap limiting distributions for $KS^{\ast}_n$ and $CvM_n^{\ast}$ remain unchanged under both the null and local alternative hypothesis\footnote{The convergence results stated here follow directly from the continuous mapping theorem.}. By contrast, the limiting distribution of the test statistic under local alternatives includes an additional deterministic shift, which shifts its value away from the null distribution, making it more likely to exceed the bootstrap critical value and thus leading to the rejection of the null hypothesis. This ensures the validity of the bootstrap approach with respect to both size control and nontrivial local power. As a consequence of the above analysis, the asymptotic critical value at the significance level $\alpha$ is $c^\ast_{\alpha}=\inf\{c_\alpha\in[0,\infty):\lim_{n\to\infty}\mathbb{P}_n^\ast\{KS_n^\ast>c_\alpha\}=\alpha\}$, where we take $KS_n^\ast$ as an example and we also note that the bootstrap procedure for $CvM_n$ can be implemented in an analogous way. In practice, $c^\ast_{\alpha}$ can be approximated as $c^\ast_{n,\alpha} = \{KS_n^\ast\}_{B(1-\alpha)}$, the $B(1-\alpha)$-th order statistic for $B$ replicates $\{KS_n^\ast\}_{b=1}^B$ and we reject $H_0$ if $KS_n>c^\ast_{n,\alpha}$.
For the case where the measurement error distribution is unknown, empirical processes constructed in the manner above are not valid, as they fail to account for the additional terms arising from the estimation of the unknown characteristic function. Specifically, in the ordinary smooth case, directly replacing the deconvolution kernel in $S_{n}^{\ast}(\cdot,\theta_n,\sigma^2_n)$ by $\hat{\mathcal{K}}_b(\cdot)$ causes the zero-mean multipliers to eliminate the estimation effect brought by $\hat{f}^{\mathrm{ft}}_\epsilon(\cdot)$. Consequently, the bootstrapped empirical process converges to $S_\infty^{os}(\cdot,\theta_0,\sigma_0^2)$ rather than to the expected null limit $\hat{S}_\infty^{os}(\cdot,\theta_0,\sigma_0^2)$. This discrepancy invalidates the bootstrap approximation, and a similar issue persists in the supersmooth case. Fortunately, motivated by dong2022nonparametric, we can perturb the estimator of the characteristic function using a set of multipliers with unit mean and unit variance. Specifically, we introduce a new sequence of multipliers $\{V^\ast_i\}_{i=1}^n$ satisfying $\mathbb{E}(V_i^\ast)=1$ and $\text{Var}(V_i^\ast)=1$ (e.g., the standard exponential distribution) and use them to construct perturbed analogs of the characteristic function estimator and the deconvolution kernel as in equation (ref),
The resulting kernel estimator $\hat{\mathcal{K}}^\ast_b(\cdot)$ is then substituted into equations (ref) and (ref) to obtain the bootstrap version of the variance estimator ${\hat{\sigma}_n^2}{}^\ast$ and the bootstrap empirical process, as defined by
and
The following theorem confirms the validity of the bootstrap procedure even in the absence of information about the measurement error distribution, as in the case with a known measurement error distribution.
Under the null or the local alternatives, the bootstrap versions of the empirical processes with a centering adjustment for the ordinary smooth and supersmooth cases converge to the same limiting processes $\hat{S}^{os}_\infty(\cdot,\theta_0,\sigma_0^2)$ and $\hat{S}^{ss}_\infty(\cdot,\theta_0,\sigma_0^2)$ as those of the original statistics constructed in Section (ref), respectively. A key implication of the above theorem is that the bootstrap versions of statistics can be constructed based on centralized correction of $\hat{S}^{\ast}_{n}(\cdot,\hat\theta_n,{\hat{\sigma}_n^2}{}^\ast)$,
For the same reasons discussed earlier, these bootstrap statistics converge, under both the null hypothesis and local alternatives, to the supremum norm and the squared $L_2$-norm of the corresponding limiting Gaussian processes $\hat{S}_{\infty}^{os}(\cdot,\theta_0,\sigma_0^2)$ for the ordinary smooth case and $\hat{S}_{\infty}^{ss}(\cdot,\theta_0,\sigma_0^2)$ for the supersmooth case, respectively. Combined with the fact established in Section (ref) that, under local alternatives, the limiting distributions of the statistics are perturbed by a deterministic shift function, this result ensures the validity of the proposed bootstrap procedure.
It is worth highlighting the projection-based approach proposed in sant2019specification, yang2024model, and song2025unified, which intuitively eliminates the \textquotedblleft parameter estimation effect\textquotedblright term by the orthogonal projection of the weighting function onto the tangent space of nuisance parameters in probability space. From this perspective, the centering method employed in this paper corresponds to projecting the weighting function onto the constant function $1$. It is also worth noting that, compared with the wild bootstrap methods discussed in other studies, our approach offers the significant advantage of greater computational efficiency. Moreover, as will be discussed in detail in Section (ref), the selection of tuning parameters such as the bandwidth does not pose practical difficulties when using the multiplier bootstrap.
In the simulations that follow, we aim to examine several key aspects. First, we evaluate whether the proposed tests, under both known and unknown measurement error distributions, maintain size close to the nominal level and exhibit reasonable power. Second, we investigate whether test performance differs substantially between the ordinary smooth and supersmooth cases. Third, we explore how different model specifications influence the validity of heteroskedasticity diagnostics. Finally, we assess whether the choice of bandwidth has a noticeable impact on the size and power characteristics of the tests.
We start with a linear model specification, hereafter referred to as model $1$. In this setting, data are generated according to $Y=\alpha_0+\alpha_1X+\sigma U$, where we set $\alpha_0=\alpha_1=1$, and assume that both $X$ and $U$ follow independent normal distributions with mean zero and unit variance. Heteroskedasticity is introduced through $\sigma$. To investigate the performance of the proposed tests, we consider three distinct data-generating processes (DGPs) that represent both the null and alternative hypotheses.
Specifically, DGP(0) corresponds to data generated under the null hypothesis, while DGP(1) and DGP(2) represent alternative scenarios with relatively high-frequency and low-frequency deviations from the null, respectively. In practical measurement error problems, the covariate $X$ is unobservable. We introduce additive noise $\epsilon\sim N(0,1/3)$, independent of $X$, such that the signal-to-noise ratio, defined as $Var(U)/Var(\epsilon)$, equals $3$. The observed contaminated variable is then given by $W=X + \epsilon$. We employ an infinite-order kernel defined by its Fourier transformation that is widely used in the measurement error literature and has been shown to be effective in a variety of applications (see, for example, mcmurry2004nonparametric, dong2021average, and dong2022nonparametric):
We estimate the parameters $\alpha_0$ and $\alpha_1$ using the method proposed in cheng1998polynomial. It is worth noting that under the null hypothesis, in which the model exhibits homoskedasticity, the estimators can be applied directly. Under the alternative hypothesis, the omission of heteroskedasticity in the estimation procedure does not compromise the performance of the test, since, as discussed in Section (ref), the estimation of $\alpha_0$ and $\alpha_1$ does not introduce any additional \textquotedblleft parameter estimation effect\textquotedblright on the asymptotic behavior of the test statistics. Subsequently, the bandwidth is selected based on the rule-of-thumb proposed in delaigle2008deconvolution and evaluated over a grid of candidate values, with the aim of minimizing pointwise mean squared error of the estimated characteristic function of $f^{\mathrm{ft}}_\epsilon(\cdot)$ in the unknown measurement error setting. Specifically, we set $b=c(5\sigma^4/n)^{1/27}$ for the ordinary smooth case and $b=c(4\sigma^2/\log(n))^{1/2}$ for the supersmooth case, where $c$ takes integer values from $1$ to $10$. To evaluate the robustness of the proposed test, we consider sample sizes $n\in\{250,500,1000\}$ and nominal significance levels $\alpha\in\{0.01,0.05,0.1\}$. Critical values are obtained using $199$ bootstrap replications, and $1000$ times Monte Carlo experiments. Due to space limitations, we present only the results for $\alpha=0.05$ and $c\in\{0.1,0.5,1\}$ with the remaining results reported in the Appendix (ref) on the online supplementary appendix.
Tables (ref) and (ref) report the simulation results under known and unknown measurement error distributions, respectively. The proposed tests exhibit desirable size accuracy and power properties across all considered cases. Specifically, as the sample size increases, the empirical size approaches the nominal level, and the power improves accordingly. A comparison between the results under DGP(1) and DGP(2) reveals that the test tends to be more powerful against low-frequency alternatives, which is consistent with results reported in the literature employing global test procedures. In our test, the convergence rate of the test statistic under the supersmooth case is $\sqrt{n}$, the same as that under the ordinary smooth case, which allows the supersmooth case to avoid the typical drawback of reduced power often reported in the literature. Interestingly, we observe that the test exhibits slightly higher power in the supersmooth case. We conjecture that this improvement stems from the test statistic and the associated plug-in estimators being exactly unbiased under the supersmooth setting, rather than merely asymptotically unbiased in the ordinary smooth case. To examine this conjecture, we further consider a constant model, where the data are generated from $Y=\alpha_0+\sigma U$, and conduct simulations under the same experimental settings as in the linear model.
As shown in Tables (ref) and (ref), the test exhibits comparable power across the ordinary smooth and supersmooth cases, which supports our conjecture that, in small samples and under the parametric model setting, the exact unbiasedness of the plug-in estimators in the supersmooth case leads to improved performance. Finally, different bandwidth choices across the grid points do not substantially affect the size and power of the test, confirming the desirable robustness to tuning-parameter selection as discussed in Sections (ref) and (ref).
Overall, the simulation results demonstrate that our proposed test is particularly sensitive to low-frequency alternatives, exhibits improved power for the supersmooth case, and is robust to bandwidth selection.
In this section, we illustrate the practical utility of our proposed tests by applying them to classic datasets from the measurement error literature. We first focus on detecting heteroskedasticity in the relationship between yields of corn and determinations of available soil nitrogen collected on Marshall soil in Iowa, a context where measurement error typically arising because only a small soil sample is taken from each plot and because of noise in the chemical analysis used to determine the nitrogen level in the soil, is often taken into account (see fuller2009measurement), whereas heteroskedasticity in the regression equation caused by differences in sampling times and conditions is often ignored. We use the dataset originally presented in fuller2009measurement (Tables $1.2.1$ and $3.1.1$). This dataset comprises observations on corn yield $Y$ and available soil nitrogen $X$ for a series of experimental plots. The primary econometric objective is to estimate the production relationship:
where $X_i$ is the true soil nitrogen content but is contaminated by measurement error, so that the researcher only observes the variable $W_i$, which is generated as $W_i = X_i + \epsilon_i$, with $\epsilon_i$ denoting the measurement error.
Standard analyses of this dataset (e.g., Ordinary Least Squares) typically ignore the measurement error, leading to attenuation bias in the slope coefficient. More sophisticated errors-in-variables (EIV) methods account for the bias but often maintain the strong assumption of homoskedasticity, $(Var(\epsilon_i\vert X_i)=\sigma_0^2)$. If the variance of the error term depends on the level of nitrogen (e.g., higher variability in yield at higher nitrogen levels), the standard errors of these EIV estimators may be invalid, leading to incorrect conclusions mentioned in Section (ref). Therefore, testing for heteroskedasticity in the presence of measurement error is a crucial diagnostic step. Specifically, we consider three scenarios corresponding to the data structures available in fuller2009measurement: known measurement error variance where we make use of the information $\sigma_{\epsilon_i}=57$ provided in Example $1.2.1$ of fuller2009measurement and specify the measurement error distribution to be Laplace and normal, representing ordinary smooth and supersmooth errors, respectively, and unknown measurement error distribution where we utilize the replicate measurements provided in Table $3.1.1$. We apply our proposed KS and CVM-type tests to this dataset. To evaluate the sensitivity of our procedure to tuning parameters and distributional assumptions, we report $p-$values across a range of bandwidth parameters $c$, and the results are summarized in Table (ref).
First, evidence of heteroskedasticity can be drawn from the empirical results presented in Table (ref). Specifically, if we adopt a standard significance level of $\alpha = 0.10$, the null hypothesis of homoskedasticity is rejected in most experimental settings. This suggests that the conditional variance of corn yield varies with soil nitrogen levels. Ignoring this heteroskedasticity, as is common in standard linear EIV regressions applied to this data, could lead to inefficient estimation and misleading confidence intervals. Our test provides a robust means of detecting latent structures otherwise obscured by measurement error. Second, we find robustness to distributional assumptions, bandwidth choice, and statistics selection. In practical applications, the precise distribution of the measurement error is rarely known. Fortunately, a comparison between the \textquotedblleft Laplace\textquotedblright and \textquotedblleft Normal\textquotedblright columns reveals that the resulting $p$-values are largely insensitive to this misspecification, although the theoretical derivation of the deconvolution kernel depends on the error distribution. Furthermore, consistent with the Monte Carlo simulations in Section (ref), the proposed tests exhibit robustness with respect to a wide range of bandwidth constants $c$. Both statistics yield comparable inferences, reinforcing the reliability of the testing procedure. Third, the results in the \textquotedblleft Unknown\textquotedblright columns demonstrate that our test remains powerful even when the error variance is not known and must be estimated via repeated measurements (as detailed in Section (ref)). However, it is worth noting that the admissible range for the bandwidth choice appears slightly narrower than in the known-variance case. This is an expected consequence of the additional estimation noise introduced by estimating the characteristic function of the error term, suggesting that practitioners should exercise caution in bandwidth selection when relying on replicates.
To further demonstrate the applicability of our test for detecting latent heteroskedasticity, we consider a second application that focuses on estimating Engel curves. This analysis utilizes data from the 2023 Consumer Expenditure Survey (CES), a dataset widely used in demand analysis (e.g., hausman1995nonlinear). We adopt the standard Leser--Working functional form (leser1963forms), which relates the budget share of a commodity to the logarithm of total expenditure. For household $i$ and commodity $j$ in quarter $t$, the structural equation is specified as:
where $Y_{ijt}$ represents the budget share of commodity $j$ (e.g., Food, Clothing), $X_{it}$ denotes the true log total expenditure and the error term $U_{ijt}$ captures heterogeneity in preferences. However, total expenditure is notoriously difficult to measure accurately in survey data. Following the literature (e.g., hausman1995nonlinear), we assume that the observed log total expenditure, $W_{it}$, measures the true latent log total expenditure $X_{it}$ with an additive error:
where $\epsilon_{it}$ represents the measurement error. Since the variance of the measurement error is unknown, we adopt the repeated measurement outlined in Section (ref). Specifically, we treat the log total expenditure reported in the subsequent quarter $x_{i, t+1}$ as a replicate measurement for the current quarter $t$. This allows us to estimate the characteristic function of the measurement error and construct the test statistic without making parametric assumptions about the error distribution. We analyze five major expenditure categories: Food, Clothing, Recreation, Health care, and Transportation. The tests are conducted separately for the second, third, and fourth quarters (Q2, Q3, Q4) of 2023 to check for temporal stability. The results of the $p$-values for the KS and CvM test statistics are presented in Table (ref).
The empirical results in Table (ref) reveal distinct patterns across commodity groups. Most notably, \textquotedblleft Transportation\textquotedblright exhibits strong evidence of heteroskedasticity. The $p$-values are consistently close to zero across all quarters for both KS and CvM statistics. This rejection of the null hypothesis aligns with the findings in hausman1995nonlinear, who argued that the simple linear-in-log Leser-Working specification is insufficient for certain goods, particularly Transportation. They found statistically significant coefficients for higher-order polynomial terms (quadratic in log expenditure). If the true relationship contains a quadratic term (e.g., $\alpha_2 X_{it}^2$) but is modeled linearly, the omitted nonlinear component becomes part of the error term, causing the error variance to vary with the level of expenditure. Therefore, our robust rejection of the homoskedasticity null hypothesis corroborates the need for higher-order Engel curve specifications as suggested by hausman1995nonlinear.
In this online supplementary appendix, we report additional simulation results in Section (ref), introduce necessary notations and definitions in Section (ref), and provide proofs of the main theoretical results in Section (ref). Finally, auxiliary lemmas and their proofs are collected in Section (ref).