EconBase
← Back to paper

Bias correction for quantile regression estimators

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.

48,071 characters · 14 sections · 45 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Bias correction for quantile regression estimators

abstract\linespread{1.2} We study the bias of classical quantile regression and instrumental variable quantile regression estimators. While being asymptotically first-order unbiased, these estimators can have non-negligible second-order biases. We derive a higher-order stochastic expansion of these estimators using empirical process theory. Based on this expansion, we derive an explicit formula for the second-order bias and propose a feasible bias correction procedure that uses finite-difference estimators of the bias components. The proposed bias correction method performs well in simulations. We provide an empirical illustration using Engel's classical data on household food expenditure. JEL Classification: C21, C26. Keywords: instrumental variables, higher-order stochastic expansion, Bahadur-Kiefer expansion, finite-difference estimators, mixed integer linear programming (MILP), Engel curve.

\onehalfspacing

Introduction

Many interesting empirical applications of classical quantile regression (QR) koenker1978regression and instrumental variable quantile regression (IVQR) chernozhukov2005iv,chernozhukov2006instrumental feature small sample sizes, which can arise either as a result of a limited number of observations or when estimating tail quantiles, or both chernozhukov2005annals,elsner2008increasing,chernozhukov2011inference,adrian2016covar,adrian2019vulnerable. QR and IVQR estimators are nonlinear and can thus exhibit substantial biases in small samples. In this paper, we theoretically characterize these biases and develop a feasible bias correction procedure.

To study the biases, we start by deriving a higher-order stochastic expansion of the classical QR and exact IVQR estimators.\footnote{We define exact IVQR estimators as estimators that exactly minimize a norm of the sample moment conditions. Such estimators can be obtained using mixed-integer programming (MIP) methods chen2018exact,zhu2019learning. See Appendix (ref).} Such an expansion is needed because the higher-order terms contribute nonzero biases while the first-order term does not. We derive explicit expressions and uniform (in the quantile level) rates for the components in the expansion, building on the empirical process arguments of ota2019quantile. This expansion can be thought of as a refined Bahadur-Kiefer (BK) representation of the estimator that decomposes the nonlinear component into terms up to the order $O_p\bigPar{n^{-1}}$ and a $O_p\bigPar{n^{-5/4}\sqrt{\log n}}$ remainder. We also derive a uniform BK representation for generic IVQR estimators after a feasible 1-step Newton correction (see Appendix (ref)).

Using the stochastic expansion, we study the bias of QR and exact IVQR estimators. We derive a bias formula based on the leading terms up to order $O_p\bigPar{n^{-1}}$ in the expansion, which we refer to as the second-order (asymptotic) bias. This approach of focusing on the moments of the leading terms in the stochastic expansions is standard in the literature.\footnote{See, for example, nagar1959bias,newey2004higher,kato2012asymptotics,galvao2016smoothed, kaplan2017smoothed,hahn2023efficient, among others.} The second-order bias formula provides an approximation of the actual bias that yields a feasible correction. Our results explicitly account for the bias due to nonzero sample moments at the estimator. Our proof strategy is different from the generalized function heuristic used in the existing literature phillips1991shortcut,lee2017second,lee2018second, which does not account for all the terms in the second-order bias (see Section (ref) and Appendix (ref) for further discussion and examples). The missing terms can be important bias contributors in practice, as we document in the empirical application in Section (ref).

A feasible bias correction procedure then follows from the second-order bias formula. We propose finite-difference estimators of all the components in the formula. These estimators admit higher-order expansions that allow us to select bandwidth rates. In particular, our finite-difference estimator for the Jacobian coincides with powell1986censored's classical estimator, and the bandwidth rate we derive coincides with the AMSE optimal bandwidth choice in kato2012asymptotic. We show that the resulting (analytically) bias-corrected estimator has zero second-order bias. This result is in contrast to the commonly-used bootstrap bias correction approaches horowitz2001bootstrap, which may not capture the higher-order bias terms of quantile estimators knight2003second.

We evaluate the performance of our bias correction procedure in a Monte Carlo simulation study. The simulations show that the theoretical (infeasible) bias formula describes well the second-order bias of classical QR and exact IVQR. We find that the proposed feasible bias correction can effectively reduce the bias in many cases. The gains from bias correction are particularly prominent in settings with large bias such as for the IVQR estimators under endogeneity. The impact of the bias correction on the root MSE (RMSE) is rather small and ambiguous.

