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.
136,511 characters · 18 sections · 77 citation commands
Smoothed Estimating Equations for Instrumental Variables Quantile Regression
\paragraph{\sc Abstract}
The moment conditions or estimating equations for instrumental variables quantile regression involve the discontinuous indicator function. We instead use smoothed estimating equations (SEE), with bandwidth $h$. We show that the mean squared error (MSE) of the vector of the SEE is minimized for some $h>0$, leading to smaller asymptotic MSE of the estimating equations and associated parameter estimators. The same MSE-optimal $h$ also minimizes the higher-order type I error of a SEE-based $\chi ^{2}$ test and increases size-adjusted power in large samples. Computation of the SEE estimator also becomes simpler and more reliable, especially with (more) endogenous regressors. Monte Carlo simulations demonstrate all of these superior properties in finite samples, and we apply our estimator to JTPA data. Smoothing the estimating equations is not just a technical operation for establishing Edgeworth expansions and bootstrap refinements; it also brings the real benefits of having more precise estimators and more powerful tests.
\allowdisplaybreaks[4]
Many econometric models are specified by moment conditions or estimating equations. An advantage of this approach is that the full distribution of the data does not have to be parameterized. In this paper, we consider estimating equations that are not smooth in the parameter of interest. We focus on instrumental variables quantile regression (IV-QR), which includes the usual quantile regression as a special case. Instead of using the estimating equations that involve the nonsmooth indicator function, we propose to smooth the indicator function, leading to our smoothed estimating equations (SEE) and SEE estimator.
Our SEE estimator has several advantages. First, from a computational point of view, the SEE estimator can be computed using any standard iterative algorithm that requires smoothness. This is especially attractive in IV-QR where simplex methods for the usual QR are not applicable. In fact, the SEE approach has been used in ChenPouzo2009,ChenPouzo2012 for computing their nonparametric sieve estimators in the presence of nonsmooth moments or generalized residuals. However, a rigorous investigation is currently lacking. Our paper can be regarded as a first step towards justifying the SEE approach in nonparametric settings. Relatedly, FanLiao2014 have employed the same strategy of smoothing the indicator function to reduce the computational burden of their focused GMM approach. Second, from a technical point of view, smoothing the estimating equations enables us to establish high-order properties of the estimator. This motivates Horowitz1998, for instance, to examine a smoothed objective function for median regression, to show high-order bootstrap refinement. Instead of smoothing the objective function, we show that there is an advantage to smoothing the estimating equations. This point has not been recognized and emphasized in the literature. For QR estimation and inference via empirical likelihood, Otsu2008 and Whang2006 also examine smoothed estimators. To the best of our knowledge, no paper has examined smoothing the estimating equations for the usual QR estimator, let alone IV-QR. Third, from a statistical point of view, the SEE estimator is a flexible class of estimators that includes the IV/OLS mean regression estimators and median and quantile regression estimators as special cases. Depending on the smoothing parameter, the SEE estimator can have different degrees of robustness in the sense of Huber1964. By selecting the smoothing parameter appropriately, we can harness the advantages of both the mean regression estimator and the median/quantile regression estimator. Fourth and most importantly, from an econometric point of view, smoothing can reduce the mean squared error (MSE) of the SEE, which in turn leads to a smaller asymptotic MSE of the parameter estimator and to more powerful tests. We seem to be the first to establish these advantages.
In addition to investigating the asymptotic properties of the SEE estimator, we provide a smoothing parameter choice that minimizes different criteria: the MSE of the SEE, the type I error of a chi-square test subject to exact asymptotic size control, and the approximate MSE of the parameter estimator. We show that the first two criteria produce the same optimal smoothing parameter, which is also optimal under a variant of the third criterion. With the data-driven smoothing parameter choice, we show that the statistical and econometric advantages of the SEE estimator are reflected clearly in our simulation results.
There is a growing literature on IV-QR. For a recent review, see ChernozhukovHansen2013. Our paper is built upon ChernozhukovHansen2005, which establishes a structural framework for IV-QR and provides primitive identification conditions. Within this framework, ChernozhukovHansen2006 and ChernozhukovEtAl2009 develop estimation and inference procedures under strong identification. For inference procedures that are robust to weak identification, see ChernozhukovHansen2008 and Jun2008, for example. IV-QR can also reduce bias for dynamic panel fixed effects estimation as in Galvao2011. None of these papers considers smoothing the IV-QR estimating equations; that idea (along with minimal first-order theory) seemingly first appeared in an unpublished draft by MaCurdyHong1999, although the idea of smoothing the indicator function in general appears even earlier, as in Horowitz1992 for the smoothed maximum score estimator. An alternative approach to overcome the computational obstacles in the presence of a nonsmooth objective function is to explore the asymptotic equivalence of the Bayesian and classical methods for regular models and use the MCMC approach to obtain the classical extremum estimator; see ChernozhukovHong2003, whose Example 3 is IV-QR. As a complement, our approach deals with the computation problem in the classical framework directly.
The rest of the paper is organized as follows. Section (ref) describes our setup and discusses some illuminating connections with other estimators. Sections (ref), (ref), and (ref) calculate the MSE of the SEE, the type I and type II errors of a chi-square test, and the approximate MSE of the parameter estimator, respectively. Section (ref) applies our estimator to JTPA data, and Section (ref) presents simulation results before we conclude. Longer proofs and calculations are gathered in the appendix.
We are interested in estimating the instrumental variables quantile regression (IV-QR) model
where $\mathbb E\mathopen{}\mathclose \bgroup \originalleft[Z_{j}\mathopen{}\mathclose \bgroup \originalleft( 1\{U_{j}<0\}-q\aftergroup \egroup \originalright)\aftergroup \egroup \originalright] =0$ for instrument vector $ Z_{j}\in \mathbb{R}^{d}$ and $1\{ \cdot \}$ is the indicator function. Instruments are taken as given; this does not preclude first determining the efficient set of instruments as in Newey2004 or NeweyPowell1990, for example. We restrict attention to the “just identified” case $X_{j}\in \mathbb{R} ^{d}$ and iid data for simpler exposition; for the overidentified case, see (ref) below.
A special case of this model is exogenous QR with $Z_{j}=X_{j}$, which is typically estimated by minimizing a criterion function:
where $\rho _{q}(u)\equiv \mathopen{}\mathclose \bgroup \originalleft( q-1\{u<0\}\aftergroup \egroup \originalright) u$ is the check function. Since the objective function is not smooth, it is not easy to obtain a high-order approximation to the sampling distribution of $\hat{\beta}_{Q}$. To avoid this technical difficulty, Horowitz1998 proposes to smooth the objective function to obtain
where $G(\cdot )$ is a smooth function and $h$ is the smoothing parameter or bandwidth. Instead of smoothing the objective function, we smooth the underlying moment condition and define $\hat{\beta}$ to be the solution of the vector of smoothed estimating equations (SEE) $m_{n}(\hat{\beta})=0$, where\footnote{ It suffices to have $m_{n}(\hat{\beta})=o_{p}(1)$, which allows for a small error when $\hat{\beta}$ is not the exact solution to $m_{n}(\hat{\beta})=0$. }
Our approach is related to kernel-based nonparametric conditional quantile estimators. The moment condition there is $\mathbb E\mathopen{}\mathclose \bgroup \originalleft[ 1\{X=x\} \mathopen{}\mathclose \bgroup \originalleft( 1\{Y<\beta \}-q\aftergroup \egroup \originalright) \aftergroup \egroup \originalright] =0$. Usually the $1\{X=x\}$ indicator function is “smoothed” with a kernel, while the latter term is not. This yields the nonparametric conditional quantile estimator $\hat{\beta}_{q}(x)= \mathop{\rm arg\,min}_{b}\sum_{i=1}^{n}\rho _{q}(Y_{i}-b)K[(x-X_{i})/h]$ for the conditional $q$-quantile at $X=x$, estimated with kernel $K(\cdot)$ and bandwidth $h$. Our approach is different in that we smooth the indicator $ 1\{Y<\beta \}$ rather than $1\{X=x\}$. Smoothing both terms may help but is beyond the scope of this paper.
Estimating $\hat{\beta}$ from the SEE is computationally easy: $d$ equations for $d$ parameters, and a known, analytic Jacobian. Computationally, solving our problem is faster and more reliable than the IV-QR method in ChernozhukovHansen2006, which requires specification of a grid of endogenous coefficient values to search over, computing a conventional QR estimator for each grid point. This advantage is important particularly when there are more endogenous variables.
If the model is overidentified with $\dim (Z_{j})>\dim (X_{j})$, we can use a $\dim (X_{j})\times \dim (Z_{j})$ matrix $\mathbb{W}$ to transform the original moment conditions $\mathbb E\mathopen{}\mathclose \bgroup \originalleft[Z_{j}\mathopen{}\mathclose \bgroup \originalleft(q-1\mathopen{}\mathclose \bgroup \originalleft\{Y_{j}<X_{j}^{\prime }\beta \aftergroup \egroup \originalright\}\aftergroup \egroup \originalright)\aftergroup \egroup \originalright]=0$ into
Then we have an exactly identified model with transformed instrument vector $ \tilde{Z}_{j}$, and our asymptotic analysis can be applied to (ref).
By the theory of optimal estimating equations or efficient two-step GMM, the optimal $\mathbb{W}$ takes the following form:
where $f_{U|Z,X}(0)$ is the conditional PDF of $U$ evaluated at $U=0$ given $ \mathopen{}\mathclose \bgroup \originalleft( Z,X\aftergroup \egroup \originalright) $ and $\sigma ^{2}\mathopen{}\mathclose \bgroup \originalleft( Z\aftergroup \egroup \originalright) =\mathrm{Var} \mathopen{}\mathclose \bgroup \originalleft(1\mathopen{}\mathclose \bgroup \originalleft \{ U<0\aftergroup \egroup \originalright \} \mid Z\aftergroup \egroup \originalright)$. The standard two-step approach requires an initial estimator of $\beta _{0}$ and nonparametric estimators of $f_{U|Z,X}(0)$ and $\sigma ^{2}\mathopen{}\mathclose \bgroup \originalleft( Z\aftergroup \egroup \originalright)$. The underlying nonparametric estimation error may outweigh the benefit of having an optimal weighting matrix. This is especially a concern when the dimensions of $X$ and $Z$ are large. The problem is similar to what HwangSun2015 consider in a time series GMM framework where the optimal weighting matrix is estimated using a nonparametric HAC approach. Under the alternative and more accurate asymptotics that captures the estimation error of the weighting matrix, they show that the conventionally optimal two-step approach does not necessarily outperform a first-step approach that does not employ a nonparametric weighting matrix estimator. While we expect a similar qualitative message here, we leave a rigorous analysis to future research.
In practice, a simple procedure is to ignore $ f_{U|Z,X}(0) $ and $\sigma ^{2}\mathopen{}\mathclose \bgroup \originalleft( Z\aftergroup \egroup \originalright) $ (or assume that they are constants) and employ the following empirical weighting matrix,
This choice of $\mathbb{W}_{n}$ is in the spirit of the influential work of LiangZeger1986 who advocate the use of a working correlation matrix in constructing the weighting matrix. Given the above choice of $\mathbb{W} _{n}$, $\tilde{Z}_{j}$ is the least squares projection of $X_{j}$ on $Z_{j}$ . It is easy to show that with some notational changes our asymptotic results remain valid in this case.
An example of an overidentified model is the conditional moment model
In this case, any measurable function of $Z_{j}$ can be used as an instrument. As a result, the model could be overidentified. According to Chamberlain1987 and Newey1990, the optimal set of instruments in our setting is given by
Let $F_{U|Z,X}\mathopen{}\mathclose \bgroup \originalleft( u\mid z,x\aftergroup \egroup \originalright) $ and $ f_{U|Z,X}\mathopen{}\mathclose \bgroup \originalleft( u\mid z,x\aftergroup \egroup \originalright) $ be the conditional distribution function and density function of $U$ given $\mathopen{}\mathclose \bgroup \originalleft( Z,X\aftergroup \egroup \originalright) =(z,x)$. Then under some regularity conditions,
The optimal instruments involve the conditional density $f_{U|Z,X}\mathopen{}\mathclose \bgroup \originalleft( u\mid z,x\aftergroup \egroup \originalright) $ and a conditional expectation. In principle, these objects can be estimated nonparametrically. However, the nonparametric estimation uncertainty can be very high, adversely affecting the reliability of inference. A simple and practical strategy\footnote{We are not alone in recommending this simple strategy for empirical work. ChernozhukovHansen2006 make the same recommendation in their Remark 5 and use this strategy in their empirical application. See also Kwak2010.} is to construct the optimal instruments as the OLS projection of each $ X_{j}$ onto some sieve basis functions $\Phi ^{K}\mathopen{}\mathclose \bgroup \originalleft( Z_{j}\aftergroup \egroup \originalright) \equiv \mathopen{}\mathclose \bgroup \originalleft[ \Phi _{1}(Z_{j}),\ldots,\Phi _{K}(Z_{j})\aftergroup \egroup \originalright] ^{\prime }$, leading to
as the instruments. Here $\mathopen{}\mathclose \bgroup \originalleft \{ \Phi _{i}\mathopen{}\mathclose \bgroup \originalleft( \cdot \aftergroup \egroup \originalright) \aftergroup \egroup \originalright \} $ are the basis functions such as power functions. Since the dimension of $\tilde{Z}_{j}$ is the same as the dimension of $ X_{j}$, our asymptotic analysis can be applied for any fixed value of $K$.\footnote{A theoretically efficient estimator can be obtained using the sieve minimum distance approach. It entails first estimating the conditional expectation $\mathbb E \mathopen{}\mathclose \bgroup \originalleft[ \mathopen{}\mathclose \bgroup \originalleft( 1\{Y_{j}<X_{j}\beta \}-q\aftergroup \egroup \originalright) \mid Z_{j}\aftergroup \egroup \originalright] $ using $\Phi ^{K}\mathopen{}\mathclose \bgroup \originalleft( Z_{j}\aftergroup \egroup \originalright) $ as the basis functions and then choosing $\beta $ to minimize a weighted sum of squared conditional expectations. See, for example, ChenPouzo2009,ChenPouzo2012. To achieve the semiparametric efficiency bound, $K$ has to grow with the sample size at an appropriate rate. In work in progress, we consider nonparametric quantile regression with endogeneity and allow $K$ to diverge, which is necessary for both identification and efficiency. Here we are content with a fixed $K$ for empirical convenience at the cost of possible efficiency loss.}
For the special case $Z_{j}=X_{j}$, we compare the SEE with the estimating equations derived from smoothing the criterion function as in Horowitz1998. The first order condition of the smoothed criterion function, evaluated at the true $ \beta _{0}$, is
The first term agrees with our proposed SEE. Technically, it should be easier to establish high-order results for our SEE estimator since it has one fewer term. Later we show that the absolute bias of our SEE estimator is smaller, too. Another subtle point is that our SEE requires only the estimating equation $\mathbb E\mathopen{}\mathclose \bgroup \originalleft[X_{j}\mathopen{}\mathclose \bgroup \originalleft( 1\{U_{j}<0\}-q\aftergroup \egroup \originalright)\aftergroup \egroup \originalright] =0$, whereas Horowitz1998 has to impose an additional condition to ensure that the second term in the FOC is approximately mean zero.
When $h\rightarrow \infty $, $G(\cdot )$ only takes arguments near zero and thus can be approximated well linearly. For example, with the $G(\cdot )$ from Whang2006 and Horowitz1998, $G(v)=0.5+(105/64)v+O(v^{3}) $ as $v\rightarrow 0$. Ignoring the $O(v^{3})$, the corresponding estimator $ \hat{\beta}_{\infty }$ is defined by
where $e_{1}=(1,0,\ldots ,0)^{\prime }$ is $d\times 1$, $\mathbf{1} _{n,1}=(1,1,\ldots ,1)^{\prime }$ is $n\times 1$, $X$ and $Z$ are $n\times d$ with respective rows $X_{i}^{\prime }$ and $Z_{i}^{\prime }$, and using the fact that the first column of $X$ is $\mathbf{1}_{n,1}$ so that $Xe_{1}= \mathbf{1}_{n,1}$. It then follows that
As $h$ grows large, the smoothed QR estimator approaches the IV estimator plus an adjustment to the intercept term that depends on $q$, the bandwidth, and the slope of $G(\cdot )$ at zero. In the special case $Z_{j}=X_{j}$, the IV estimator is the OLS estimator.\footnote{ This is different from ZhouEtAl2011, who add the $d$ OLS moment conditions to the $d$ median regression moment conditions before estimation; our connection to IV/OLS emerges naturally from smoothing the (IV)QR estimating equations.}
The intercept is often not of interest, and when $q=0.5$, the adjustment is zero anyway. The class of SEE estimators is a continuum (indexed by $h$) with two well-known special cases at the extremes: unsmoothed IV-QR and mean IV. For $q=0.5$ and $Z_j=X_j$, this is median regression and mean regression (OLS). Well known are the relative efficiency advantages of the median and the mean for different error distributions. Our estimator with a data-driven bandwidth can harness the advantages of both, without requiring the practitioner to make guesses about the unknown error distribution.
With $Z_j=X_j$, the result that our SEE can yield OLS when $h\to \infty$ or median regression when $h=0$ calls to mind robust estimators like the trimmed or Winsorized mean (and corresponding regression estimators). Setting the trimming/Winsorization parameter to zero generates the mean while the other extreme generates the median. However, our SEE mechanism is different and more general/flexible; trimming/Winsorization is not directly applicable to $q\ne0.5$; our method to select the smoothing parameter is novel; and the motivations for QR extend beyond (though include) robustness.
With $X_{i}=1$ and $q=0.5$ (population median estimation), our SEE becomes
If $G'(u)=1\{-1\le u\le 1\}/2$ (the uniform kernel), then $H(u)\equiv 2G(u)-1=u$ for $u\in[-1,1]$, $H(u)=1$ for $u>1$, and $H(u)=-1$ for $u<-1$. The SEE is then $0=\sum_{i=1}^{n}\psi\mathopen{}\mathclose \bgroup \originalleft(Y_i;\beta\aftergroup \egroup \originalright)$ with $\psi\mathopen{}\mathclose \bgroup \originalleft(Y_i;\beta\aftergroup \egroup \originalright)=H\mathopen{}\mathclose \bgroup \originalleft(\mathopen{}\mathclose \bgroup \originalleft(\beta-Y_i\aftergroup \egroup \originalright)/h\aftergroup \egroup \originalright)$. This produces the Winsorized mean estimator of the type in Huber1964.\footnote{ For a strict mapping, multiply by $h$ to get $\psi (Y_{i};\beta )=hH[(\beta -Y_{i})/h]$. The solution is equivalent since $\sum h\psi (Y_{i};\beta )=0$ is the same as $\sum \psi (Y_{i};\beta )=0$ for any nonzero constant $h$.}
Further theoretical comparison of our SEE-QR with trimmed/Winsorized mean regression (and the IV versions) would be interesting but is beyond the scope of this paper. For more on robust location and regression estimators, see for example Huber1964, KoenkerBassett1978, and RuppertCarroll1980.
Since statistical inference can be made based on the estimating equations (EEs), we examine the mean squared error (MSE) of the SEE. An advantage of using EEs directly is that inference can be made robust to the strength of identification. Our focus on the EEs is also in the same spirit of the large literature on optimal estimating equations. For the historical developments of EEs and their applications in econometrics, see BeraEtAl2006. The MSE of the SEE is also related to the estimator MSE and inference properties both intuitively and (as we will show) theoretically. Such results may provide helpful guidance in contexts where the SEE MSE is easier to compute than the estimator MSE, and it provides insight into how smoothing works in the QR model as well as results that will be used in subsequent sections.
We maintain different subsets of the following assumptions for different results. We write $f_{U|Z}(\cdot \mid z)$ and $F_{U|Z}(\cdot \mid z)$ as the conditional PDF and CDF of $U$ given $Z=z$. We define $f_{U|Z,X}(\cdot \mid z,x)$ and $F_{U|Z,X}(\cdot \mid z,x) $ similarly.
Assumption (ref) describes the sampling process. Assumption (ref) is analogous to Assumption 3 in both Horowitz1998 and Whang2006. As discussed in these two papers, the boundedness assumption for $Z_{j}$, which is a technical condition, is made only for convenience and can be dropped at the cost of more complicated proofs.
Assumption (ref)(i) allows us to use the law of iterated expectations to simplify the asymptotic variance. Our qualitative conclusions do not rely on this assumption. Assumption (ref)(ii) is critical. If we are not willing to make such an assumption, then smoothing will be of no benefit. Inversely, with some small degree of smoothness of the conditional error density, smoothing can leverage this into the advantages described here. Also note that Horowitz1998 assumes $r\geq 4$, which is sufficient for the estimator MSE result in Section (ref).
Assumptions (ref)(i--ii) are analogous to the standard high-order kernel conditions in the kernel smoothing literature. The integral condition in (i) ensures that smoothing reduces (rather than increases) variance. Note that
using the evenness of $G'(u)$. When $r=2$, we can use any $G(u)$ such that $G'(u)$ is a symmetric PDF on $[-1,1]$. In this case, $1-\int_{-1}^1 G^2(u)du>0$ holds automatically. When $r>2$, $G'(u)<0$ for some $u$, and $G(u)$ is not monotonic. It is not easy to sign $1-\int_{-1}^1 G^2(u)du$ generally, but it is simple to calculate this quantity for any chosen $G(\cdot)$. For example, consider $r=4$ and the $G(\cdot)$ function in Horowitz1998 and Whang2006 shown in Figure (ref):
The range of the function is outside $[0,1]$. Simple calculations show that $1-\int_{-1}^1 G^2(u)du>0$.
Assumption (ref)(iii) is needed for the Edgeworth expansion. As Horowitz1998 and Whang2006 discuss, Assumption (ref)(iii) is a technical assumption that (along with Assumption (ref)) leads to a form of Cram\'{e }r's condition, which is needed to justify the Edgeworth expansion used in Section (ref). Any $G(u)$ constructed by integrating polynomial kernels in Muller1984 satisfies Assumption (ref)(iii). In fact, $G(u)$ in (ref) is obtained by integrating a fourth-order kernel given in Table 1 of Muller1984. Assumption (ref) ensures that the bias of the SEE is of smaller order than its variance. It is needed for the asymptotic normality of the SEE as well as the Edgeworth expansion.
Assumption (ref) is an identification assumption. See Theorem 2 of ChernozhukovHansen2006 for more primitive conditions. It ensures the consistency of the SEE estimator. Assumption (ref) is necessary for the $\sqrt{n}$-consistency and asymptotic normality of the SEE estimator.
Define
and abbreviate $m_{n}\equiv m_{n}(\beta _{0})=n^{-1/2}\sum_{j=1}^{n}W_{j}$. The theorem below gives the first two moments of $W_{j}$ and the first-order asymptotic distribution of $m_{n}$.
Compared with the EE derived from smoothing the criterion function as in Horowitz1998, our SEE has smaller bias and variance, and these differences affect the bias and variance of the parameter estimator. The former approach only applies to exogenous QR with $Z_j=X_j$. The EE derived from smoothing the criterion function in (ref) for $Z_j=X_j$ can be written
Consequently, as calculated in the appendix,
The dominating term of the bias of our SEE in (ref) is $r+1$ times smaller in absolute value than that of the EE derived from a smoothed criterion function in (ref). A larger bias can lead to less accurate confidence regions if the same variance estimator is used. Additionally, the smoothed criterion function analog of $\mathbb E(W_jW_j')$ in (ref) has a positive $O(h)$ term instead of the negative $O(h)$ term for SEE. The connection between these terms and the estimator's asymptotic mean squared error (AMSE) is shown in Section (ref) to rely on the inverse of the matrix in equation (ref). Here, though, the sign of the $O(h)$ term is indeterminant since it depends on a PDF derivative. (A negative $O(h)$ term implies higher AMSE since this matrix is inverted in the AMSE expression, and positive implies lower.) If $U=0$ is a mode of the conditional (on $X$) distribution, then the $O(h)$ term is zero and the AMSE comparison is driven by $\mathbb E(W_j)$ and $\mathbb E(W_jW_j')$. Since SEE yields smaller $\mathbb E(W_jW_j')$ and smaller absolute $\mathbb E(W_j)$, it will have smaller estimator AMSE in such cases. Simulation results in Section (ref) add evidence that the SEE estimator usually has smaller MSE in practice.
The first-order asymptotic variance $V$ is the same as the asymptotic variance of
the scaled EE of the unsmoothed IV-QR. The effect of smoothing to reduce variance is captured by the term of order $h$, where $1- \int_{-1}^{1}G^2(u)du>0$ by Assumption (ref)(i). This reduction in variance is not surprising. Replacing the discontinuous indicator function $ 1\{U<0\}$ by a smooth function $G(-U/h)$ pushes the dichotomous values of zero and one into some values in between, leading to a smaller variance. The idea is similar to Breiman1994's (Breiman1994) bagging (bootstrap aggregating), among others.
Define the MSE of the SEE to be $\mathbb E\mathopen{}\mathclose \bgroup \originalleft(m_{n}^{\prime }V^{-1}m_{n}\aftergroup \egroup \originalright)$. Building upon (ref) and (ref), and using $W_{i} \mathpalette{\protect \independenT}{\perp}W_{j}$ for $i\neq j$, we have:
where
Ignoring the $o(\cdot )$ term, we obtain the asymptotic MSE of the SEE. We select the smoothing parameter to minimize the asymptotic MSE:
The proposition below gives the optimal smoothing parameter $h_{\text{SEE} }^{\ast }$.
When $r=2$, the MSE-optimal $h_{\text{SEE}}^{\ast }\asymp n^{-1/(2r-1)}=n^{-1/3}$. This is smaller than $n^{-1/5}$, the rate that minimizes the MSE of estimated standard errors of the usual regression quantiles. Since nonparametric estimators of $f_{U}^{(r-1)}(0)$ converge slowly, we propose a parametric plug-in described in Section (ref).
We point out in passing that the optimal smoothing parameter $h_{\text{SEE} }^{\ast }$ is invariant to rotation and translation of the (non-constant) regressors. This may not be obvious but can be proved easily.
For the unsmoothed IV-QR, let
then the MSE of the estimating equations is $\mathbb E\mathopen{}\mathclose \bgroup \originalleft(\tilde{m}_{n}^{\prime }V^{-1} \tilde{m}_{n}\aftergroup \egroup \originalright)=d$. Comparing this to the MSE of the SEE given in (ref), we find that the SEE has a smaller MSE when $h=h_{\text{SEE} }^{\ast }$ because
In terms of MSE, it is advantageous to smooth the estimating equations. To the best of our knowledge, this point has never been discussed before in the literature.
In this section, we explore the effect of smoothing on a chi-square test. Other alternatives for inference exist, such as the Bernoulli-based MCMC-computed method from ChernozhukovEtAl2009, empirical likelihood as in Whang2006, and bootstrap as in Horowitz1998, where the latter two also use smoothing. Intuitively, when we minimize the MSE, we may expect lower type I error: the $\chi ^{2}$ critical value is from the unsmoothed distribution, and smoothing to minimize MSE makes large values (that cause the test to reject) less likely. The reduced MSE also makes it easier to distinguish the null hypothesis from some given alternative. This combination leads to improved size-adjusted power. As seen in our simulations, this is true especially for the IV case.
Using the results in Section (ref) and under Assumption (ref), we have
where we continue to use the notation $m_{n}\equiv m_{n}(\beta _{0})$. From this asymptotic result, we can construct a hypothesis test that rejects the null hypothesis $H_{0}:\beta =\beta _{0}$ when
where
is a consistent estimator of $V$ and $c_{\alpha }\equiv \chi _{d,1-\alpha }^{2}$ is the $1-\alpha $ quantile of the chi-square distribution with $d$ degrees of freedom. As desired, the asymptotic size is
Here $P\equiv P_{\beta _{0}}$ is the probability measure under the true model parameter $\beta _{0}$. We suppress the subscript $\beta _{0}$ when there is no confusion.
It is important to point out that the above result does not rely on the strong identification of $\beta _{0}$. It still holds if $\beta _{0}$ is weakly identified or even unidentified. This is an advantage of focusing on the estimating equations instead of the parameter estimator. When a direct inference method based on the asymptotic normality of $\hat{\beta}$ is used, we have to impose Assumptions (ref) and (ref).
To more precisely measure the type I error $P\mathopen{}\mathclose \bgroup \originalleft( S_{n}>c_{\alpha }\aftergroup \egroup \originalright) $, we first develop a high-order stochastic expansion of $S_{n}$. Let $ V_{n}\equiv \mathrm{Var}\mathopen{}\mathclose \bgroup \originalleft( m_{n}\aftergroup \egroup \originalright) $. Following the same calculation as in (ref), we have
where $V^{1/2}$ is the matrix square root of $V$ such that $V^{1/2}\mathopen{}\mathclose \bgroup \originalleft( V^{1/2}\aftergroup \egroup \originalright) ^{\prime }=V$. We can choose $V^{1/2}$ to be symmetric but do not have to.
Details of the following are in the appendix; here we outline our strategy and highlight key results. Letting
such that $\Lambda _{n}\Lambda _{n}^{\prime }=V_{n}$, and defining
we can approximate the test statistic as $S_{n}=S_{n}^{L}+e_{n}$, where
and $e_{n}$ is the remainder term satisfying $P\mathopen{}\mathclose \bgroup \originalleft( \mathopen{}\mathclose \bgroup \originalleft \vert e_{n}\aftergroup \egroup \originalright \vert >O\mathopen{}\mathclose \bgroup \originalleft( h^{2}\aftergroup \egroup \originalright) \aftergroup \egroup \originalright) =O\mathopen{}\mathclose \bgroup \originalleft( h^{2}\aftergroup \egroup \originalright) $.
The stochastic expansion above allows us to approximate the characteristic function of $S_{n}$ with that of $S_{n}^{L}$. Taking the Fourier--Stieltjes inverse of the characteristic function yields an approximation of the distribution function, from which we can calculate the type I error by plugging in the critical value $c_{\alpha}$.
From Theorem (ref), an approximate measure of the type I error of the SEE-based chi-square test is
and an approximate measure of the coverage probability error (CPE) is \footnote{ The CPE is defined to be the nominal coverage minus the true coverage probability, which may be different from the usual definition. Under this definition, smaller CPE corresponds to higher coverage probability (and smaller type I error).}
which is also the error in rejection probability under the null.
Up to smaller-order terms, the term $nh^{2r}\mathbb E(B)'\mathbb E(B)$ characterizes the bias effect from smoothing. The bias increases type I error and reduces coverage probability. The term $h\mathrm{tr}\mathopen{}\mathclose \bgroup \originalleft[\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright) \aftergroup \egroup \originalright]$ characterizes the variance effect from smoothing. The variance reduction decreases type I error and increases coverage probability. The type I error is $\alpha $ up to order $O\mathopen{}\mathclose \bgroup \originalleft(h+nh^{2r}\aftergroup \egroup \originalright)$. There exists some $h>0$ that makes bias and variance effects cancel, leaving type I error equal to $\alpha$ up to smaller-order terms in $R_{n}$.
Note that $nh^{2r}\mathbb E(B)'\mathbb E(B) -h\mathrm{tr}\mathopen{}\mathclose \bgroup \originalleft[\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright) \aftergroup \egroup \originalright]$ is identical to the high-order term in the asymptotic MSE of the SEE in (ref). The $h_{ \text{CPE}}^{\ast }$ that minimizes type I error is the same as $h_{\text{SEE }}^{\ast }$.
The result that $h_{\text{CPE}}^{\ast }=h_{\text{SEE}}^{\ast }$ is intuitive. Since $h_{\text{SEE}}^{\ast }$ minimizes $\mathbb E\mathopen{}\mathclose \bgroup \originalleft(m_{n}^{\prime }V^{-1}m_{n}\aftergroup \egroup \originalright)$, for a test with $c_{\alpha }$ and $\hat{V}$ both invariant to $h$, the null rejection probability $P\mathopen{}\mathclose \bgroup \originalleft(m_{n}^{\prime }\hat{V} ^{-1}m_{n}>c_{\alpha }\aftergroup \egroup \originalright)$ should be smaller when the SEE's MSE is smaller.
When $h=h_{\text{CPE}}^{\ast }$,
where $C^{+}=\mathopen{}\mathclose \bgroup \originalleft( 1-\frac{1}{2r}\aftergroup \egroup \originalright) \mathrm{tr}\mathopen{}\mathclose \bgroup \originalleft[\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright) \aftergroup \egroup \originalright] >0$. If instead we construct the test statistic based on the unsmoothed estimating equations, $\tilde{S}_{n}=\tilde{m} _{n}^{\prime }\hat{V}^{-1}\tilde{m}_{n}$, then it can be shown that
for some constant $C$, which is in general not equal to zero. Given that $ n^{-1/2}=o(h_{\text{CPE}}^{\ast })$ and $C^{+}>0$, we can expect the SEE-based chi-square test to have a smaller type I error in large samples.
To obtain the local asymptotic power of the $S_{n}$ test, we let the true parameter value be $\beta _{n}=\beta _{0}-\delta /\sqrt{n}$, where $\beta _{0} $ is the parameter value that satisfies the null hypothesis $H_{0}$. In this case,
In the proof of Theorem (ref), we show that
When $\delta=0$, which implies $\tau=0$, the expansion in Theorem (ref) reduces to that in Theorem (ref).
When $h=h_{\text{SEE}}^{\ast }$, it follows from Theorem (ref) that
To remove the error in rejection probability of order $h_{\text{SEE}}^{\ast } $, we make a correction to the critical value $c_{\alpha }$. Let $ c_{\alpha }^{\ast }$ be a high-order corrected critical value such that $ P_{\beta_{0}}\mathopen{}\mathclose \bgroup \originalleft( S_{n}>c_{\alpha }^{\ast }\aftergroup \egroup \originalright) =\alpha +o(h_{\text{SEE} }^{\ast})$. Simple calculation shows that
meets the requirement.
To approximate the size-adjusted power of the $S_{n}$ test, we use $ c_{\alpha }^{\ast }$ rather than $c_{\alpha }$ because $c_{\alpha }^{\ast }$ leads to a more accurate test in large samples. Using Theorem (ref) , we can prove the following corollary.
In the asymptotic expansion of the local power function in (ref), $1-\mathcal{G}_{d}\mathopen{}\mathclose \bgroup \originalleft( c_{\alpha };\tau ^{2}\aftergroup \egroup \originalright) $ is the usual first-order power of a standard chi-square test. The next term of order $O(h_{\text{SEE}}^{\ast })$ captures the effect of smoothing the estimating equations. To sign this effect, we plot the function $Q_{d}\mathopen{}\mathclose \bgroup \originalleft( c_{\alpha },\tau ^{2},r\aftergroup \egroup \originalright) $ against $\tau ^{2}$ for $r=2$, $\alpha =10\%$, and different values of $d$ in Figure (ref). Figures for other values of $r$ and $\alpha $ are qualitatively similar. The range of $\tau ^{2}$ considered in Figure (ref) is relevant as the first-order local asymptotic power, i.e.,\ $1-\mathcal{G} _{d}\mathopen{}\mathclose \bgroup \originalleft( c_{\alpha };\tau ^{2}\aftergroup \egroup \originalright) $, increases from $10\%$ to about $ 94\%$, $96\%$, $97\%$, and $99\%$, respectively for $d=1,2,3,4$. It is clear from this figure that $Q_{d}\mathopen{}\mathclose \bgroup \originalleft( c_{\alpha },\tau ^{2},r\aftergroup \egroup \originalright) >0$ for any $\tau ^{2}>0$. This indicates that smoothing leads to a test with improved power. The power improvement increases with $r$. The smoother the conditional PDF of $U$ in a neighborhood of the origin is, the larger the power improvement is.
In this section, we examine the approximate MSE of the parameter estimator. The approximate MSE, being a Nagar-type approximation Nagar1959, can be motivated from the theory of optimal estimating equations, as presented in Heyde1997, for example.
The SEE estimator $\hat{\beta}$ satisfies $m_{n}(\hat{\beta})=0$. In Lemma (ref) in the appendix, we show that
and
Consequently, the approximate MSE (AMSE) of $\sqrt{n}\mathopen{}\mathclose \bgroup \originalleft( \hat{\beta} -\beta _{0}\aftergroup \egroup \originalright) $ is\footnote{ Here we follow a common practice in the estimation of nonparametric and nonlinear models and define the AMSE to be the MSE of $\sqrt{n}\mathopen{}\mathclose \bgroup \originalleft( \hat{ \beta}-\beta _{0}\aftergroup \egroup \originalright) $ after dropping some smaller-order terms. So the asymptotic MSE we define here is a Nagar-type approximate MSE. See Nagar1959.}
where
The first term of $\mathrm{AMSE}_{\beta }$ is the asymptotic variance of the unsmoothed QR estimator. The second term captures the higher-order effect of smoothing on the AMSE of $\sqrt{n}(\hat{\beta}-\beta _{0})$. When $ nh^{r}\rightarrow \infty $ and $n^{3}h^{4r+1}\rightarrow \infty$, we have $ h^{r}=o\mathopen{}\mathclose \bgroup \originalleft( nh^{2r}\aftergroup \egroup \originalright) $ and $1/\sqrt{nh}=o\mathopen{}\mathclose \bgroup \originalleft( nh^{2r}\aftergroup \egroup \originalright) $, so the terms of order $O_{p}(1/\sqrt{nh})$ in (ref) and of order $O\mathopen{}\mathclose \bgroup \originalleft( h^{r}\aftergroup \egroup \originalright) $ in (ref) are of smaller order than the $O(nh^{2r})$ and $O(h)$ terms in the AMSE. If $h\asymp n^{-1/(2r-1)} $ as before, these rate conditions are satisfied when $r>2$.
The optimal $h^{\ast }$ that minimizes the high-order AMSE satisfies
in the sense that the difference between the two sides is nonpositive definite for all $h$. This is equivalent to
This choice of $h$ can also be motivated from the theory of optimal estimating equations. Given the estimating equations $m_{n}=0$, we follow Heyde1997 and define the standardized version of $m_{n}$ by
We include $h$ as an argument of $m_{n}^{s}$ to emphasize the dependence of $ m_{n}^{s}$ on $h$. The standardization can be motivated from the following considerations. On one hand, the estimating equations need to be close to zero when evaluated at the true parameter value. Thus we want $ \mathbb E(m_{n}m_{n}^{\prime })$ to be as small as possible. On the other hand, we want $m_{n}\mathopen{}\mathclose \bgroup \originalleft( \beta +\delta \beta \aftergroup \egroup \originalright) $ to differ as much as possible from $m_{n}\mathopen{}\mathclose \bgroup \originalleft( \beta \aftergroup \egroup \originalright) $ when $\beta $ is the true value. That is, we want $\mathbb E\frac{\partial }{\partial \beta ^{\prime }}m_{n}\mathopen{}\mathclose \bgroup \originalleft( \beta _{0}\aftergroup \egroup \originalright) $ to be as large as possible. To meet these requirements, we choose $h$ to maximize
More specifically, $h^{\ast }$ is optimal if
is nonnegative definite for all $h\in \mathbb{R}^{+}$. But $ \mathbb E\mathopen{}\mathclose \bgroup \originalleft[m_{n}^{s}\mathopen{}\mathclose \bgroup \originalleft( m_{n}^{s}\aftergroup \egroup \originalright) ^{\prime }\aftergroup \egroup \originalright]=\mathopen{}\mathclose \bgroup \originalleft( \mathrm{AMSE}_{\beta }\aftergroup \egroup \originalright) ^{-1}$, so maximizing $\mathbb E\mathopen{}\mathclose \bgroup \originalleft[m_{n}^{s}\mathopen{}\mathclose \bgroup \originalleft( m_{n}^{s}\aftergroup \egroup \originalright) ^{\prime }\aftergroup \egroup \originalright]$ is equivalent to minimizing $\mathrm{AMSE}_{\beta }$.
The question is whether such an optimal $h$ exists. If it does, then the optimal $h^{\ast }$ satisfies
for all $u\in \mathbb{R}^{d}$, by the definition of nonpositive definite plus the fact that the above yields a unique minimizer for any $u$. Using unit vectors $e_1=(1,0,\ldots,0)$, $e_2=(0,1,0,\ldots,0)$, etc.,\ for $u$, and noting that $\mathrm{tr}(A)=e_1^{\prime }Ae_1+\cdots+e_d^{\prime }Ae_d$ for $d\times d$ matrix $A$, this implies that
In view of (ref), $h_\text{SEE}^*=h^*$ if $h^*$ exists. Unfortunately, it is easy to show that no single $h$ can minimize the objective function in (ref) for all $u\in \mathbb{R}^{d}$. Thus, we have to redefine the optimality with respect to the direction of $u$. The direction depends on which linear combination of $\beta $ is the focus of interest, as $u^{\prime }\mathopen{}\mathclose \bgroup \originalleft[ nh^{2r}\mathbb E(B)\mathbb E(B') - h\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright) \aftergroup \egroup \originalright] u$ is the high-order AMSE of $c^{\prime }\sqrt{n}(\hat{\beta}-\beta _{0})$ for $c=\Sigma_{XZ} \mathopen{}\mathclose \bgroup \originalleft( V^{-1/2}\aftergroup \egroup \originalright)^{\prime} u$.
Suppose we are interested in only one linear combination. Let $h_{c}^{\ast }$ be the optimal $h$ that minimizes the high-order AMSE of $c^{\prime }\sqrt{n} (\hat{\beta}-\beta _{0})$. Then
for $u=\mathopen{}\mathclose \bgroup \originalleft( V^{1/2}\aftergroup \egroup \originalright) ^{\prime } \Sigma_{XZ}^{-1} c$. Some algebra shows that
So although $h_{c}^{\ast }$ depends on $c$ via $u$, it is nevertheless greater than zero.
Now suppose without loss of generality we are interested in $d$ directions $ \mathopen{}\mathclose \bgroup \originalleft( c_{1},\ldots ,c_{d}\aftergroup \egroup \originalright) $ jointly where $c_{i}\in \mathbb{R}^{d}$. In this case, it is reasonable to choose $h_{c_{1},\ldots ,c_{d}}^{\ast }$ to minimize the sum of direction-wise AMSEs, i.e.,
where $u_{i}=\mathopen{}\mathclose \bgroup \originalleft( V^{1/2}\aftergroup \egroup \originalright) ^{\prime }\Sigma _{XZ}^{-1}c_{i}$. It is easy to show that
As an example, consider $u_{i}=e_{i}=\mathopen{}\mathclose \bgroup \originalleft( 0,\ldots ,1,\ldots ,0\aftergroup \egroup \originalright) $, the $i$th unit vector in $\mathbb{R}^{d}$. Correspondingly
It is clear that
so all three selections coincide with each other. A special case of interest is when $Z=X$, non-constant regressors are pairwise independent and normalized to mean zero and variance one, and $U \mathpalette{\protect \independenT}{\perp} X$. Then $u_{i}=c_{i}=e_{i}$ and the $d$ linear combinations reduce to the individual elements of $\beta $.
The above example illustrates the relationship between $h_{c_{1},\ldots ,c_{d}}^{\ast }$ and $h_{\text{SEE}}^{\ast }$. While $h_{c_{1},\ldots ,c_{d}}^{\ast }$ is tailored toward the flexible linear combinations $ \mathopen{}\mathclose \bgroup \originalleft(c_{1},\ldots ,c_{d}\aftergroup \egroup \originalright)$ of the parameter vector, $h_{\text{SEE} }^{\ast }$ is tailored toward the fixed $\mathopen{}\mathclose \bgroup \originalleft( \tilde{c}_{1},\ldots ,\tilde{ c}_{d}\aftergroup \egroup \originalright) $. While $h_{c_{1},\ldots ,c_{d}}^{\ast }$ and $h_{\text{SEE} }^{\ast }$ are of the same order of magnitude, in general there is no analytic relationship between $h_{c_{1},\ldots ,c_{d}}^{\ast }$ and $h_{ \text{SEE}}^{\ast }$.
To shed further light on the relationship between $h_{c_{1},\ldots ,c_{d}}^{\ast }$ and $h_{\text{SEE}}^{\ast }$, let $\mathopen{}\mathclose \bgroup \originalleft \{ \lambda _{k},k=1,\ldots ,d\aftergroup \egroup \originalright \} $ be the eigenvalues of $nh^{2r}\mathbb E(B)\mathbb E(B') - h\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright)$ with the corresponding orthonormal eigenvectors $\mathopen{}\mathclose \bgroup \originalleft \{ \ell _{k},k=1,\ldots ,d\aftergroup \egroup \originalright \} $. Then we have $nh^{2r}\mathbb E(B)\mathbb E(B') - h\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright) =\sum_{k=1}^{d}\lambda _{k}\ell _{k}\ell _{k}^{\prime }$ and $u_{i}=\sum_{j=1}^{d}u_{ij}\ell _{j}$ for $ u_{ij}=u_{i}^{\prime }\ell _{j}$. Using these representations, the objective function underlying $h_{c_{1},\ldots ,c_{d}}^{\ast }$ becomes
That is, $h_{c_{1},\ldots ,c_{d}}^{\ast }$ minimizes a weighted sum of the eigenvalues of $nh^{2r}\mathbb E(B)\mathbb E(B') - h\mathbb E\mathopen{}\mathclose \bgroup \originalleft( AA'\aftergroup \egroup \originalright)$ with weights depending on $c_{1},\ldots ,c_{d}$. By definition, $h_{\text{SEE}}^{\ast }$ minimizes the simple unweighted sum of the eigenvalues, viz.\ $\sum_{j=1}^{d}\lambda _{j}$. While $h_{\text{SEE} }^{\ast }$ may not be ideal if we know the linear combination(s) of interest, it is a reasonable choice otherwise.
In empirical applications, we can estimate $h_{c_{1},\ldots ,c_{d}}^{\ast }$ using a parametric plug-in approach similar to our plug-in implementation of $h_{\text{SEE}}^{\ast }$. If we want to be agnostic about the directional vectors $c_{1},\ldots ,c_{d}$, we can simply use $h_{\text{SEE}}^{\ast }$.
We revisit the IV-QR analysis of Job Training Partnership Act (JTPA) data in AbadieEtAl2002, specifically their Table III.\footnote{ Their data and Matlab code for replication are helpfully provided online in the Angrist Data Archive, \url{http://economics.mit.edu/faculty/angrist/data1/data/abangim02}.} They use 30-month earnings as the outcome, randomized offer of JTPA services as the instrument, and actual enrollment for services as the endogenous treatment variable. Of those offered services, only around 60 percent accepted, so self-selection into treatment is likely. Section 4 of AbadieEtAl2002 provides much more background and descriptive statistics.
We compare estimates from a variety of methods.\footnote{Code and data for replication is available on the first author's website.} “AAI” is the original paper's estimator. AAI restricts $X$ to have finite support (see condition (iii) in their Theorem 3.1), which is why all the regressors in their example are binary. Our fully automated plug-in estimator is “SEE ($\hat{h}$).” “CH” is ChernozhukovHansen2006. Method “tiny $h$” uses $h=400$ (compared with our plug-in values on the order of $10\,000$), while “huge $h$” uses $h=5\times 10^{6}$. 2SLS is the usual (mean) two-stage least squares estimator, put in the $q=0.5$ column only for convenience of comparison.
Table (ref) shows results from the sample of $5102$ adult men, for a subset of the regressors used in the model. Not shown in the table are coefficient estimates for dummies for Hispanic, working less than 13 weeks in the past year, five age groups, originally recommended service strategy, and whether earnings were from the second follow-up survey. CH is very close to “tiny $h$”; that is, simply using the smallest possible $h$ with SEE provides a good approximation of the unsmoothed estimator in this case. Demonstrating our theoretical results in Section (ref), “huge $h$” is very close to 2SLS for everything except the constant term for $q\neq 0.5$. The IVQR-SEE estimator using our plug-in bandwidth has some economically significant differences with the unsmoothed estimator. Focusing on the treatment variable (“Training”), the unsmoothed median effect estimate is below $400$ (dollars), whereas SEE$(\hat{h})$ yields $1080$, both of which are smaller than AAI's $1544$ (AAI is the most positive at all quantiles). For the $0.15$-quantile effect, the unsmoothed estimates are actually slightly negative, while SEE$(\hat{h})$ and AAI are slightly positive. For $q=0.85$, though, the SEE$(\hat{h})$ estimate is smaller than the unsmoothed one, and the two are quite similar for $q=0.25$ and $q=0.75$; there is no systematic ordering.
Computationally, our code takes only one second total to calculate the plug-in bandwidths and coefficient estimates at all five quantiles. Using the fixed $h=400$ or $h=5\times10^6$, computation is immediate.
Table (ref) shows estimates of the endogenous coefficient when various “continuous” control variables are added, specifically a quartic polynomial in age (replacing the age range dummies), baseline weekly hours worked, and baseline hourly wage.\footnote{Additional JTPA data downloaded from the W.E.\ Upjohn Institute at \url{http://upjohn.org/services/resources/employment-research-data-center/national-jtpa-study}; variables are named age, bfhrswrk, and bfwage in file expbif.dta.} The estimates do not change much; the biggest difference is for the unsmoothed estimate at the median. Our code again computes the plug-in bandwidth and SEE coefficient estimates at all five quantiles in one second. Using the small $h=400$ bandwidth now takes nine seconds total (more iterations of fsolve are needed); $h=5\times10^6$ still computes almost immediately.
For our simulation study,\footnote{ Code to replicate our simulations is available on the first author's website.} we use $G(u)$ given in (ref) as in Horowitz1998 and Whang2006. This satisfies Assumption (ref) with $r=4$. Using (the integral of) an Epanechnikov kernel with $r=2$ also worked well in the cases we consider here, though never better than $r=4$. Our error distributions always have at least four derivatives, so $r=4$ working somewhat better is expected. Selection of optimal $r$ and $G(\cdot)$, and the quantitative impact thereof, remain open questions.
We implement a plug-in version ($\hat h$) of the infeasible $h^{\ast }\equiv h_{\text{ SEE}}^{\ast }$. We make the plug-in assumption $U \mathpalette{\protect \independenT}{\perp} Z$ and parameterize the distribution of $U$. Our current method, which has proven quite accurate and stable, fits the residuals from an initial $h=(2nr)^{-1/(2r-1)}$ IV-QR to Gaussian, $t$, gamma, and generalized extreme value distributions via maximum likelihood. With the distribution parameter estimates, $f_{U}(0)$ and $f_{U}^{(r-1)}(0)$ can be computed and plugged in to calculate $\hat{h}$. With larger $n$, a nonparametric kernel estimator may perform better, but nonparametric estimation of $f_U^{(r-1)}(0)$ will have high variance in smaller samples. Viewing the unsmoothed estimator as a reference point, potential regret (of using $\hat h$ instead of $h=0$) is largest when $\hat h$ is too large, so we separately calculate $\hat{h}$ for each of the four distributions and take the smallest. Note that this particular plug-in approach works well even under heteroskedasticity and/or misspecification of the error distribution: DGPs 3.1--3.6 in Section (ref) have error distributions other than these four, and DGPs 1.3, 2.2, 3.3--3.6 are heteroskedastic, as are the JTPA-based simulations. For the infeasible $h^{\ast }$, if the PDF derivative in the denominator is zero, it is replaced by $0.01$ to avoid $h^{\ast }=\infty $.
For the unsmoothed IV-QR estimator, we use code based on ChernozhukovHansen2006 from the latter author's website. We use the option to let their code determine the grid of possible endogenous coefficient values from the data. This code in turn uses the interior point method in rq.m (developed by Roger Koenker, Daniel Morillo, and Paul Eilers) to solve exogenous QR linear programs.
We use two DGPs based on the JTPA data examined in Section (ref). The first DGP corresponds to the variables used in the original analysis in AbadieEtAl2002. For individual $i$, let $Y_{i}$ be the scalar outcome (30-month earnings), $X_{i}$ be the vector of exogenous regressors, $D_{i}$ be the scalar endogenous training dummy, $Z_{i}$ be the scalar instrument of randomized training offer, and $U_{i}\sim \text{Unif} (0,1)$ be a scalar unobservable term. We draw $X_{i}$ from the joint distribution estimated from the JTPA data. We randomize $Z_{i}=1$ with probability $0.67$ and zero otherwise. If $Z_{i}=0$, then we set the endogenous training dummy $D_{i}=0$ (ignoring that in reality, a few percent still got services). If $Z_{i}=1$, we set $D_{i}=1$ with a probability increasing in $U_{i}$. Specifically, $P(D_{i}=1\mid Z_{i}=1,U_{i}=u)=\min \{1,u/0.75\}$, which roughly matches the $P(D_{i}=1\mid Z_{i}=1)=0.62$ in the data. This corresponds to a high degree of self-selection into treatment (and thus endogeneity). Then, $Y_{i}=X_{i}\beta _{X}+D_{i}\beta _{D}(U_{i})+G^{-1}(U_{i})$, where $\beta _{X}$ is the IVQR-SEE $\hat{\beta} _{X}$ from the JTPA data (rounded to the nearest $500$), the function $\beta _{D}(U_{i})=2000U_{i}$ matches $\hat{\beta}_{D}(0.5)$ and the increasing pattern of other $\hat{\beta}_{D}(q)$, and $G^{-1}(\cdot )$ is a recentered gamma distribution quantile function with parameters estimated to match the distribution of residuals from the IVQR-SEE estimate with JTPA data. In each of $1000$ simulation replications, we generate $n=5102$ iid observations.
For the second DGP, we add a second endogenous regressor (and instrument) and four exogenous regressors, all with normal distributions. Including the intercept and two endogenous regressors, there are $20$ regressors. The second instrument is $Z_{2i}\overset{iid}{\sim}N(0,1)$, and the second endogenous regressor is $D_{2i}=0.8Z_{2i}+0.2\Phi^{-1}(U_i)$. The coefficient on $D_{2i}$ is $1000$ at all quantiles. The new exogenous regressors are all standard normal and have coefficients of $500$ at all quantiles. To make the asymptotic bias of 2SLS relatively more important, the sample size is increased to $n=50\,000$.
Table (ref) shows results for the first JTPA-based DGP, for three estimators of the endogenous coefficient: ChernozhukovHansen2006; SEE with our data-dependent $\hat h$; and 2SLS. The first and third can be viewed as limits of IVQR-SEE estimators as $h\to0$ and $h\to\infty$, respectively. We show median bias and “robust MSE,” which is squared median bias plus the square of the interquartile range divided by $1.349$, $\textrm{Bias}_{\textrm{median}}^{2}+(\textrm{IQR}/1.349)^2$. We report these “robust” versions of bias and MSE since the (mean) IV estimator does not even possess a first moment in finite samples Kinal1980. We are unaware of an analogous result for IV-QR but remain wary of presenting bias and MSE results for IV-QR, too, especially since the IV estimator is the limit of the SEE IV-QR estimator as $h\to\infty$. At all quantiles, for all methods, the robust MSE is dominated by the IQR rather than bias. Consequently, even though the 2SLS median bias is quite large for $q=0.15$, it has less than half the robust MSE of SEE($\hat h$), which in turn has half the robust MSE of the unsmoothed estimator. With only a couple exceptions, this is the ordering among the three methods' robust MSE at all quantiles. Although the much larger bias of 2SLS than that of SEE($\hat h$) or the unsmoothed estimator is expected, the smaller median bias of SEE($\hat h$) than that of the unsmoothed estimator is surprising. However, the differences are not big, and they may be partly due to the much larger variance of the unsmoothed estimator inflating the simulation error in the simulated median bias, especially for $q=0.15$. The bigger difference is the reduction in variance from smoothing.
Table (ref) shows results from the second JTPA-based DGP. The first estimator is now a nearly-unsmoothed SEE estimator instead of the unsmoothed ChernozhukovHansen2006 estimator. Although in principle ChernozhukovHansen2006 can be used with multiple endogenous coefficients, the provided code allows only one, and Tables (ref) and (ref) show that SEE with $h=400$ produces very similar results in the JTPA data. For the binary endogenous regressor's coefficient, the 2SLS estimator now has the largest robust MSE since the larger sample size reduces the variance of all three estimators but does not reduce the 2SLS median bias (since it has first-order asymptotic bias). The plug-in bandwidth yields smaller robust MSE than the nearly-unsmoothed $h=400$ at four of five quantiles. At the median, for example, compared with $h=400$, $\hat h$ slightly increases the median bias but greatly reduces the dispersion, so the net effect is to reduce robust MSE. This is consistent with the theoretical results. For the continuous endogenous regressor's coefficient, the same pattern holds for the nearly-unsmoothed and $\hat h$-smoothed estimators. Since this coefficient is constant across quantiles, the 2SLS estimator is consistent and very similar to the SEE estimators with $q=0.5$.
For exogenous QR, smoothing the criterion function (SCF) is a different approach, as discussed. The following simulations compare the MSE of our SEE estimator with that of the SCF estimator. All DGPs have $n=50$, $X_{i} \overset{iid}{\sim }\text{Unif}(1,5)$, $U_{i}\overset{iid}{\sim }N(0,1)$, $ X_{i}\mathpalette{\protect \independenT}{\perp}U_{i}$, and $ Y_{i}=1+X_{i}+\sigma (X_{i})\mathopen{}\mathclose \bgroup \originalleft( U_{i}-\Phi ^{-1}(q)\aftergroup \egroup \originalright) $. DGP 1 has $ q=0.5$ and $\sigma (X_{i})=5$. DGP 2 has $q=0.25$ and $\sigma (X_{i})=1+X_{i} $. DGP 3 has $q=0.75$ and $\sigma (X_{i})=1+X_{i}$. In addition to using our plug-in $\hat{h}$, we also compute the estimators for a much smaller bandwidth in each DGP: $h=1$, $h=0.8$, and $h=0.8$, respectively. Each simulation ran $1000$ replications. We compare only the slope coefficient estimators.
Table (ref) shows MSE and bias for the SEE and SCF estimators, for our plug-in $\hat h$ as well as the small, fixed $h$ mentioned above. The SCF estimator can have slightly lower MSE, as in the third DGP ($q=0.75$ with heteroskedasticity), but the SEE estimator has more substantially lower MSE in more DGPs, including the homoskedastic conditional median DGP. The differences are quite small with the small $h$, as expected. Deriving and implementing an MSE-optimal bandwidth for the SCF estimator could shrink the differences, but based on these simulations and the theoretical comparison in Section (ref), such an effort seems unlikely to yield improvement over the SEE estimator.
We tried additional data generating processes (DGPs). The first three DGPs are for exogenous QR, taken directly from Horowitz1998. In each case, $q=0.5$, $Y_i=X_i'\beta_0+U_i$, $\beta_0=(1,1)'$, $X_i=(1,x_i)'$ with $x_i\stackrel{iid}{\sim}\textrm{Uniform}(1,5)$, and $n=50$. In DGP 1.1, the $U_i$ are sampled iid from a $t_3$ distribution scaled to have variance two. In DGP 1.2, the $U_i$ are iid from a type I extreme value distribution again scaled and centered to have median zero and variance two. In DGP 1.3, $U_i=(1+x_i)V/4$ where $V_i\stackrel{iid}{\sim}N(0,1)$.
DGPs 2.1, 2.2, and 3.1--3.6 are shown in the working paper version; they include variants of the Horowitz1998 DGPs with $q\ne0.5$, different error distributions, and another regressor.
DGPs 4.1--4.3 have endogeneity. DGP 4.1 has $q=0.5$, $n=20$, and $\beta _{0}=(0,1)^{\prime }$. It uses the reduced form equations in CattaneoEtAl2012 with $\gamma _{1}=\gamma _{2}=1$, $x_{i}=1$, $z_{i}\sim N(0,1)$, and $\pi =0.5$. Similar to their simulations, we set $\rho =0.5$, $(\tilde{v}_{1i},\tilde{v}_{2i})$ iid $ N(0,1)$, and $(v_{1i},v_{2i})^{\prime }=\mathopen{}\mathclose \bgroup \originalleft(\tilde{v}_{1i},\sqrt{1-\rho ^{2}} \tilde{v}_{2i}+\rho \tilde{v}_{1i}\aftergroup \egroup \originalright)'$. DGP 4.2 is similar to DGP 4.1 but with $(\tilde v_{1i},\tilde v_{2i})'$ iid Cauchy, $n=250$, and $\beta_0=\mathopen{}\mathclose \bgroup \originalleft(0,\mathopen{}\mathclose \bgroup \originalleft[\rho-\sqrt{1-\rho^2}\aftergroup \egroup \originalright]^{-1}\aftergroup \egroup \originalright)'$. DGP 4.3 is the same as DGP 4.1 but with $q=0.35$ (and consequent re-centering of the error term) and $n=30$.
We compare MSE for our SEE estimator using the plug-in $\hat h$ and estimators using different (fixed) values of $h$. We include $h=0$ by using unsmoothed QR or the method in ChernozhukovHansen2006 for the endogenous DGPs. We also include $h=\infty$ (although not in graphs) by using the usual IV estimator. For the endogenous DGPs, we consider both MSE and the “robust MSE” defined in Section (ref) as $\textrm{Bias}_{\textrm{median}}^{2}+(\textrm{IQR}/1.349)^2$.
For “size-adjusted” power (SAP) of a test with nominal size $\alpha$, the critical value is picked as the $(1-\alpha)$-quantile of the simulated test statistic distribution. This is for demonstration, not practice. The size adjustment fixes the left endpoint of the size-adjusted power curve to the null rejection probability $\alpha $. The resulting size-adjusted power curve is one way to try to visualize a combination of type I and type II errors, in the absence of an explicit loss function. One shortcoming is that it does not reflect the variability/uniformity of size and power over the space of parameter values and DGPs.
Regarding notation in the size-adjusted power figures, the vertical axis in the size-adjusted power figures shows the simulated rejection probability. The horizontal axis shows the magnitude of deviation from the null hypothesis, where a randomized alternative is generated in each simulation iteration as that magnitude times a random point on the unit sphere in $\mathbb{R}^{d}$, where $\beta \in \mathbb{R}^{d}$. As the legend shows, the dashed line corresponds to the unsmoothed estimator ($h=0$), the dotted line to the infeasible $h_{ \text{SEE}}^{\ast }$, and the solid line to the plug-in $\hat{h}$.
For the MSE graphs, the flat horizontal solid and dashed lines are the MSE of the intercept and slope estimators (respectively) using feasible plug-in $ \hat{h}$ (recomputed each replication). The other solid and dashed lines (that vary with $h$) are the MSE when using the value of $h$ from the horizontal axis. The left vertical axis shows the MSE values for the intercept parameter; the right vertical axis shows the MSE for slope parameter(s); and the horizontal axis shows a log transformation of the bandwidth, $\log _{10}(1+h)$.
Our plug-in bandwidth is quite stable. The range of $\hat{h}$ values over the simulation replications is usually less than a factor of $10$, and the range from $0.05$ to $0.95$ empirical quantiles is around a factor of two. This corresponds to a very small impact on MSE; note the log transformation in the x-axis in the MSE graphs.
In DGPs 1.1--1.3, SEE($\hat h$) has smaller MSE than either the unsmoothed estimator or OLS, for both the intercept and slope coefficients. Figure (ref) shows MSE for DGPs 1.1 and 1.3. It shows that the MSE of SEE($\hat h$) is very close to that of the best estimator with a fixed $h$. In principle, a data-dependent $\hat h$ can attain MSE even lower than any fixed $h$. SAP for SEE($\hat h$) is similar to that with $h=0$; see Figure (ref) for DGPs 1.1 and 1.3.
Figures (ref) and (ref) show MSE and “robust MSE” for two DGPs with endogeneity. Graphs for the other endogenous DGP (4.1) are similar to those for the slope estimator in DGP 4.3 but with larger MSE; they may be found in the working paper. The MSE graph for DGP 4.2 is not as informative since it is sensitive to very large outliers that occur in only a few replications. However, as shown, the MSE for SEE($\hat h$) is still better than that for the unsmoothed IV-QR estimator, and it is nearly the same as the MSE for the mean IV estimator (not shown: $1.1\times10^6$ for $\beta_1$, $2.1\times10^5$ for $\beta_2$). For robust MSE, SEE($\hat h$) is again always better than the unsmoothed estimator. For DGP 4.3 with normal errors and $q=0.35$, it is similar to the IV estimator, slightly worse for the slope coefficient and slightly better for the intercept, as expected. Also as expected, for DGP 4.2 with Cauchy errors, SEE($\hat h$) is orders of magnitude better than the mean IV estimator. Overall, using $\hat{h}$ appears to consistently reduce the MSE of all estimator components compared with $h=0$ and with IV ($h=\infty $). Almost always, the exception is cases where MSE is monotonically decreasing with $h$ (mean regression is more efficient), in which $\hat{h}$ is much better than $h=0$ but not quite large enough to match $h=\infty $.
Figure (ref) shows SAP for DGPs 4.1 and 4.3. The gain from smoothing is more substantial than in the exogenous DGPs, close to $10$ percentage points for a range of deviations. Here, the randomness in $\hat h$ is not helpful. In DGP 4.2 (not shown), the SAP for $\hat h$ is actually a few percentage points below that for $h=0$ (which in turn is below the infeasible $h^*$), and in DGP 4.1, the SAP improvement from using the infeasible $h^*$ instead of $\hat h$ is similar in magnitude to the improvement from using $\hat h$ instead of $h=0$. Depending on one's loss function of type I and type II errors, the SEE-based test may be preferred or not.
We have presented a new estimator for quantile regression with or without instrumental variables. Smoothing the estimating equations (moment conditions) has multiple advantages beyond the known advantage of allowing higher-order expansions. It can reduce the MSE of both the estimating equations and the parameter estimator, minimize type I error and improve size-adjusted power of a chi-square test, and allow more reliable computation of the instrumental variables quantile regression estimator especially when the number of endogenous regressors is larger. We have given the theoretical bandwidth that optimizes these properties, and simulations show our plug-in bandwidth to reproduce all these advantages over the unsmoothed estimator. Links to mean instrumental variables regression and robust estimation are insightful and of practical use.
The strategy of smoothing the estimating equations can be applied to any model with nonsmooth estimating equations; there is nothing peculiar to the quantile regression model that we have exploited. For example, this strategy could be applied to censored quantile regression, or to select the optimal smoothing parameter in Horowitz2002's (Horowitz2002) smoothed maximum score estimator. The present paper has focused on parametric and linear IV quantile regression; extensions to nonlinear IV quantile regression and nonparametric IV quantile regression along the lines of ChenPouzo2009,ChenPouzo2012 are currently under development.
\theendnotes