The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
48,071 characters
Bias correction for quantile regression estimators
\maketitle
\begin{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.
\medskip
\noindent \textbf{JEL Classification:} C21, C26.
\medskip
\noindent \textbf{Keywords:} instrumental variables, higher-order stochastic expansion, Bahadur-Kiefer expansion, finite-difference estimators, mixed integer linear programming (MILP), Engel curve.
\end{abstract}
\newpage
\onehalfspacing
\section{Introduction}\label{sec:intro}
Many interesting empirical applications of classical quantile regression (QR) \citep{koenker1978regression} and instrumental variable quantile regression (IVQR) \citep{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 \citep[e.g.,][]{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 \citep[e.g.,][]{chen2018exact,zhu2019learning}. See Appendix \ref{app:LP_MILP}.}
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 \citet{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{app: Stochastic expansion of general quantile regression estimators}).
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 \emph{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, \citet{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 \citep{phillips1991shortcut,lee2017second,lee2018second}, which does not account for all the terms in the second-order bias (see Section \ref{sec: Bias formula for exact estimators} and Appendix \ref{app: Illustration of approximate bias formula in univariate case} 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{sec: empirical application}.
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 \citet{powell1986censored}'s classical estimator, and the bandwidth rate we derive coincides with the AMSE optimal bandwidth choice in \citet{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 \citep{horowitz2001bootstrap}, which may not capture the higher-order bias terms of quantile estimators \citep[see][]{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 \citet{engel1857} data \citep[e.g.,][]{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 \cite{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{sec:setup_background} describes the model and the estimators.
Section \ref{sec: Asymptotic theory for bias correction} provides our main theoretical results.
Section \ref{sec:MC} presents the Monte Carlo simulation results.
Section \ref{sec: empirical application} contains the empirical application.
Section \ref{sec:conclusion} concludes. All the proofs and some additional details are given in the Appendix.
\section{Model and estimators}
\label{sec:setup_background}
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,
\begin{equation}
\expect[(1\{Y\leq W^{\prime}\thetaTrue\}-\tau) Z ] = 0, \quad \tau\in (0,1).\label{eq:unconditionalMoment}
\end{equation}
We consider two cases: (i) classical QR, where $ Z=W$ \citep{koenker1978regression}, and (ii) linear IVQR, where $Z\neq W$ in general \citep{chernozhukov2006instrumental,chernozhukov2008instrumental}.
The classical QR estimator of $\thetaTrue$ is a solution to the following convex minimization problem,
\begin{equation}
\hat\theta_{\tau,QR} \in \operatorname{argmin}_{\theta\in\Theta} \emp \rho_\tau(Y-W^\prime \theta ), \label{eq:QR_check_function_minimization}
\end{equation}
where $\rho_\tau(u)=u(\tau-1\{u<0\})$ is the check function \citep{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,
\begin{equation}
\thetaLp \in \operatorname{argmin}_{\theta\in \Theta} ||\gthat(\theta)||_p, \label{eq:GMM_Lp}
\end{equation}
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 \citet{chen2018exact} for just-identified models, and the estimator proposed by \citet{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{app:LP_MILP}).
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.
\begin{assumption}[Identification] \label{ass:identification}
\text{ }
\begin{enumerate}
\item $\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$.\label{ass:identification_global}
\item The Jacobian $G(\thetaTrue)$ has full rank for all $\tau\in (0,1)$.\label{ass:identification_full_rank}
\end{enumerate}
\end{assumption}
As noted by \citet{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 \citep[e.g.,][]{ota2019quantile}.
\begin{assumption}[Conditional density]\label{ass: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)$.
\end{assumption}
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$.
\label{REPLY:R1.1.part1}
Finally, we assume that the regressors and the instruments have bounded higher-order moments.
\begin{assumption}[Regressors and instruments]\label{ass:regressors_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$.
\end{assumption}
Assumption \ref{ass:regressors_instruments} 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.
\section{Asymptotic theory for bias correction}
\label{sec: Asymptotic theory for bias correction}
To derive a bias correction procedure, we follow the approach of \cite{nagar1959bias} and focus on the bias of the leading terms in the asymptotic stochastic expansion of the estimator.
\subsection{Stochastic expansion of quantile regression estimators}
\label{sec:BK_IVQR}
The classical first-order asymptotic theory for quantile regression estimators \citep[e.g.,][]{koenker1978regression,angrist2006quantile,chernozhukov2006instrumental,kaplan2017smoothed,kaido2021decentralization} is based on the following leading term,
\begin{equation}
\hat\xi_\tau \bydef \thetaTrue - G^{-1}(\thetaTrue) \emp Z\bigPar{1\{Y\le W^\prime\thetaTrue\} - \tau }. \label{eq:definitionOfTheta_1}
\end{equation}
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$.
\begin{thm} \label{thm:stochasticExpansion}
Suppose that Assumptions~\ref{ass:identification}--\ref{ass:regressors_instruments} hold.
Consider $\thetahat = \thetaLp$ obtained from program \eqref{eq:GMM_Lp} 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 }}, \text{ if } \thetahat=\thetaQR,\\
O_p\bigPar{\frac{ \log n }{ n^{1-\frac{2}{\gamma}} }}, \text{ if } \thetahat=\thetaLp,
\end{cases} \label{eq:ghat_supbound} \\
\sup_{\tau \in [\varepsilon,1-\varepsilon]}\| B_n(\thetahat) \|&= O_p\bigPar{\frac{\sqrt{\log n}}{n^{1/4}}} \label{eq:Bn_supbound}, \\
\sup_{\tau \in [\varepsilon,1-\varepsilon]}\| \hat\xi_\tau-\thetaTrue\|
&= O_p\bigPar{\frac{ 1}{n^{1/2}}},\label{eq:DonskerLinTerm}\\
\sup_{\tau \in [\varepsilon,1-\varepsilon]}\|R_{n,\tau}\|&=O_p\bigPar{\frac{\sqrt{\log n}}{n^{5/4}}}. \notag
\end{align}
\end{thm}
The proof of Theorem \ref{thm:stochasticExpansion} builds on the empirical process arguments in Lemma 3 of \citet{ota2019quantile} and uses the maximal inequality in Corollary 5.1 of \citet{chernozhukov2014gaussian}.
The rates in Theorem \ref{thm:stochasticExpansion} 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{thm:stochasticExpansion} 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{thm:stochasticExpansion} applies to the pseudo-true value defined via moment condition \eqref{eq:unconditionalMoment}. For classical QR, this pseudo-true value can be interpreted as the minimizer of a weighted mean-squared error loss function \citep{angrist2006quantile}. Note further that the result in Theorem \ref{thm:stochasticExpansion} 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{thm:stochasticExpansion} 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 \citep[as opposed to, for example,][]{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 \eqref{eq:Bn_supbound}).
Theorem 3 in \citet{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.
\begin{rem}[Alternative approach for deriving stochastic expansions]
\citet{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{thm:stochasticExpansion} yields a bias formula for both QR and IVQR estimators that admits a feasible implementation.
The results in \citet{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
\end{rem}
\begin{rem}[General IVQR estimators]
While we focus on exact IVQR estimators in the main text, the results in Theorem \ref{thm:stochasticExpansion} can be used to obtain a uniform BK expansion for general 1-step corrected IVQR estimators. See Appendix \ref{app: Stochastic expansion of general quantile regression estimators} for details. \qed
\end{rem}
\subsection{Bias formula for exact estimators}
\label{sec: Bias formula for exact estimators}
Following common practice \citep[e.g.,][]{nagar1959bias,kaplan2017smoothed}, for a generic estimator $\hat\gamma$, we define the \emph{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{ass:density}, the moment condition \eqref{eq:unconditionalMoment} is equivalent to
\begin{equation*}
\expect[(1\{- Y\leq W^{\prime}(-\thetaTrue)\}-(1-\tau)) Z ] = 0.
\end{equation*}
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{sec:BK_IVQR}, we define $\partial_\theta G_j(\theta) \bydef \partial_{\theta } \partial_{ \theta} g_j(\theta) $ for all $\theta\in\Theta$.
\begin{thm}\label{thm:bias} Suppose that Assumptions~\ref{ass:identification}--\ref{ass:regressors_instruments} hold. Consider $\thetahat = \thetaLp$ obtained from program \eqref{eq:GMM_Lp} 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 ) },\label{eq:biasFormula}
\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 &\text{ 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*}
\end{thm}
The second-order bias formula \eqref{eq:biasFormula} 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{app: Illustration of approximate bias formula in univariate case}.
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{thm:stochasticExpansion}.
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{thm:stochasticExpansion}.
Similar terms are typically present in most nonlinear estimators with nonzero Hessian of the score function \citep[see, for example,][]{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 \citep[e.g.,][]{ahsanullah2013introduction}.
We show in Appendix \ref{app: Illustration of approximate bias formula in univariate case} 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{thm:stochasticExpansion}.
Figure \ref{fig:exact_asymptotic_uniform_quantile} illustrates the precision of the asymptotic formula by comparing it to the actual bias of the order statistic.
\begin{figure}[ht]
\begin{center}
\caption{Comparison to exact formula in univariate case}
\includegraphics[width=0.7\textwidth]{sections/figures/biasAnalytic.pdf}
\label{fig:exact_asymptotic_uniform_quantile}
\end{center}
\footnotesize{\textit{Notes:} Exact (circles) and second-order (crosses) biases, scaled by $n$, as functions of quantile level $\tau$ for $\thetahat = Y_{(\intPart{\tau n})}$, where $Y\sim \text{Uniform}(0,1)$, $n=10$.}
\end{figure}
It is interesting to compare our results to the higher-order bias analysis of non-smooth estimators based on the generalized functions heuristic \citep[e.g.,][]{phillips1991shortcut}.
In recent work, \citet{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{app: Illustration of approximate bias formula in univariate case} that this term is non-negligible even in simple cases (see also Figure~\ref{fig:empirical_application} in Section~\ref{sec: empirical application}).
\subsection{Feasible bias correction} \label{sec:feasibleCorrection}
The bias formula suggests the following feasible bias-corrected estimator,
\begin{equation}
\hat\theta_{bc}= \thetahat - \frac{1}{2}\hat{G}^{-1} \bigBra{\gthat(\thetahat)-\gthat^*(-\thetahat)} + \frac{1}{n} \hat G^{-1} \bigBra{ \hat \kappa_{\tau } + \frac{1}{2} \hat Q^\prime vec( \hat \Omega ) } , \notag
\end{equation}
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.
\begin{assumption}[Consistency of component estimators]
\label{ass:consistent_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}}$.
\end{assumption}
In Section \ref{sec: finite difference estimators of bias components}, we propose finite-difference estimators for which Assumption \ref{ass:consistent_estimators} holds under Assumptions \ref{ass:identification}--\ref{ass:regressors_instruments}.
Assumption \ref{ass:consistent_estimators} 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.
\begin{thm}\label{thm:biascorrection} Suppose that Assumptions \ref{ass:identification}--\ref{ass:consistent_estimators} hold. Consider $\thetahat = \thetaLp$ obtained from program \eqref{eq:GMM_Lp} 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.$
\end{thm}
Note that the requirement that $\hat G-G = o_p\bigPar{(n^{1/3}\log n)^{-1}}$ in Assumption \ref{ass:consistent_estimators} 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{ass:regressors_instruments}) by Theorem \ref{thm:stochasticExpansion}, 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{thm:biascorrection}.
\subsection{Finite difference estimators of bias components}
\label{sec: 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{ass:consistent_estimators}.
The variance matrix $\Omega$ can be estimated using the analogy principle,
\begin{equation*}
\hat \Omega_{\tau} \bydef \emp [Z(1\{Y\le W'\thetahat\}-\tau)-\emp Z(1\{Y\le W'\thetahat\}-\tau)]^2.
\end{equation*}
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{ass:identification} and \ref{ass:density}, 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 \citet{powell1986censored}'s estimator
\begin{equation}
\hat G_{i,j}= \emp \left[ \frac{1\{Y\le W'\thetahat+ h_{1,n}\} - 1\{Y\le W'\thetahat- h_{1,n}\} }{2h_{1,n}}Z_i W_j \right], \label{eq:G-hat-ij}
\end{equation}
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 \eqref{eq:G-hat-ij},
\begin{align*}
\widehat{({\partial_\theta G}_{i,j})}_\ell &= \emp \left[ \frac{1\{Y\le W'\thetahat+ h_{2,n}\} - 2 \cdot 1\{Y\le W'\thetahat\} + 1\{Y\le W'\thetahat- h_{2,n}\}}{ h_{2,n}^2} Z_i W_j W_\ell \right],
\end{align*}
where $h_{2,n} \to 0$ is a (potentially different) bandwidth. For $\kappa_{\tau}$, the finite difference sample analog is
\begin{equation*}
\hat \kappa_{\tau} = \bigPar{\tau-\frac{1}{2}} \emp \left[ \frac{1\{Y\le W'\thetahat+ h_{3,n}\} - 1\{Y\le W'\thetahat- h_{3,n}\} }{2h_{3,n}} Z W^{\prime} \hat G^{-1}Z \right],
\end{equation*}
where $h_{3,n} \to 0$ is a bandwidth. Finally, $Q$ can be estimated by the sample analog matrix $\hat Q$ with columns
\begin{equation*}
\hat Q_j \bydef vec\bigBra{ ( \hat G^{-1})^\prime \widehat{\partial_\theta G_j} \hat G^{-1} }, \quad j=1,\dots,k,
\end{equation*}
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{ass:consistent_estimators}.
Moreover, it provides \emph{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).
\begin{lemMT}\label{thm: finite difference}
Suppose that Assumptions~\ref{ass:identification}--\ref{ass:regressors_instruments} 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]$.
\end{lemMT}
Proposition 1 in \cite{kato2012asymptotic} shows that the remainder rate for $\hat{G}$ in Lemma \ref{thm: finite difference}, $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 \citep[e.g.][]{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\}$,
\begin{equation*}
A_j= \tilde{A}_j \cdot 1.48 \cdot \widehat{\operatorname{MAD}}_\tau ,
\end{equation*}
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{\cite{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{subsec:BandwithChoice}.
\section{Simulation evidence}
\label{sec:MC}
In this section, we evaluate the performance of our feasible bias correction procedure in a Monte Carlo simulation study.
\subsection{Simulation design}
We consider data-generating processes (DGPs) inspired by the simulations in \citet{andrews2016conditional}.
The outcome is generated according to the following location-scale model
\begin{align}
Y_i= W_i+( 0.5 + W_i)U_i,\quad i=1,\dots,n, \notag
\end{align}
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.
\begin{center}
\begin{tabular}{ l l l}
\hline
\hline
DGP1 (Uniform, exogenous) & $F_U(u)=\int_{-\infty}^u 1\{t\in[0,1]\}dt$ & $\Sigma_{12}=1$, $\Sigma_{13}=0$ \\
\hline
DGP2 (Triangular, exogenous) & $F_U(u)=\int_{-\infty}^u 2 t 1\{t\in[0,1]\}dt$& $\Sigma_{12}=1$, $\Sigma_{13}=0$ \\
\hline
DGP3 (Cauchy, exogenous) & $F_U(u)=\int_{-\infty}^u \frac{1}{\pi(1+(4 t)^2)}dt$ & $\Sigma_{12}=1$, $\Sigma_{13}=0$\\
\hline
DGP4 (Uniform, endogenous) &$F_U(u)=\int_{-\infty}^u 1\{t\in[0,1]\}dt$ & $\Sigma_{12}=0.75$, $\Sigma_{13}=0.25$\\
\hline
DGP5 (Triangular, endogenous) &$F_U(u)=\int_{-\infty}^u 2 t 1\{t\in[0,1]\}dt$ & $\Sigma_{12}=0.75$, $\Sigma_{13}=0.25$\\
\hline
DGP6 (Cauchy, endogenous) &$F_U(u)=\int_{-\infty}^u \frac{1}{\pi(1+(4 t)^2)}dt$ & $\Sigma_{12}=0.75$, $\Sigma_{13}=0.25$\\
\hline
\hline
\end{tabular}
\end{center}
In Appendix \ref{sec:additionalFigures}, we consider two additional DGPs to assess the impact of the strength of the instrument.
\subsection{Bandwidth choice}\label{subsec:BandwithChoice}
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{sec: finite difference estimators of bias components}. 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$.
\subsection{Results}
We focus on the performance of bias correction for $\tau\in \{0.25, 0.5, 0.75\}$. Figure \ref{fig:bias_dgp13} shows the impact of bias correction for the exogenous DGP1--DGP3 with $n=100$. Figure \ref{fig:bias_dgp13_iv} 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{app:LP_MILP}, for DGP1--DGP3 and IVQR, implemented using the MILP formulation in Appendix \ref{app:LP_MILP}, 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{sec:additionalFigures} 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{fig:rmse14} 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{fig:rmse25} and \ref{fig:RMSEstable} show the results for the other DGPs, and Appendix Figure \ref{fig:mad} 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{fig:coverageDGP1DGP4} 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{fig:coverageDG25} and \ref{fig:coverageDG36}, and Appendix Figure \ref{fig:coverageDGP1_DGP4_n200} presents the results for $n=200$.
\begin{figure}[H]
\caption{Bias (multiplied by $n$) before and after correction for DGP1--DGP3}
\begin{center}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEX/QR/opt_feasible.eps}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/triangular_locScaleEX/QR/opt_feasible.eps}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/stable_locScaleEX/QR/opt_feasible.eps}
\end{subfigure}
\end{center}
\footnotesize{\textit{Notes:} The panels display the bias (multiplied by $n$) of the intercept and the slope for classical QR without bias correction (blue dots), QR with feasible bias correction based on the rule-of-thumb bandwidth (gold squares), and QR with infeasible bias correction (gold dashed line) for DGP1--DGP3. All results are based on 5,000 simulation repetitions.}
\label{fig:bias_dgp13}
\end{figure}
\begin{figure}[H]
\caption{Bias (multiplied by $n$) before and after correction for DGP4--DGP6}
\begin{center}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEND/IV/opt_feasible.eps}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/triangular_locScaleEND/IV/opt_feasible.eps}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/stable_locScaleEND/IV/opt_feasible.eps}
\end{subfigure}
\end{center}
\footnotesize{\textit{Notes:} The panels display the bias (multiplied by $n$) of the intercept and the slope for IVQR (implemented via the MILP formulation in Appendix \ref{app:LP_MILP}) without bias correction (blue dots), IVQR with feasible bias correction based on the rule-of-thumb bandwidth (gold squares), and IVQR with infeasible bias correction (gold dashed line) for DGP4--DGP6. All results are based on 5,000 simulation repetitions. The infeasible bias correction is based on the feasible formula applied to a simulated sample of 10,000,000 observations.}
\label{fig:bias_dgp13_iv}
\end{figure}
\begin{figure}[H]
\caption{RMSE comparison of raw and bias-corrected estimators }\label{fig:RMSEuniform}
\begin{center}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEX/QR/opt_rmse.eps}
\caption{}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEND/IV/opt_rmse.eps}
\caption{}
\end{subfigure}
\end{center}
\footnotesize{\textit{Notes:}
The panels compare the RMSE for estimators without bias correction (blue), with infeasible bias correction (grey) and with feasible bias correction based on the rule-of-thumb bandwidth choice (gold) for (a) DGP1, classical QR and (b) DGP4, IVQR. All results are based on 5,000 simulation repetitions.
}
\label{fig:rmse14}
\end{figure}
\begin{figure}[H]
\caption{Confidence interval coverage before and after correction for DGP1 and DGP4}
\begin{center}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEX/QR/opt_coverage.eps}
\caption{}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\includegraphics[width=\textwidth]{sections/output/uniform_locScaleEND/IV/opt_coverage.eps}
\caption{}
\end{subfigure}
\end{center}
\footnotesize{\textit{Notes:}
The panels display the coverage probability of the $90\% $ confidence intervals for the intercept and the slope without bias correction (blue dots) and with the feasible bias correction based on the rule-of-thumb bandwidth choice (gold squares) for DGP1 (classical QR) and DGP4 (IVQR). All results are based on 5,000 simulation repetitions.}
\label{fig:coverageDGP1DGP4}
\end{figure}
\section{Empirical application}
\label{sec: 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 \citet{engel1857}, analyzed by \citet{koenker1982robust} and \citet{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 \texttt{R} package \texttt{quantreg} \citep{quantreg_package}.
One feature of these data is the growing dispersion of the outcome variable (food expenditure) as a function of the regressor (income) \citep{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.
\begin{figure}[H]
\caption{Quantile regression of annual food expenditure on income}
\begin{center}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{sections/figuresRR2/fig4_engel1.eps}
\caption{Impact of bias correction}
\end{subfigure}
\begin{subfigure}[b]{0.45\textwidth}
\includegraphics[width=\textwidth]{sections/figuresRR2/fig4_engel2.eps}
\caption{Composition of second-order bias}
\end{subfigure}
\end{center}
\footnotesize{\textit{Notes:} Panel (a) compares the classical QR estimates (blue dots) to the bias-corrected QR estimates (gold squares) with 90\% confidence intervals (bars). The bias correction is based on the rule-of-thumb bandwidth choice with $(\tilde{A}_G,\tilde{A}_Q,\tilde{A}_\kappa)=(2,1.5,2)$. Panel (b) shows the contributions of the different bias components to the overall second-order bias.
}
\label{fig:empirical_application}
\end{figure}
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{fig:empirical_application} 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:
\begin{align*}
\thetahat-\hat\theta_{bc}= \underbrace{\frac{1}{2}\hat{G}^{-1} \bigBra{\gthat(\thetahat) -\gthat^*(-\thetahat)}}_{\text{(i)}}
-\underbrace{ \frac{1}{n} \hat G^{-1} \hat \kappa_{\tau }}_{\text{(ii)}}
-\underbrace{ \frac{1}{2n}G^{-1}\hat Q^\prime vec( \hat \Omega )}_{\text{(iii)}} .
\end{align*}
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.
\section{Conclusion}
\label{sec: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.
\section*{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.
\bibliographystyle{ecta}
\bibliography{main}
\newpage
\part*{Online appendix}