Endogenous Quantile Regression with Measurement Error in Dependent Variable

Endogenous Quantile Regression with Measurement Error in Dependent Variable

\onehalfspacing

\microtypesetup{ activate=false, protrusion=false, expansion=false }

Introduction

Quantile regression (QR) has been extensively studied in theoretical econometrics koenker1978regression, koenker2005quantile and has gained substantial traction in empirical work. It helps capture the distributional effects of covariates on the outcome variable, not only at the center but also heterogeneously across different quantiles.

In observational data, however, two practical challenges frequently arise. First, variables of interest are often endogenous, rendering standard QR estimators inconsistent for marginal quantile effects. Second, variables used in empirical analyses are commonly subject to measurement error, especially in self-reported survey data, such as income bound1991extent, bound2001measurement, health status bound1999dynamic, or receipt of government transfers meyer2009under, among many others. While left-hand-side errors in variables (LHS EIV) do not bias slope estimates in linear models, such robustness does not carry over to nonlinear models like quantile regression. hausman2021errors examine LHS EIV in quantile regression under the assumption of independence between the treatment and unobserved heterogeneity. This paper extends their analysis to an endogenous quantile regression framework, allowing for dependence between individual heterogeneity and treatment selection, as is typical when the treatment arises as an equilibrium outcome of individual decisions.

We address endogeneity using a control-function approach within a triangular simultaneous equations model. The model consists of an outcome equation and a reduced-form first stage. The outcome equation takes the form of a linear conditional quantile specification with an additive measurement error in the dependent variable. This setup retains the structure in hausman2021errors, while allowing a continuously distributed endogenous regressor to be correlated with the unobserved heterogeneity (i.e., the latent rank). The endogenous regressor is generated in the first-stage equation, where instruments are assumed independent of the reduced-form scalar disturbance. In this framework, the conditional distribution of the endogenous variable given the instruments, which coincides with the distribution of the first-stage disturbance, functions as the control variable imbens2009identification. This control variable thereby restores conditional independence between the endogenous regressor and the latent rank in the outcome equation.

Although the control-function approach effectively isolates the endogenous component, it also introduces additional nuisance parameters due to the dependence between the latent rank and the control variable. In our model, this dependence is naturally captured through a copula with uniform margins. Exploiting this marginal uniformity property, this paper establishes nonparametric identification of the structural quantile coefficient functions, as in the presence of both additive measurement error in the outcome variable and endogeneity. This extends the analysis in hausman2021errors by accommodating endogeneity without imposing any parametric restrictions on the stochastic relationship between the latent rank and the first-stage disturbance.

For estimation, we propose a two-stage sieve maximum likelihood estimator (2SSMLE). We parametrize the nuisance distributions, while sieving over the main functional parameters, namely the quantile coefficient functions. The first stage constructs the control variable. The second stage applies a sieve MLE to recover the distributional effects conditional on both the endogenous regressor and the generated control variable. Under a $\sqrt{n}$-consistent parametric first stage and a mildly ill-posed second stage, the resulting estimator achieves pointwise asymptotic normality with convergence rate of $(n\kappa_{J_n})^{-1/2}$, where $\{\kappa_{J_n}\}$ denotes the decay rate of the minimum eigenvalue of the Hessian matrix as the sieve dimension $J_n\to\infty$, capturing the degree of ill-posedness. Both the convergence rate and the shape of the asymptotic variance are determined by the second stage, while the first-stage error propagates into the variance at first order. When the first stage is estimated nonparametrically, for example using a series estimator, an additional regularization bias arises due to series approximation error. Under an undersmoothing condition on the series dimension $K_n$, which ensures that the first stage approximation bias is asymptotically negligible relative to the second stage rate, the 2SSMLE retains asymptotic normality, thereby permitting inference by nonparametric bootstrap.

The proposed estimator applies to a wide range of empirical settings involving continuous treatments, where endogeneity can be addressed through control functions while the outcome variable is potentially subject to measurement error. A prominent example is the estimation of individual Engel curves, as studied in blundell2007semi and imbens2009identification. In this setting, the latent rank captures heterogeneous consumption preferences, total expenditure is endogenous and survey-reported expenditure shares are prone to misreporting. An empirical investigation along this line will be reported in future work. In this paper, Monte Carlo simulations show that neglecting measurement error can lead to substantial bias, even after correcting for endogeneity, and that the magnitude of this bias depends on both the scale and the shape of the error distribution.