We illustrate the bias correction approach by revisiting the relationship between food expenditure and income based on the original engel1857 data koenker1982robust,koenker2001quantile.\footnote{The QR approach was used more recently to analyze changes in heterogeneity in the income elasticity over time of various consumer spending categories in taylor2009consumer.} Our results highlight the importance of bias correction in empirical applications with small sample sizes. Specifically, we find that the second-order bias of classical QR is empirically relevant: it can be larger than 50% of the standard error, which is substantial given that this is the second-order bias.

\paragraph{Roadmap.} The remainder of the paper is organized as follows. Section (ref) describes the model and the estimators. Section (ref) provides our main theoretical results. Section (ref) presents the Monte Carlo simulation results. Section (ref) contains the empirical application. Section (ref) concludes. All the proofs and some additional details are given in the Appendix.

Model and estimators

Consider a setting with a continuous outcome variable $Y$, a $(k\times 1)$ vector of covariates $W$, and a $(k\times 1)$ vector of instruments $Z$. We assume throughout that $k$ is fixed. Every observation $(Y_i,W_i,Z_i)$, $i=1,\dots,n$, is jointly drawn from a distribution $P$. We assume that $(Y_i,W_i,Z_i)$ is i.i.d., and we will sometimes suppress the index $i$ to lighten up the notation. The parameter of interest $\thetaTrue\in\Theta\subset \R^k$ is defined as a solution to the following unconditional quantile moment restrictions,

equation[equation omitted — 122 chars of source]

We consider two cases: (i) classical QR, where $ Z=W$ koenker1978regression, and (ii) linear IVQR, where $Z\neq W$ in general chernozhukov2006instrumental,chernozhukov2008instrumental.

The classical QR estimator of $\thetaTrue$ is a solution to the following convex minimization problem,

equation[equation omitted — 164 chars of source]

where $\rho_\tau(u)=u(\tau-1\{u<0\})$ is the check function koenker2005quantile and $\emp$ denotes the sample average, i.e., the expectation with respect to the empirical measure. For IVQR, we consider estimators that exactly minimize the $p$-norm of the sample moments,

equation[equation omitted — 110 chars of source]

where $p\in[1,\infty]$ and $\gthat(\theta)\bydef \emp (1\{Y\leq W^{\prime}\theta\}-\tau)Z$. This class of exact IVQR estimators includes GMM, which corresponds to $p=2$ as in chen2018exact for just-identified models, and the estimator proposed by zhu2019learning, which corresponds to $p=\infty$. The cases $p=1$ and $p=\infty$ have computationally convenient mixed integer linear programming (MILP) representations, while the MILP formulation for $p=2$ has many more decision variables. In our Monte Carlo simulations, we use $p=1$ for computational convenience (see Appendix (ref)).

We use the notation $\gt(\theta) \bydef \expect(1\{Y\leq W^{\prime}\theta\}-\tau) Z$ for the unconditional moment restrictions as a function of $\theta\in\Theta$, and write $G(\theta)\bydef \partial_\theta \expect Z1\{Y\leq W^{\prime}\theta\}=\partial_\theta g_{\tau}(\theta)$ for its derivative. We maintain the following standard identification assumptions.

assumption[Identification] \begin{enumerate} • $\thetaTrue$ is the unique solution to $\gt(\theta)=0$ over a compact set $\Theta \subset \R^k$, and $\thetaTrue$ is in the interior of $\Theta$ for all $\tau\in [\varepsilon,1-\varepsilon]$ for some $\varepsilon>0$. • The Jacobian $G(\thetaTrue)$ has full rank for all $\tau\in (0,1)$. \end{enumerate}

As noted by chernozhukov2006instrumental, “compactness [of the parameter space $\Theta$] is not restrictive in micro-econometric applications” (p.502). Throughout the paper, we use the short notation $G$ for $G(\thetaTrue)$ whenever it does not lead to ambiguity.

We impose the following smoothness assumptions on the conditional density and its derivatives. Such assumptions are standard in the literature on higher-order properties of quantile estimators ota2019quantile.

assumption[Conditional density] The conditional density of $Y_i$ given $(W_i,Z_i)$, $f_{Y}(y|w,z)$, exists, is a.s.\ three times continuously differentiable on $supp(Y)$, and there exists a constant $\bar{f}$ such that $|f^{(r)}_{Y}(y|w,z)|\le \bar{f}$ for all $(y,w,z)\in supp(Y)\times supp(W)\times supp(Z)$, where $r=0,1$ and $f_Y^{(r)}(\cdot|w,z)$ is the $r$-th derivative of $f_Y(\cdot|w,z)$.