This paper adds to the literature on measurement error in nonlinear models. Although mismeasurement in independent variables has received considerable attention, errors in the dependent variable remain relatively understudied. For discrete outcomes, hausman1998misclassification analyze misclassification in discrete-response models. For continuous outcomes, abrevaya1999semiparametric examines mismeasurement within a general linear index framework, while cosslett2004efficient examine censored regressions. The present work specifically advances the literature on quantile regression. Most existing studies in this domain address measurement error in covariates and restore identification through instruments or repeated measures, including schennach2008quantile, chesher2017understanding, wei2009quantile, firpo2017measurement, and Song2026. In contrast, left-hand-side measurement error in quantile regression has received far less attention. hausman2001mismeasured demonstrates that such errors can produce substantial bias in quantile regression estimates when conditional heteroskedasticity is present. More recently, hausman2021errors establishes that the parameters of a linear quantile regression model with left-hand-side measurement error remain identified and can be consistently estimated when the regressors are independent of the unobserved heterogeneity.\footnote{Alternatively, if there is a collection of sufficiently rich included exogenous variables $\mathbf{Z}$, such that the endogenous regressors $\mathbf{X}$ are independent of the latent rank $U$ conditional on $\mathbf{Z}$, then one can also resort to hausman2021errors. However, such a set $\mathbf{Z}$ is always unavailable in empirical data, thus we need additional information, typically provided by either instrumental variables or by control functions.} Our paper extends their analysis by incorporating endogeneity and introducing a two-stage estimation procedure based on a control-function approach. It further highlights how copula structures can be used to integrate control variables into random-coefficient models. Relatedly, DotySong2023 consider nonparametric identification and estimation of conditional quantiles with LHS measurement error in a dynamic production function framework. Their setting imposes a form of Hicks-neutral productivity, under which the transitory shock and productivity are independent conditional on inputs. This structure effectively renders latent quantile heterogeneity exogenous after conditioning on the control function, placing the model closer in spirit to hausman2021errors after conditioning. Finally, callaway2021distributional study two-sided measurement error in QR models without endogeneity, and propose a three-step EM-type algorithm for a broad class of distributional effect parameters.

The literature on identifying heterogeneous treatment effects with endogenous regressors is extensive and largely follows two approaches: instrumental variables methods, such as chernozhukov2005iv, and control function methods, such as imbens2009identification and newey1999nonparametric. Both strategies, however, break down in the presence of left-hand-side measurement error, as the convolution of the error with unobserved heterogeneity distorts the observed conditional quantiles. We address this identification challenge by combining a control variable with the scalar unobservable through a copula representation. Copulas have been widely used in econometrics to model dependence structures. For instance, among many others, han2017identification employ single-parameter copulas to characterize dependence among unobservables, and callaway2021distributional also adopt copula to model the joint distribution of treatment and outcome. In this paper, we impose no parametric restrictions on the copula in identification. Instead, identification is restored through large-support instruments and the monotonicity inherent in the linear quantile outcome specification.

The rest of the paper is organized as follows. Section (ref) introduces the model specification and develops the nonparametric identification results. Section (ref) presents the two-stage sieve ML estimator. Section (ref) derives its asymptotic properties and establishes the validity of the bootstrap for inference. Section (ref) reports Monte Carlo evidence. Section (ref) concludes the study and outlines directions for future research. Appendix (ref) collects all proofs, additional simulation results, and implementation details of the estimator.

We adopt the following notation: We use upper-case Latin letters for random variables and the corresponding lower-cases for their realizations. Bold Latin letters denote matrices. For a random vector $X$ with support $\mathcal{X}$, $X_p$ denotes its $p$th dimension, and $X_{-p}$ denotes the subvector of $X$ corresponding to all but the $p$th dimension. $\mathbb{E}_n(X_i)=\frac1n\sum_{i=1}^n X_i$ denotes the empirical analogue of $\mathbb{E}[X_i]$. We write $F(\cdot)$ for a CDF, $F_{\cdot|\cdot}(\cdot|\cdot)$ for a conditional CDF, and $f(\cdot)$ for a PDF. For a parameter space $\Theta$, $\Theta^{\circ}$ denotes its interior, and $\mathcal{N}(\theta)\subseteq\Theta$ is a neighborhood of $\theta\in\Theta$. For any $\tau\in(0,1)$, the check function is $\rho_\tau(w)=w(\tau-\mathds{1}(w<0))$.

Model and Identification

Let $Y^*$ be the latent outcome and $Y$ be the observed outcome measured with error. Let $X$ be a random vector with dimension $d_x$ and support $\mathcal{X}$. The model we consider has an outcome equation

where $\varepsilon$ is a scalar measurement error. We adopt a classical mismeasured assumption: $\varepsilon$ is mean-zero, i.i.d., and independent of $Y^*$ and $X$. All unobserved heterogeneity in $Y^*$ is absorbed by a scalar rank $U\sim\mathcal{U}[0,1]$, and $x'\beta_0(\cdot)$ is monotonically increasing for each value $x\in\mathcal{X}$. When $X$ and $U$ are independent, the $\tau$th conditional quantile of $Y^*$ on $X=x$ is $Q_{Y^*|X}(x;\tau)=x'\beta_0(\tau)$, so that $U$ naturally represents the unobserved quantile level of $Y^*$ conditional on $X$.

In practice, however, $U$ may be correlated with $X$, which often incorporates equilibrium outcomes partially determined by individual heterogeneity. Consider a triangular system with a single continuously distributed endogenous variable $X_1$ included in $X=(X_1,Z_1')'$, where $Z_1$ is a vector of exogenous covariates. In general, if $X_1$ becomes independent of $U$ after conditioning on a sufficiently rich set $Z_1$, then endogeneity ceases to be a concern, and $U$ retains the same interpretation. Such conditioning variables are rarely available in practice, and additional information is therefore required.

We obtain this additional information from a control variable. Let $Z=(Z_1',Z_2')'$, where $Z_2$ contains the excluded instrumental variables. The reduced form for $X_1$ is given by

where $h(Z,\eta)$ is strictly monotonic in the scalar disturbance $\eta$ with probability one (w.p.1), and $Z$ is independent of $(\eta,U)$. While we allow for discrete instruments in $Z_2$, $\eta$ is a continuously distributed scalar with a strictly increasing CDF on its support. In this setting, imbens2009identification demonstrate that the uniformly distributed

serves as a control variable: $X$ and $U$ are independent conditional on $V$.\footnote{Our framework is not restricted to a single endogenous regressor. Appendix (ref) presents an extension to multiple endogenous variables within a triangular system. Identification and estimation proceed analogously by employing a vector of control variables $\mathbf{V}$ and a corresponding multivariate copula $F_{U|\mathbf{V}}$.}

In this framework, the coefficient function $\beta_0(\cdot)$ admits multiple interpretations. One perspective, based on the potential outcome framework with a continuous treatment (e.g., chernozhukov2005iv, d2023nonparametric,callaway2024difference), views $\beta_{0,1}(\tau)$ as the $\tau$-th conditional quantile treatment effect (QTE) of $X_1$, that is, $\partial_{x_1} q_{Y^*}(\tau,x_1,z_1)=\beta_{0,1}(\tau)$, holding $X_1=x_1$ and conditioning on $Z_1=z_1$. Alternatively, in a triangular simultaneous equations model, the quantile structural function (QSF) $q_{Y^*}(\tau,x)$, defined as the $\tau$-th quantile of $Y^*_x$, equals $x'\beta_0(\tau)$ under a monotonicity condition lehmann2006nonparametrics,imbens2009identification. In this representation, $\beta_0(\cdot)$ is the quantile coefficient function that parameterizes the structural quantile response of $Y^*$ to covariates.

However, in the presence of both endogeneity and measurement error, pointwise identification of $\beta_0(\tau)$ fails. Consequently, simple quantile regression of $Y$ on $X$ may lead to significant bias. Namely, note that the minimization problem

is generally no longer minimized at the true $\beta_0(\tau)$. First, endogeneity shifts the $\tau$th conditional quantile of $Y^*|X=x$ to $x'\beta_0(\tilde\tau)$, where $\tilde\tau=F^{-1}_{U|V}(\tau)$. Moreover, measurement error further moves the $\tau$th conditional quantile of $Y|X=x$ away from $x'\beta_0(\tilde\tau)$ through convolution with the distribution of $\varepsilon$. Hence, instead of assigning a moment condition for each $\beta_0(\tau)$, we characterize $\beta_0(\cdot)$ using the conditional likelihood function given in ((ref)).

Technically, proposition (ref) shows that the observed conditional distribution $F_{Y|X,V}$ is a convolution of $F_\varepsilon$ and $F_{U|V}$. A formal derivation is provided in Appendix (ref). In the absence of endogeneity (i.e., when $U\perp\!\!\!\!\perp V$), $F_{U|V}(du|v)=F_U(du)=du$, and (ref) reduces to the representation in hausman2021errors:

With an endogenous regressor, in contrast, $F_\varepsilon$ is weighted by the copula $F_{U|V}$, introducing an additional nuisance component, which captures the dependence between the latent heterogeneity $U$ and the endogenous variation in $X$. While hausman2021errors establish nonparametric identification of $(\beta_0(\cdot),f_\varepsilon)$, we proceed to show that $(\beta_0(\cdot),f_\varepsilon,f_{U|V})$ can also be nonparametrically identified under the following assumptions in Section (ref).

Nonparametric Identification

With the above conditions on the parameters, covariates, and the two distributions, we state our main nonparametric identification result as follows.

See Appendix (ref) for the proof. As suggested by hausman2021errors, monotonicity of $x'\beta_0(\cdot)$, together with independence of $\varepsilon$, is critical for securing nonparametric identification. Although the unobservable measurement error effectively blur each quantile curve through a common kernel, the quasi-linear quantile structure, combined with the presence of one continuous regressor that enters strictly monotonically, transforms this into an identifiable deconvolution problem. In our setting, the common support assumption allows us to integrate the copula weights $f_{U|V}$ over $v\in(0,1)$. Because both $U$ and $V$ have uniform marginal distributions, this integration effectively collapses the problem to the exogenous case. In other words, identification of $(\beta_0(\cdot),f_\varepsilon)$ does not depend on specifying a parametric form for $F_{U|V}$. Once $\beta_0(\cdot)$ is identified, the copula $F_{U|V}$ can then be recovered.

Estimation

Based on the nonparametric identification result, this section describes our estimator. We adopt a two-step estimation strategy based on independent and identically distributed (i.i.d.) data $\{Y_i,X_i,Z_i\}_{i=1}^n$. The first step estimates the control variable $V_i$, yielding $\hat V_i$. The second step plugs in the generated variable $\hat V_i$ to construct a sieve maximum likelihood estimator (SMLE). While Theorem (ref) demonstrates nonparametric identification of $(\beta_0(\cdot),f_{0,\varepsilon},f_{0,U|V})$, for estimation, we assume that $f_{0,\varepsilon}$ and $f_{0,U|V}$ are known up to a finite set of parameters,\footnote{In practice, the number of mixture components in each distribution can be increased flexibly to ensure adequate approximation.} and construct a semi-nonparametric estimator targeting $\beta_0(\cdot)$. Specifically, we impose the following assumptions on $(f_{0,\varepsilon},f_{0,U|V})$ in addition to Assumptions (ref) and (ref).

Assumption (ref) (b) holds for a broad class of mean-zero error distributions in the exponential family. Condition (c) accommodates many commonly used copulas, including Gaussian, Student-$t$, Clayton, and Frank. Specifically, the trimming sequence $\{c_n\}$ imposes a lower bound condition in each finite sample as the copula density approaches the boundary of the unit square. Condition (d) ensures uniformly bounded fourth moments of the score functions, so that the information matrix is locally nonsingular when applying the triangular central limit theorem. Finally, following hausman2021errors, conditions (e)-(f) ensure local identification of $(\sigma,\gamma)$ through characteristic-function-based transforms.

Under such parametrization, we adopt the following parameter space and sieve spaces for approximation. Let $\theta_0:=(\beta_0(\cdot),\sigma_0,\gamma_0)\in\Theta_0=M\times\Sigma\times\Gamma$, for $(M,\Sigma,\Gamma)$ satisfying Assumptions (ref)-(ref). Moreover, for each $k\in\{1,\dots,d_x\}$, let $M_k\subseteq \Lambda_{c_k}^{p_k}([0,1])$ with $p_k>\frac12$.\footnote{Following chen2007large, $\Lambda^p_c(\mathcal{X})$ denotes the class of all $p$-smooth real-valued functions on $\mathcal{X}$. A real-valued function $h$ is $p$-smooth $(p=m+\gamma)$ if it is $m$ times continuously differentiable on $\mathcal{X}$ and its differential operator $D^\alpha h$ satisfies a H\"{o}lder condition with exponent $\gamma$ for all $\alpha$ with $[\alpha]=m$ and a dilation $c$.} For any $\theta,\tilde\theta\in \Theta_0$, define the metric $d(\theta,\tilde\theta):=\sqrt{\|\beta-\tilde\beta\|_{2,2}^2+\|\sigma-\tilde\sigma\|_2^2+\|\gamma-\tilde\gamma\|_2^2}$, where $\|\cdot\|_2$ denotes the Euclidean norm and $\|\cdot\|^2_{2,2}$ is the $l_2$-norm of $\|\beta(u)-\tilde\beta(u)\|^2_2$ across $u\in[0,1]$. The sieve spaces are defined as follows.

By definition, the sieve spaces satisfy $\Theta_n\subseteq\Theta_{n+1}\subseteq \Theta_0$ with each $\Theta_n$ compact under the metric $d(\cdot,\cdot)$ and $\inf_{\theta\in\Theta_n} d(\theta,\theta_0)\to0$. Define $\theta_J^*:=\pi_n\theta_0=(\beta_J^*(\cdot),\sigma_0,\gamma_0)$, where $\pi_n$ denotes the projection mapping from $\Theta_0$ to $\Theta_n$. Denote the sieve dimension as $J_n=\tilde{J}_n+r+1$. For each $k=1,\dots,d_x$, the sieve approximation is $\beta_k(u)=b_k' S(u)$, $b_k\in\mathbb{R}^{{J}_n}$, with basis vector $S(u)=(S_1(u),\dots,S_{{J}_n}(u))'$. Let $\mathbf{b} := (b_1',\dots,b_{d_x}') \in \mathbb{R}^{d_x \times{J}_n}$ denote the stacked coefficient vectors, then $x'\beta(u) = \sum_{k=1}^{d_x} x_k b_k' S(u)=x'\cdot \mathbf{b} \cdot S(u)$. The directional derivative of $f(y|x,v;\theta)$ w.r.t $b_k$ along direction of $S(u)$ is given by

Let $\partial_{{b}} f$ denote the stacked gradient in $\mathbb{R}^{d_x {J}_n}$. The score function is given by

with the corresponding information matrix defined as $\mathcal{I}:=\mathbb{E}[(\partial_\theta\log f)(\partial_\theta\log f)']$ and the Hessian matrix as $\mathcal{H}:=\mathbb{E}[\partial_{\theta\theta'}\log f]$.

Our estimation proceeds in two steps. The first step constructs the control variable $V=F_{X_1|Z}=F_\eta$. Following hahn2013asymptotic, we consider both parametric and nonparametric approaches for this step, which clarifies how first-stage estimation error propagates into the second-stage criterion, even when the first-stage estimator achieves $\sqrt{n}$-consistency. In the parametric specification, suppose

where $F_\eta$ is known. Let $\hat\pi$ denote the OLS estimator of $\pi_0$, satisfying $\hat\pi-\pi_0=O_p(n^{-1/2})$. The control variable is then constructed as

Such parametric specifications are often justified when researchers are willing to impose structural assumptions; see, for example, petrin2010control. More generally, we estimate the control variable nonparametrically by approximating the conditional distribution function $F_{X_1|Z}$. Following imbens2009identification, we employ series estimators, which are less sensitive to edge effects that often present in triangular models.\footnote{In particular, the joint density of $(X_1,V)$ may approach zero near the boundary of the support of $V$, thereby limiting the amount of information available in the tails imbens2009identification.} Let $p^{K_n}(z) = (p_{1K_n}(z),\dots,p_{K_nK_n}(z))'$ denote a $K_n\times1$ vector of spline basis functions, then

where $A^{-}$ denotes a generalized inverse of the matrix $A$.

The second stage implements sieve maximum likelihood estimation by replacing the unobserved control variable $V$ in the conditional log-likelihood with its estimate $\hat{V}$. Denote $W_i=(Y_i,X_i',V_i)'$ and $\hat{W}_i=(Y_i,X_i',\hat{V}_i)'$. Let $l(W_i;\theta)=\log f(Y_i|X_i,V_i;\theta)$. The expected conditional log-likelihood function is

By Theorem (ref), the true parameter $\theta_0$ uniquely maximizes $Q(\theta)$ over $\Theta_0$. Let $\hat{Q}(\theta)=\mathbb{E}[l(\hat{W}_i;\theta)]$, with $Q_n(\theta)$ and $\hat{Q}_n(\theta)$ be the corresponding empirical counterparts. The feasible estimator $\hat\theta_n$ is then defined to maximize $\hat{Q}_n(\theta)$ over $\Theta_n$, i.e.,\footnote{For any compact $\hat{\mathcal{V}}_n\subseteq\mathcal{V}$, this integral can be well approximated using standard quadrature rules such as the Gauss-Legendre or Tanh-sinh quadrature.}

Asymptotic Properties

This section establishes consistency and asymptotic normality of the two-step sieve maximum likelihood estimator (2SSMLE), thus permitting nonparametric bootstrap for inference. We distinguish between two cases: (i) the first step is parametric, yielding $\hat\theta_n^{p}$ in Section (ref); (ii) the first step is nonparametric, yielding $\hat\theta_n^{np}$ in Section (ref).

Parametric First Stage

We begin by imposing some sufficient regularity conditions on the parametric first stage.

Under the preceding conditions, the generated control variable in ((ref)) satisfies $|\hat{V}_i-V_i|=O_p(n^{-1/2})$, i.e., it is pointwise $\sqrt{n}$-consistent. Combined with consistency of the sieve approximation in the second stage when $V$ is observed, this yields consistency of $\hat\theta_n^{p}$.

The proof is given in Appendix (ref). Next, we establish the $L_2$-convergence rate of $\hat\theta_n^{p}$ to the finite-dimensional projection of the true parameter, i.e., $\|\hat\theta_n^{p}-\theta_J^*\|_2^2$. This rate determines how quickly the sieve dimension $J_n$ may grow in the presence of an ill-posed inverse problem. As noted by chen2006efficient, the convergence rate of sieve MLE depends on the smoothness level of the likelihood function. Following fan1991optimal, we characterize the smoothness of the error distribution through the tail behavior of its characteristic function $\phi_\varepsilon(s)$ as $|s|\to\infty$.\footnote{Although in our model the error distribution is convolved with a copula component, the copula density is uniformly bounded over the trimmed domain $\mathcal{V}_n$, thus does not affect the effective smoothness classification.} Namely, the distribution is supersmooth of order $\lambda>0$ if $|\phi_\varepsilon(s)|\asymp |s|^\lambda\exp(-|s|^\lambda/\varsigma)$ as $s\to\infty$, for some constant $\varsigma > 0$. It is ordinary smooth of order $\lambda>0$ if $|\phi_\varepsilon(s)|\asymp |s|^{-\lambda}$ as $s\to\infty$. Although the smoothness of the error distribution does not affect consistency of $\hat\theta_n^p$, we require it to be ordinary smooth to establish the asymptotic normality result. While this condition rules out the exactly normal distribution (which is supersmooth of order 2), hausman2021errors (hausman2021errors) demonstrate that even small perturbations from normality typically render the characteristic function ordinary smooth of finite order. Consequently, a sieve based on mixtures of ordinary-smooth densities remains sufficiently flexible to approximate supersmooth distributions while effectively avoiding ill-posedness in finite samples.

We introduce one final assumption, adapted from hausman2021errors, which facilitates the derivation of convergence rate for the finite-dimensional distributional parameters.

Assumption (ref) requires that, when weighted by the true copula density, there is sufficient variation in $(x,v)$ within the local neighborhood so that the characteristic function is nondegenerate. We now establish the $L_2$-convergence rate of $\|\hat\theta_n^p - \theta_J^*\|_2^2$ as follows.

The proof is provided in Appendix (ref). The derived convergence rates impose restrictions on how quickly the sieve dimension $J_n$ may diverge, because the second-stage estimation involves an ill-posed inverse problem, under which the minimum eigenvalue of the Hessian can approach zero as the number of basis functions increases. To ensure the validity of the high dimensional central limit theorem, the convergence rates established above must dominate the decay rate of this eigenvalue (See more details in Theorem (ref)). This requirement motivates a careful analysis of the rates in Lemma (ref). In particular, combining Lemma (ref) with consistency of $\hat\theta_n^p$ yields $\|\hat\sigma-\sigma_0\|_2=o_p(n^{-1/4})$ for the error distribution parameter. By contrast, the convergence rates for the quantile coefficient functions and the copula parameter are slower and depend directly on the first-stage estimation, since the true $V$ is unobserved. As a result, a slow-converging first stage may limit the range of admissible growth rates of $J_n$ required for weak convergence.\footnote{The distinction between $(a)$ and $(b$-$c)$ stems from independence between $\varepsilon$ and $x'\beta$, resulting in the factorization $\phi_{y\mid x,v}(s)=\phi_\varepsilon(s) \cdot \phi_{x'\beta\mid v}(s)$, which is then used to determine the convergence rates.} We now present the asymptotic normality of $\hat\theta_n^p$.

The proof is provided in Appendix (ref). Theorem (ref) shows that the 2SSMLE achieves pointwise asymptotic normality provided that the sieve dimension $J_n$ grows at an appropriate rate. In particular, $J_n^{2r+2}/n\to\infty$ requires $J_n$ to grow sufficiently fast so that the sieve approximation bias is adequately controlled. Meanwhile, $J_n^{\frac{\lambda}{\lambda+1}}n^{\frac{-1}{4(\lambda+1)}}/\kappa_{J_n}\to0$ restricts $J_n$ from growing too quickly, thereby preserving weak convergence in the growing dimensional setting and preventing the problem from becoming severely ill-posed. By the information identity, $\mathcal{I}(\theta_0,v) = -\mathcal{H}(\theta_0,v)$, so the minimum eigenvalues of $\mathcal{H}(\theta_0,v)$, $\mathcal{H}(\theta_J^*,v)$, and $\mathcal{I}(\theta_0,v)$ share the same decay rate. Consequently, the estimator is $(n\kappa_{J_n})^{-1/2}$ consistent, where $\kappa_{J_n}$ captures the degree of ill-posedness in the inverse problem. Using the same measure of ill-posedness as in hausman2021errors, if the problem is mildly ill-posed, so that $\kappa_{J_n} \asymp J_n^{-\delta}$, then weak convergence holds under the above conditions if $(1+\delta)\lambda + \delta < \frac{1}{2}(r+1)$. By contrast, if $\kappa_{J_n} \asymp \exp(-J_n^\delta)$, corresponding to the severely ill-posed case, then weak convergence fails even though consistency is maintained.

Theorem (ref) also indicates that the contribution of the first stage is not simply negligible, even though it attains a parametric rate that dominates the mildly ill-posed second stage. Instead, the generated regressor error propagates through the second step in a first-order way: It enters through a $J_n\times d_z$ loading, which embeds the first-stage randomness into the same high-dimensional sieve score space. More specifically, the asymptotic variance $\Omega_J$ can be decomposed into three components, corresponding to the randomness arising from both stages as well as their interaction. The first-and second-stage randomness are given by

where $G_\pi(\theta_0,\pi_0):=\mathbb{E}\big[\partial_v\psi(W_i;\theta_0,\pi_0)\cdot \partial_\pi V(X_i,Z_i;\pi_0)\big]$ and $\Sigma_\pi:=\mathbb{E}[Z_iZ_i']^{-1}\mathbb{E}[\eta_i^2]$. In general, the lack of orthogonality renders $G_\pi\neq0$. However, in the special case without endogeneity, where $f_c(u|v) = f_c(u) = 1$ and $\partial_v \psi = 0$, we have $G_\pi(\theta_0,\pi_0) = 0$. In this case, the variance reduces to that found in hausman2021errors, i.e., $\Omega_{\text{sec}}$, and the first-stage contribution disappears. Moreover, by the definition of $\kappa_{J_n}$, the eigenvalues of $\Omega_{\text{sec}}$ are bounded. We show that the full asymptotic variance satisfies $0\le\Omega_{\text{sec}}\le \Omega_J\le (1+C)\Omega_{\text{sec}}$, for some uniform constant $C>0$. Hence, while the first-stage contribution is not asymptotically negligible, the second stage determines both the convergence rate $\sqrt{n\kappa_{J_n}}$ and the shape of the asymptotic variance.

Nonparametric First Stage

We now present the asymptotic theory for the series first-stage estimator defined in ((ref)). As emphasised by imbens2009identification, the convergence rate of $\hat V$ depends on the smoothness of the conditional distribution $F_{X_1\mid Z}$. We hereby adopt the following assumption.

Under Assumptions (ref) and (ref), Lemma 11 of imbens2009identification shows that

This rate depends on both the series dimension $K_n$ and the smoothness level of $F_{X_1|Z}$. It consists of a variance term, $(K_n/n)$, and a squared bias term, $(K_n^{1-2d_1/d_z})$.\footnote{The additional factor $K_n$ in the squared bias arises because the predicted values $\hat{V}_i$ are obtained from regressions whose dependent variables vary across observations.} Let $\alpha_n:=(K_n/n+K_n^{1-2d_1/d_z})^{1/2}$. As in the parametric case in Section (ref), when $\alpha_n \to 0$, the estimator $\hat\theta_n^{np}$ is consistent and satisfies the following convergence rates.

The proofs of the two results above parallel those of Theorem (ref) and Lemma (ref). For the series basis functions, we impose Assumption 2 of newey1997convergence. Specifically, let $p^K(z)$ denote the original $K$-dimensional vector of basis functions, and let $P^K(z)$ denote its normalized transformation used in estimation. Assume there exists a nonsingular constant matrix $B$ such that $P^K(z)=Bp^K(z)$ for all $K$; the smallest eigenvalue of $\mathbb{E}[p^K(Z_i)p^K(Z_i)']$ is bounded away from zero uniformly in $K$, and $\sup_{z\in\mathcal{Z}}\|P^K(z)\|\le \zeta_0(K)$ for a sequence $\zeta_0(K)$ satisfying $\zeta_0(K)^2 K/n\to0$ as $n\to\infty$. It is well acknowledged that for power series $\zeta_0(K)\le CK$ and for splines $\zeta_0(K)\le C K^{1/2}$, for some generic positive constant $C$. Denote the conditional distribution function $F_{X_1|Z}$ as $F$, with the true distribution as $F_0$ and the series estimates as $\hat{F}$. Then the pointwise asymptotic distribution of $\hat\theta_n^{np}$ is as follows.

See Appendix (ref) for the proof. The above theorem extends the asymptotic normality result to the case of a nonparametric first stage. We show that the first-stage error can be governed so that it does not break down the limiting distribution, while it introduces an additional regularization bias term arising from series approximation error, which enters linearly through the score derivative $\partial_v\psi$. The resulting influence function includes a U-statistic projection that accounts for the dependence between the nonparametric control residuals and the second-stage score. As in Theorem (ref), the influence of the first stage on both the asymptotic variance and the bias term disappears only in the absence of endogeneity, where $\partial_v\psi=0$. In practice, the above result guarantees flexible adoption of nonparametric control-function estimators without sacrificing valid inference, provided the regularity and undersmoothing conditions are satisfied. The empirical selection of the optimal series dimension $K_n$ and sieve dimension $J_n$ is beyond the scope of this paper.

Inference via Bootstrap

Although Theorem (ref) and (ref) establishes asymptotic normality of the two-stage sieve MLE, direct variance estimation can be cumbersome depending on the distributional specification. In practice, we therefore recommend inference using nonparametric pairs bootstrap, following the same principle as in chen2003estimation. Specifically, for each bootstrap replication $b=1,\dots,B$, draw a bootstrap sample $\{(Y_i^{*(b)},X_i^{*(b)},Z_i^{*(b)})\}_{i=1}^n$ from the empirical distribution of the data. We re-estimate the full two-step procedure to obtain bootstrap control variables $\hat V_i^{*(b)}$ and the second-stage sieve MLE $\hat\theta_n^{*(b)}$. We then construct pointwise $(1-\alpha)$ confidence intervals using the bootstrap standard errors: $\hat\beta_k(\tau)\pm z_{1-\alpha/2}\,\widehat{\mathrm{se}}_{k}(\tau)$, where $z_{1-\frac{\alpha}{2}}$ denotes the $(1-\frac{\alpha}{2})$ standard normal quantile and \[ \widehat{\mathrm{se}}_{k}(\tau) = \left[ \frac{1}{B-1} \sum_{b=1}^B \left( \hat\beta_k^{*(b)}(\tau) - \bar\beta_k^*(\tau) \right)^2 \right]^{1/2}, \qquad \bar\beta_k^*(\tau) = \frac{1}{B} \sum_{b=1}^B \hat\beta_k^{*(b)}(\tau), \] for each quantile index $\tau$ of interest. Analogous intervals can be constructed for finite-dimensional parameters. Such procedure automatically accounts for estimation uncertainty in both stages. We establish its validity in Appendix (ref).

Simulation

In this section, we examine the finite-sample performance of the proposed estimator using Monte Carlo simulations. The data-generating-process (DGP) is given by

The latent rank is generated as $U=\Phi(A)$, where $A\sim\mathcal{N}(0,1)$ and $\Phi$ denotes the standard normal CDF, so that $U\sim\mathcal{U}[0,1]$. The coefficient functions are specified as

The model includes one endogenous regressor $X_1$ and one exogenous control $X_2$, such that $X_1\perp\!\!\!\!\perp X_2$ and $X_2\sim \text{LN}(0,1)$. We generate $X_1\sim\text{LN}(0,1)$ by:

where $Z,\eta\sim \mathcal{N}(0,1)$ and $Z\perp\!\!\!\!\perp \eta$. The parameter $\delta\in (0,1)$ controls instrument strength and is set to $0.5$ in our baseline DGP. In the baseline scenario, endogeneity is introduced through a Gaussian copula linking the two latent variables $U$ and $V$. Specifically,

where $B\sim \mathcal{N}(0,1)$, $B\perp\!\!\!\!\perp A$, and $\rho\in(-1,1)$ governs the degree of endogeneity. Note that $X_1$ is strictly increasing in $\eta$, and $V$ is the CDF of $\eta$. Finally, we draw mean-zero measurement error $\varepsilon$ from a mixed Gaussian distribution following hausman2021errors:

We illustrate the two-step sieve maximum likelihood estimator (2SSMLE) using both parametric and nonparametric first-stage estimation. In the parametric case, assume the true specification is fully known, so that $V$ can be estimated by directly plugging in OLS estimates. More generally, in the nonparametric case, we run a spline regression of $X_1$ on $(Z,X_2)$ and then construct $\hat{V}_i=\hat{F}(X_{1i}|Z_i,X_{2i})$ as in ((ref)). The second step plugs in first-step estimates $\hat{v}$ and solves the empirical log-likelihood maximization problem specified in ((ref)). We use B-splines of order $r$ with $J_n$ grid knots to approximate the unknown coefficient functions at $\{\beta_k(\tau_j)\}_{j=1}^{J_n}$. Given the properties of the measurement error and copula distributions, the integral within the log-likelihood can be well approximated by Gauss-Legendre quadrature with $q$ nodes per interval on the same partition used for the B-spline basis.

Despite the sophistication of the second-stage likelihood and optimization problem, the proposed estimator remains computationally efficient under our implementation. Specifically, we solve the optimization problem using a gradient-based constrained interior-point method, with analytical gradients supplied to stabilize the descent path and facilitate convergence. Stochastic gradient descent may also be employed to mitigate the risk of convergence to local stationary points. Additional implementation details are provided in Appendix (ref).

Our 2SSMLE estimator simultaneously corrects the bias from endogeneity and measurement error. In contrast, standard quantile regression (QR) is biased from both sources. The control-function QR (hereafter denoted CFQR), such as lee2007endogeneity, addresses only endogeneity, and the sieve MLE (hereafter denoted SMLE) of hausman2021errors addresses only measurement error but ignores possible endogeneity. For CFQR in our simulations, we implement an adaptation of lee2007endogeneity: the control variable $\hat V$ is estimated from our first stage, and the second runs a partial linear QR of $Y$ on $(X_1,X_2,\hat{V})$.\footnote{Specifically, define $P_k(w)=[x_1,x_2,p_1(v),\dots,p_k(v)]'$ for a power-series basis $\{p_k:k=1,2,\dots\}$. With a trimming function $t(w)=\mathds{1}\{w\in\mathcal{W}\}$ to limit unduly values lee2007endogeneity,blundell2007censored, CFQR solves $\min_\theta S_{nk}(\theta)=\frac1n\sum_{i=1}^n t(\hat{W}_i)\rho_\tau[Y_i-P_k(\hat{W}_i)'\theta]$ for any $\tau\in(0,1)$ and $\hat{W}_i=(X_{1i},X_{2i},\hat{V}_i)'$.}

We first illustrate how standard estimators are biased across quantiles when both endogeneity and measurement error are present. We then shut down one source of bias at a time to assess the cost of accounting for features that may be absent. Finally, we examine the robustness of 2SSMLE to misspecification and evaluate bootstrap validity.

Figure (ref) plots the average point estimates $\hat\beta(\tau_j)$ at 20 grid knots to visually compare the bias patterns across methods. We first focus on the slope coefficient function $\beta_1(\cdot)$ (Figure (ref), middle panel), while the intercept function $\beta_0(\cdot)$ shown in the left panel exhibits analogous patterns. CFQR (red squares) is biased due to the measurement error $\varepsilon$: lower quantiles biased upwards and upper quantiles downwards, reflecting the influence of the two spikes of the $\varepsilon$-distribution at $-3$ and $3$, as shown in Figure (ref). SMLE (yellow diamonds), unaffected by EIV, tracks the shape of the true $\beta_1$ but is upward biased due to endogeneity. Simple QR (blue circles) suffers from both issues. By comparison, 2SSMLE (purple crosses) exhibits much smaller bias at all quantiles: its average absolute bias is about 4% of the true $\beta_1$, versus 43% for QR, 24% for CFQR, and 35% for SMLE. For the slope coefficient of the exogenous covariate $X_2$ (Figure (ref), right panel), there is no endogeneity issue. Consequently, QR and CFQR closely align with each other and are affected only by measurement error. Likewise, SMLE and 2SSMLE closely overlap--both track the true coefficient function well.

Such bias worsens as endogeneity exacerbates: In Figure (ref), both QR and SMLE results in larger deviation as the correlation between $U$ and $V$ increases from $\rho=0.5$ to $\rho=0.9$. Note that the impact of EIV depends on its magnitude relative to the true outcome. In Figure (ref), when $\varepsilon\sim\mathcal N(0,1)$, the EIV-induced distortion is small, making SMLE and QR almost indistinguishable, both primarily reflecting endogeneity rather than measurement error. Finally, we compare 2SSMLE with the default estimators in settings where either only endogeneity or measurement error is absent, thereby illustrating the cost of conservatively accounting for features that may not actually be present. As shown in Table (ref), 2SSMLE generally outperforms SMLE when there is no endogeneity and CFQR when there is no EIV, although it produces slightly larger bias in the slope coefficient function than CFQR due to misspecification in the likelihood model.\footnote{Although our identification framework does not formally allow $\varepsilon\stackrel{a.s.}{=}0$, the 2SSMLE can approximate a kernel-smoothed indicator function through the error distribution, making it applicable to multivalued discrete measurement error as well.}

In the absence of knowledge about the true distributional families, we illustrate the flexibility of the Gaussian–mixture specification by showing that it can approximate a variety of alternative EIV distributions following hausman2021errors, including student-$t$ and Laplace. Numerical details of the aforementioned results can be found in Tables (ref) to (ref). Table (ref) shows the bootstrap standard errors and average coverage of our bootstrap confidence interval. At a significance level $\alpha=0.05$, the overall coverage is close to $95\%$. Without further specification, all the above implementation of 2SSMLE uses a parametric first stage. Finally, we repeat the above exercises using series-estimated control functions. Table (ref) shows that results parallel the pattern seen in the parametric framework. More details on mean bias and MSE of different methods, bootstrap standard errors, bootstrap coverage and additional simulation results are provided in Appendix (ref).

Conclusion

This paper studies quantile regression in the presence of an endogenous regressor and additive measurement error in the dependent variable. After isolating the endogenous treatment from the latent rank with a control variable, we establish nonparametric identification of the conditional quantile coefficient function. Assuming the nuisance distributions are known up to some finite-dimensional parameters, we propose a plug-in two-step sieve maximum likelihood estimator (2SSMLE), which combines a control-function first stage with a sieve ML second stage, incorporating the generated control variable through copula weights. The first-stage control function can be estimated flexibly by series regression. When the series dimension and the sieve dimension grow at appropriate rates, the estimator is consistent and asymptotically normal, with convergence rates governed by the degree of ill-posedness. Monte Carlo simulations show that our estimator substantially outperform existing methods in settings where both endogeneity and measurement error may be present. In future work, a substantive empirical application—such as revisiting Engel curve estimation in blundell2007semi—would provide an opportunity to examine how outcome mismeasurement interacts with endogeneity in practice, and to assess the empirical consequences of our proposed corrections for distributional treatment effects.