In our theoretical analysis of the bias, we will often work with a related object, the conditional density $f_{\varepsilon_{\tau}}(e|W,Z) \bydef f_{Y}(e+W'\thetaTrue|W,Z)$ of the quantile residual $ \varepsilon_{\tau} \bydef Y- W^\prime\thetaTrue$.

Finally, we assume that the regressors and the instruments have bounded higher-order moments.

assumption[Regressors and instruments] There exists constants $m<\infty$ and $\gamma\geq6$ such that $\expect |W_j|^\gamma \leq m$ and $\expect |Z_j|^\gamma \leq m$ for all $j=1,\dots,k$.

Assumption (ref) guarantees the existence of all relevant moments of the terms involving $W$ and $Z$ in the higher-order derivatives of the moment conditions and the bias correction (e.g., $\expect Z_\ell W_j W_q W_r$). The power parameter $\gamma$ (as we show below) determines the upper bound on the rate at which the norm of the sample moment functions converges to zero --- higher $\gamma$ implies faster convergence to zero.

Asymptotic theory for bias correction

To derive a bias correction procedure, we follow the approach of nagar1959bias and focus on the bias of the leading terms in the asymptotic stochastic expansion of the estimator.

Stochastic expansion of quantile regression estimators

The classical first-order asymptotic theory for quantile regression estimators koenker1978regression,angrist2006quantile,chernozhukov2006instrumental,kaplan2017smoothed,kaido2021decentralization is based on the following leading term,

equation[equation omitted — 162 chars of source]

For correctly specified models, $\hat\xi_\tau$ is an infeasible unbiased estimator of $\thetaTrue$. However, because feasible quantile estimators are nonlinear, they generally have a nonzero higher-order bias. The following theorem provides a characterization of the terms in a stochastic expansion of $\thetahat$ up to order $O_p\bigPar{n^{-1}}$ (ignoring logarithmic terms).

To state the result, we introduce some additional notation. For $\theta\in\Theta$, define the auxiliary functions $g^\circ(\theta) \bydef \expect1\{Y\leq W^{\prime}\theta\}Z$, $B^\circ_n(\theta) \bydef \sqrt{n} (\emp 1\{Y\leq W^{\prime} \theta \} Z - g^\circ (\theta))$, and $B_n(\theta) \bydef B^\circ_n(\theta) - B^\circ_n(\thetaTrue)$.\footnote{The processes $B_n(\theta) $ and $B^\circ_n(\theta) $ take values in the space $\ell^{\infty}(\Theta)$ of bounded functions on $\Theta$.} Also, denote the Hessian of the $j$-th moment function $g_j$ by $\partial_\theta G_j(\theta) \bydef \partial_{\theta } \partial_{ \theta} g_j(\theta) $. Finally, for any $x\in \R^{k}$, denote by $x^\prime \partial_\theta G(\theta) x$ the vector with components $x^\prime \partial_\theta G_j(\theta) x$, $j=1,\dots,k$.

thmSuppose that Assumptions (ref)--(ref) hold. Consider $\thetahat = \thetaLp$ obtained from program (ref) for some $p\in [1,\infty]$ or $\thetahat = \hat\theta_{\tau,QR}$. Then \begin{equation*} \thetahat = \hat\xi_\tau+ G^{-1} (\thetaTrue)\bigBra{\gthat(\thetahat) - \frac{B_n(\thetahat)}{\sqrt{n}} - \frac{1 }{2} ( \hat\xi_\tau-\thetaTrue)^\prime \partial_\theta G(\thetaTrue) (\hat\xi_\tau-\thetaTrue) } + R_{n,\tau}, \end{equation*} where \begin{align} \sup_{\tau \in [\varepsilon,1-\varepsilon]}\|\gthat(\thetahat)\| &= \begin{cases} O_p\bigPar{\frac{ 1 }{ n }}, if \thetahat=\thetaQR,\\ O_p\bigPar{\frac{ \log n }{ n^{1-\frac{2}{\gamma}} }}, if \thetahat=\thetaLp, \end{cases} \\ \sup_{\tau \in [\varepsilon,1-\varepsilon]}\| B_n(\thetahat) \|&= O_p\bigPar{\frac{\sqrt{\log n}}{n^{1/4}}} , \\ \sup_{\tau \in [\varepsilon,1-\varepsilon]}\| \hat\xi_\tau-\thetaTrue\| &= O_p\bigPar{\frac{ 1}{n^{1/2}}},\\ \sup_{\tau \in [\varepsilon,1-\varepsilon]}\|R_{n,\tau}\|&=O_p\bigPar{\frac{\sqrt{\log n}}{n^{5/4}}}. \notag \end{align}

The proof of Theorem (ref) builds on the empirical process arguments in Lemma 3 of ota2019quantile and uses the maximal inequality in Corollary 5.1 of chernozhukov2014gaussian.

The rates in Theorem (ref) are uniform in the quantile level $\tau$. Uniformity is important in theory and practice because QR and IVQR methods are particularly powerful when used to analyze the entire quantile process.

Classical QR and linear IVQR are motivated by the linearity of the true conditional quantile function. In practice, linearity may be restrictive, especially when it is imposed at each $\tau\in (0,1)$. While the result in Theorem (ref) remains valid if linearity fails, the interpretation is more complicated in this case. Specifically, if the true conditional quantiles are nonlinear, then the expansion in Theorem (ref) applies to the pseudo-true value defined via moment condition (ref). For classical QR, this pseudo-true value can be interpreted as the minimizer of a weighted mean-squared error loss function angrist2006quantile. Note further that the result in Theorem (ref) holds for every fixed $\tau\in (0,1)$. Thus, if linearity holds at a given $\tau$, we can interpret the expansion at this $\tau$ under the correct specification without requiring linearity at the other quantile levels.

The expansion in Theorem (ref) can be thought of as a refined Bahadur-Kiefer (BK) expansion. Different from standard BK expansions, we do not bundle together all the higher-order terms zhou1996direct,ota2019quantile. Notice that the dominant nonlinear term in the BK expansion, $ n^{-1/2}G^{-1}B_n(\thetahat)$, has order $O_p\bigPar{ n^{-3/4}\sqrt{\log n} }$ (Equation (ref)). Theorem 3 in knight2002comparing shows that for classical QR with discrete covariates, this term converges in distribution to a zero mean random process. Therefore, we explicitly extract the higher-order terms up to order $O_p\bigPar{n^{-1}}$ (ignoring logarithmic terms) from the BK remainder. As we will show in the following sections, these higher-order terms admit feasible counterparts.

rem[Alternative approach for deriving stochastic expansions] portnoy2012nearly proposed an alternative approach for deriving a stochastic expansion of classical QR estimators. This approach yields bounds on the precision of a nonlinear Gaussian approximation of order $O_p\bigPar{n^{-1}\log^{5/2}n}$. As we will show, the expansion in Theorem (ref) yields a bias formula for both QR and IVQR estimators that admits a feasible implementation. The results in portnoy2012nearly are specific to classical QR, and it is not clear to us whether these results can be used for bias correction using a Nagar-style approach. \qed
rem[General IVQR estimators] While we focus on exact IVQR estimators in the main text, the results in Theorem (ref) can be used to obtain a uniform BK expansion for general 1-step corrected IVQR estimators. See Appendix (ref) for details. \qed

Bias formula for exact estimators

Following common practice nagar1959bias,kaplan2017smoothed, for a generic estimator $\hat\gamma$, we define the second-order bias $\bias(\hat\gamma)$ as the bias of the leading terms in the stochastic expansion of $\hat\gamma$ up to the order $O_p\bigPar{n^{-1}}$. This second-order bias can be interpreted as an approximation of the actual bias that works with arbitrarily high probability in large samples.

Before stating the result, we observe that under our Assumption (ref), the moment condition (ref) is equivalent to

equation*[equation* omitted — 85 chars of source]

Thus, we can characterize the QR and IVQR estimators using the moment function $g_\tau^*(\theta) \bydef \expect[(1\{- Y\leq W^{\prime}\theta\}-(1-\tau)) Z ] $ with sample analog $\gthat^\ast(\theta)$. The following theorem characterizes the second-order bias in terms of $\gthat(\hat\theta) $ and $ \gthat^*(-\hat\theta) $. As in Section (ref), we define $\partial_\theta G_j(\theta) \bydef \partial_{\theta } \partial_{ \theta} g_j(\theta) $ for all $\theta\in\Theta$.

thmSuppose that Assumptions (ref)--(ref) hold. Consider $\thetahat = \thetaLp$ obtained from program (ref) for some $p\in [1,\infty]$ or $\thetahat = \hat\theta_{\tau,QR}$. Then the second-order bias is \begin{equation} \bias(\thetahat) = G^{-1}(\thetaTrue) \bigBra{ \frac{1}{2}\expect \bigPar{ \gthat(\thetahat) - \gthat^*(-\thetahat) } - \frac{\kappa_{\tau}}{n} - \frac{1}{2n } Q^\prime vec( \Omega_\tau ) }, \end{equation} where \begin{align*} \kappa_{\tau } &\bydef \bigPar{\tau -\frac{1}{2}} \expect f_{ \varepsilon_\tau} (0 |W, Z) Z W^{\prime} G^{-1}Z, \\ \Omega_\tau &\bydef \operatorname{Var} [Z(1\{Y\le W' \thetaTrue\}-\tau)],\\ Q & is a matrix with columns Q_j \bydef \operatorname{vec}\bigBra{ ( G^{-1})^\prime \partial_\theta G_j(\thetaTrue) G^{-1} }, \,\, j=1,\dots,k. \end{align*}

The second-order bias formula (ref) has three components.

The first component, $ G^{-1}\expect \bigPar{ \hat g_\tau(\thetahat) - \hat g_\tau^*(-\thetahat) }/2 $, captures the bias from the sample moments not being zero at the estimator. This term is not equal to zero in general, as we illustrate based on a simple example in Appendix (ref).

The second component, $n^{-1}G^{-1} \kappa_\tau $, appears because of the discontinuity in the sample moment functions. The term $\kappa_\tau$ reflects the dependence between the sample moments and the linear influence of a single observation on $\thetahat$. Notice that the first two components combined correspond to the second-order bias of the terms $G^{-1} (\thetaTrue)(\gthat(\thetahat) - n^{-1/2} B_n(\thetahat)) $ in Theorem (ref).

The last component, $(2n)^{-1}G^{-1}Q^\prime vec(\Omega)$, stems from the non-uniformity of the conditional distribution of $Y$ given $(W,Z)$. This term corresponds to the term $ G^{-1} (\thetaTrue) (\hat\xi_\tau-\thetaTrue)^\prime \partial_\theta G(\thetaTrue) (\hat\xi_\tau-\thetaTrue)/2 $ in Theorem (ref). Similar terms are typically present in most nonlinear estimators with nonzero Hessian of the score function rilstone1996secondorder.

To illustrate the approximate bias formula, consider an order statistic of $Y\sim \text{Uniform}(0,1)$ (corresponding to our framework with $Z=W=1$), for which an exact bias formula is available ahsanullah2013introduction. We show in Appendix (ref) that the precision of the second-order bias formula in this case is $O\bigPar{n^{-2}}$, which is smaller than the order of the remainder term in the stochastic expansion of Theorem (ref). Figure (ref) illustrates the precision of the asymptotic formula by comparing it to the actual bias of the order statistic.

figure[figure omitted — 444 chars of source]

It is interesting to compare our results to the higher-order bias analysis of non-smooth estimators based on the generalized functions heuristic phillips1991shortcut. In recent work, lee2017second,lee2018second derived a second-order bias formula for classical QR and IVQR under the assumption that the sample moments are zero at the estimator so that the first term of the bias formula vanishes. We show in Appendix (ref) that this term is non-negligible even in simple cases (see also Figure (ref) in Section (ref)).

Feasible bias correction

The bias formula suggests the following feasible bias-corrected estimator,

equation[equation omitted — 238 chars of source]

where $\hat G$, $\hat\kappa_{\tau}$, $\hat Q$, and $\hat\Omega$ are estimators of $G$, $\kappa_{\tau}$, $Q$, and $\Omega$, respectively, satisfying the following consistency requirement.

assumption[Consistency of component estimators] The estimators $\hat G$, $\hat\kappa_{\tau}$, $\hat Q$, and $\hat\Omega$ are consistent for $G$, $\kappa_{\tau}$, $Q$, and $\Omega$, respectively. Moreover, $\hat G-G = o_p\bigPar{(n^{1/3}\log n)^{-1}}$.

In Section (ref), we propose finite-difference estimators for which Assumption (ref) holds under Assumptions (ref)--(ref). Assumption (ref) could also be verified for other nonparametric estimators of the bias components.

The next theorem shows that the second-order bias of the bias-corrected estimator is zero.

thmSuppose that Assumptions (ref)--(ref) hold. Consider $\thetahat = \thetaLp$ obtained from program (ref) for some $p\in [1,\infty]$ or $\thetahat = \hat\theta_{\tau,QR}$. Then the feasible bias correction eliminates the second-order bias, $\bias \bigPar{ \hat\theta_{bc} } = 0.$

Note that the requirement that $\hat G-G = o_p\bigPar{(n^{1/3}\log n)^{-1}}$ in Assumption (ref) is necessary to ensure that the contribution of the product of the sample moments, which are $O_p\bigPar{n^{-1+2/\gamma}\log n}$ with $\gamma\geq 6$ (as required by Assumption (ref)) by Theorem (ref), and the estimation error in $G^{-1}$ can be omitted in computing the second-order bias. This condition is only required for IVQR estimators. For classical QR estimators, the sample moments are of order $O_p\bigPar{n^{-1}}$, so that consistency of $\hat G$ at any rate of convergence suffices for Theorem (ref).

Finite difference estimators of bias components

To implement the bias correction, we need estimators of $G$, $\kappa_{\tau}$, $ Q$, and $ \Omega$ that satisfy Assumption (ref). The variance matrix $\Omega$ can be estimated using the analogy principle,

equation*[equation* omitted — 122 chars of source]

All other bias components take the form of derivatives. Therefore, we leverage our theoretical results on the properties of the sample moments to develop a unified finite-difference framework for estimating these components.

Under Assumptions (ref) and (ref), the Jacobian is $G=\expect f_{\varepsilon_{\tau}}(0|W,Z)ZW'$ and the Hessian consists of gradients of the components of $G$, i.e. $\partial_\theta G_{i,j}(\thetaTrue)=\expect f^{(1)}_{\varepsilon_{\tau}}(0|W,Z) Z_i W_j W $, where $i,j=1,\dots,k.$ This suggests the following analog estimators.

The $(i,j)$-th component of $G$ can be estimated using powell1986censored's estimator

equation[equation omitted — 165 chars of source]

where $h_{1,n} \to 0$ is a bandwidth. Denote by $e_\ell$ the $\ell$-th unit vector in $\R^k$, where $\ell=1,\dots, k$. The derivative of the $(i,j)$-th component of $G$ in the direction $e_\ell$ (i.e. the second partial derivative of $\gt$) can be estimated as the symmetric first difference of (ref),

align*[align* omitted — 220 chars of source]

where $h_{2,n} \to 0$ is a (potentially different) bandwidth. For $\kappa_{\tau}$, the finite difference sample analog is

equation*[equation* omitted — 196 chars of source]

where $h_{3,n} \to 0$ is a bandwidth. Finally, $Q$ can be estimated by the sample analog matrix $\hat Q$ with columns

equation*[equation* omitted — 135 chars of source]

where $ \widehat{\partial_\theta G_j} $ is the matrix with elements $\widehat{({\partial_\theta G}_{i,j})}_\ell$ for $i,\ell=1,\dots,k$.

The next lemma establishes the consistency of the estimators of the bias components. It implies that these estimators satisfy the high-level conditions in Assumption (ref). Moreover, it provides nearly remainder-optimal bandwidth rates, i.e., the rates that yield the fastest convergence rates of the stochastic remainder terms of the corresponding stochastic expansions (up to logarithmic terms).

lemMTSuppose that Assumptions (ref)--(ref) hold. Then the nearly remainder-optimal bandwidth rates are $h_{1,n} \propto n^{-1/5}$, $h_{2,n} \propto n^{-1/7}$, $h_{3,n} \propto n^{-1/5}$, and, under these bandwidth rates, \begin{align*} \hat G &= G + O_p\bigPar{\frac{\sqrt{\log n}}{n^{2/5}}},\\ \widehat{({\partial_\theta G}_{i,j})}_\ell &= {({\partial_\theta G}_{i,j})}_\ell + O_p\bigPar{\frac{\sqrt{\log n}}{ n^{2/7} }},\\ \hat Q_j &= Q_j + O_p \bigPar{ \frac{\sqrt{\log n}}{n^{2/7}} }, \\ \hat \kappa_{\tau} &= \kappa_{\tau} + O_p\bigPar{\frac{\sqrt{\log n}}{n^{2/5}}},\\ \hat \Omega_{\tau} &= \Omega + O_p \bigPar{ \frac{1}{\sqrt{n}} }. \end{align*} Moreover, the convergence rate for $\hat G$ is uniform in $\tau\in[\varepsilon,1-\varepsilon]$.

Proposition 1 in kato2012asymptotic shows that the remainder rate for $\hat{G}$ in Lemma (ref), $h_{1,n} \propto n^{-1/5}$, is the AMSE-optimal rate.\footnote{ We conjecture that analogous AMSE-optimality results could be established for estimators $\widehat{({\partial_\theta G}_{i,j})}_\ell$ and $\hat Q_j $, but leave this extension for future work.}

To implement the finite difference estimators in practice, one needs to choose constants in addition to the bandwidth rates. A natural approach would be to use non-parametric methods such as cross-validation to estimate the optimal bandwidth parameter. However, such non-parametric methods typically require large sample sizes simonoff1996smoothing, which makes them impractical for our purposes. Another common approach, which we take here, is to propose rule-of-thumb choices that account for the scale of the data. To be robust to heavy-tailed distributions, we suggest using a robust measure of dispersion. Specifically, we propose a rule-of-thumb choice in form $h_{1,n}=A_{G}\cdot n^{-1/5}$, $h_{2,n}=A_Q\cdot n^{-1/7}$, $h_{3,n}=A_\kappa\cdot n^{-1/5}$, where, for all $j\in \{G,Q,\kappa\}$,

equation*[equation* omitted — 86 chars of source]

and $\widehat{\operatorname{MAD}}_\tau$ is the estimated median absolute deviation of the $\tau$-quantile residuals. The constant $1.48$ is chosen because for the normal distribution 1.48 times the median absolute deviation is equal to the standard deviation.\footnote{chernozhukov2013inference suggest using $\operatorname{IQR}/1.35$ as a robust estimate of the scale parameter. For the normal distribution, this choice coincides with $1.48 \operatorname{MAD}$.} We suggest choosing $\tilde{A}_G=\tilde{A}_\kappa=2 $ and $\tilde{A}_Q=1.5 $ based on our simulation evidence, as we describe in more detail in Section (ref).

Simulation evidence

In this section, we evaluate the performance of our feasible bias correction procedure in a Monte Carlo simulation study.

Simulation design

We consider data-generating processes (DGPs) inspired by the simulations in andrews2016conditional. The outcome is generated according to the following location-scale model

align[align omitted — 69 chars of source]

where $W_i = \Phi(\tilde{W}_i)$, $Z_i = \Phi(\tilde{Z}_i)$, $U_i = F_U^{-1}\left(\Phi(\tilde{U}_i)\right)$, $(\tilde{W}_i,\tilde{Z}_i, \tilde{U}_i)\sim N(0,\Sigma)$, $\Sigma_{11}=\Sigma_{22}=\Sigma_{33}=1$, $\Sigma_{23}=0$, and $\Phi$ is the standard normal CDF. Hence, in all the designs, both the regressors and the instruments are $\operatorname{Uniform} [0,1]$. We consider six DGPs that differ with respect to the error distribution $F_U$ and whether or not $W$ is exogenous.

center[center omitted — 831 chars of source]

In Appendix (ref), we consider two additional DGPs to assess the impact of the strength of the instrument.

Bandwidth choice

An important practical issue is the choice of the bandwidths $h_{1,n}$, $h_{2,n}$, and $h_{3,n}$. We use the (data-dependent) rule-of-thumb bandwidth choice described in Section (ref). To implement this rule-of-thumb bandwidth choice, we need to choose the constants $\tilde{A}_j$, $j\in \{G,Q,\kappa\}$. We found that across the DGPs considered in our simulations, $\tilde{A}_G=\tilde{A}_\kappa=2$ and $\tilde{A}_Q=1.5$ perform well in terms of bias across all designs, especially when the bias is large.\footnote{If more information about the DGP is available, it is possible to further refine the tuning parameter. For example, under DGP1 where $Q=0$, choosing smaller values for $A_Q$ results in a better alignment with the infeasible formula. However, absent such information, we consider the proposed rule of thumb to be a reasonable compromise when one cares about the performance across different designs, especially those with larger biases (e.g., DGP3--DGP6).}

We compare the performance of the feasible bias correction based on the rule-of-thumb bandwidth choice to the performance of the corresponding infeasible bias correction based on the true $G$, $Q$, and $\kappa$. For IVQR, the exact formulas are not available. Therefore, we use numerically computed values based on 10 million observations and tuning choice ${A}_G={A}_Q={A}_\kappa=1$.

Results

We focus on the performance of bias correction for $\tau\in \{0.25, 0.5, 0.75\}$. Figure (ref) shows the impact of bias correction for the exogenous DGP1--DGP3 with $n=100$. Figure (ref) shows the corresponding results for the endogenous DGPs (DGP4--DGP6). We use classical QR of $Y$ on $W$, implemented via the linear programming formulation in Appendix (ref), for DGP1--DGP3 and IVQR, implemented using the MILP formulation in Appendix (ref), for DGP4--DGP6.

The main findings can be summarized as follows. First, the QR estimators based on the exogenous designs (DGP1--DGP3) exhibit lower biases than the IVQR estimators based on the corresponding endogenous designs (DGP4--DGP6). For both estimators, the biases tend to be the largest for the designs with heavy tails (DGP3 and DGP6). Second, the feasible bias correction based on the rule-of-thumb bandwidth choice reduces the bias of QR and IVQR estimators in many cases. The bias reductions are the most notable when the biases of the original estimators are large, which is when bias correction is most needed. Finally, the infeasible bias correction reduces the bias in most cases, underscoring the usefulness of the proposed theory. See Appendix (ref) for additional simulation evidence on the impact of the sample size and instrument strength on the performance of bias correction.

Next, we investigate the impact of bias correction on the RMSE of the estimators. The proposed bias correction approach is designed to reduce the bias but is not theoretically guaranteed to reduce the RMSE. Figure (ref) reports the results for DGP1 and DGP4 with $n=100$. Overall, the impact of bias correction on the RMSE is rather small. While the infeasible bias correction can slightly increase the RMSE for DGP1, it decreases the RMSE across all quantile levels for DGP4. The feasible bias correction with the rule-of-thumb bandwidth slightly increases the RMSE in most cases. Appendix Figures (ref) and (ref) show the results for the other DGPs, and Appendix Figure (ref) shows the corresponding results for the Mean Absolute Deviation, an alternative measure of risk.

Finally, we study the impact of the bias correction on the coverage probability of the standard confidence intervals. To focus on the impact of bias correction, we use the same standard errors based on the rule-of-thumb bandwidth for the original and the bias-corrected estimator. As a result, by construction, the bias correction does not affect the length of the confidence intervals. Figure (ref) shows the empirical coverage for DGP1 and DGP4 before and after bias correction. The feasible bias correction based on the rule-of-thumb bandwidth leads to higher coverage accuracy in the majority of cases but can lead to some undercoverage at the median. Thus, while bias correction may not improve the RMSE in small samples, it can lead to more accurate inferences. The results for the other DGPs are in Appendix Figures (ref) and (ref), and Appendix Figure (ref) presents the results for $n=200$.

figure[figure omitted — 978 chars of source]
figure[figure omitted — 1,171 chars of source]
figure[figure omitted — 813 chars of source]
figure[figure omitted — 864 chars of source]

Empirical application

The second-order bias matters most in applications with small sample sizes. We therefore illustrate our bias correction approach using the classical dataset of engel1857, analyzed by koenker1982robust and koenker2001quantile, among others. The data contain information on annual income and food expenditure (in Belgian francs) for $n=235$ Belgian working-class households and are obtained from the R package quantreg quantreg_package. One feature of these data is the growing dispersion of the outcome variable (food expenditure) as a function of the regressor (income) koenker2001quantile, which is similar to our Monte Carlo designs. We divide the values of income and food expenditure by $1000$ so that the unit of measurement becomes a thousand Belgian francs. This makes the scale of intercept and slope parameters comparable.

figure[figure omitted — 905 chars of source]

We estimate classical QRs of food expenditure ($Y$) on income ($W$) and a constant (blue dots). The bias-corrected estimates (gold squares) are obtained using the recommended rule-of-thumb bandwidth choice with $(\tilde{A}_G,\tilde{A}_Q,\tilde{A}_\kappa)=(2,1.5,2)$. Figure (ref) presents the results. Panel (a) compares the classical and the bias-corrected QR estimates. The results suggest that the impact of bias correction is more pronounced around the median and in the tails. The magnitude of the differences between the original and bias-corrected estimates can be larger than 50% of the standard errors, which is substantial given that we are focusing on the second-order bias.

Panel (b) shows the individual contributions of the different components to the overall second-order bias. We can decompose the bias correction term as follows:

align*[align* omitted — 310 chars of source]

The main takeaway from the bias decomposition is that while all three components play a role, the sample moment term (i) and especially the Hessian term (iii) can be large and account for most of the bias. The $\kappa$-term (ii) is smaller overall and only matters in the tails.

Conclusion

We demonstrate that classical QR and IVQR estimators can exhibit a non-negligible second-order bias. We characterize this bias theoretically and use this characterization to derive a novel analytical bias correction method. The proposed feasible bias correction reduces the bias of QR and IVQR estimators across a variety of settings at a very low computational cost. However, there is scope for further improving its performance when the sample size is very small and the instruments are weak. The simulation performance of the infeasible bias correction based on the population version of our theoretical bias formula suggests that exploring alternative feasible bias correction approaches is a promising direction for future research. For example, one could consider regularization approaches or explore imposing parametric assumptions to improve the estimation of the bias components.

Acknowledgments

We are grateful to the Editor (Xiaohong Chen), the Associate Editor, and three anonymous referees, as well as Victor Chernozhukov, Zheng Fang, Dalia Ghanem, Jiaying Gu, Marc Henry, Keisuke Hirano, Nail Kashaev, Roger Koenker, Vladimir Koltchinskii, Michal Kolesar, Simon Lee, Blaise Melly, Hyungsik Roger Moon, Hashem Pesaran, Joris Pinkse, Wolfgang Polonik, Stephen Portnoy, Geert Ridder, Andres Santos, Davide Viviano, Yuanyuan Wan, and seminar participants at UC Berkeley, UC Davis, UC Los Angeles, University of Toronto, and USC for valuable comments. All errors and omissions are our own.

\part*{Online appendix}