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.
80,324 characters · 22 sections · 66 citation commands
Smoothed instrumental variables quantile regression
\inserttype[st0001]{article}
{ David M.\ Kaplan\\Department of Economics\\University of Missouri\\Columbia, MO\\[email removed] }
The instrumental variables quantile regression (IVQR) model of ChernozhukovHansen2005 has become popular for modeling causal effects that differ across individuals (or firms, countries, etc.)\ due to unobserved variables (rather than interactions among observed variables), but the existing IVQR \Stata\ commands \stcmd{ivqreg} Stata.ivqreg and \stcmd{ivqreg2} Stata.ivqreg2 have several limitations. First, \stcmd{ivqreg} ChernozhukovHansen2006 only allows a single endogenous term. Besides excluding multiple endogenous regressors, this also excludes interaction or quadratic terms involving the endogenous regressor, for example. Second, \stcmd{ivqreg2} MachadoSantosSilva2019 imposes a location--scale model that can help if correct but that can result in inconsistency if misspecified. \Citet[p.\ 159]{MachadoSantosSilva2019} write, “In these linear models, the validity of [our model] depends on assumptions that are stronger than those required by the IVQR but, when these assumptions are valid, [our estimator] has some potential advantages,” like ensuring non-crossing of structural quantile functions. Ideally, \stcmd{sivqr} and \stcmd{ivqreg2} could be combined using a “model averaging” approach like in ChengLiaoShi2019 or Liu2019; meanwhile, simply running both may be insightful. Third, the existing commands' computation time can be long. Fourth, \stcmd{ivqreg} and \stcmd{ivqreg2} share some minor inconveniences, like not supporting syntax for factor variables and interaction terms. Fifth, \stcmd{ivqreg} is seemingly not actively supported and has additional issues; for example, the newest version 1.0.1 (02dec2010) fixes a problem from version 1.0.0 displaying coefficient labels when using an asterisk in the regressor list (like \stcmd{xvar*}), but it introduces an error when there is overidentification (more excluded instruments than endogenous regressors).
To address these limitations, I introduce the \stcmd{sivqr} command. The \stcmd{sivqr} command allows models with multiple endogenous terms, supports convenient syntax like for factor variables and interactions terms, and computes a consistent estimator of the IVQR parameters in a reasonable amount of time. The runtime in the section (ref) example is around $25$ seconds for \stcmd{sivqr}, as opposed to a few minutes (\stcmd{ivqreg}) or over $20$ minutes (\stcmd{ivqreg2}). Although these differences become trivial on small datasets where everything computes quickly, they become even more important when using \stcmd{bootstrap} for cluster-robust standard errors, which requires computing the estimator many times.
The \stcmd{sivqr} command implements the smoothed IVQR estimator of KaplanSun2017, using familiar syntax. Syntax is similar to the \stcmd{ivregress} command (see \rref{ivregress}) but additionally specifying a quantile level. Likewise, syntax is similar to the \stcmd{qreg} command (see \rref{qreg}) but additionally specifying instruments for the endogenous regressors. With non-i.i.d.\ sampling, the \stcmd{bootstrap} prefix can be used for standard errors; see \rref{bootstrap} and the example in section (ref).
Other \Stata\ commands that address endogeneity in quantile models are based on different models than that of ChernozhukovHansen2005. These alternative models have neither “stronger” nor “weaker” assumptions; they may be more plausible in certain settings and less plausible in others. The \stcmd{ivqte} command FrolichMelly2010 estimates the unconditional (default) or conditional (option \stcmd{aai}) “local quantile treatment effect,” analogous to the local average treatment effect model of ImbensAngrist1994. This conditional (\stcmd{aai}) estimand and the IVQR estimand are different but closely related, as explained in section 10.5 of MellyWuthrich2017 and detailed in Wuthrich2020, both of whom compare the models more generally, too. See also FrolichMelly2013 for details. Practically, compared to \stcmd{ivqte} with the \stcmd{aai} option, the main advantage of \stcmd{sivqr} is the ability to handle more than a single, binary endogenous regressor. The \stcmd{cqiv} command ChernozhukovEtAl2019,Stata.cqiv uses a control function estimator based on a triangular model Lee2007, also allowing for censoring ChernozhukovEtAl2015. The estimand is the same as for IVQR (see section (ref)). Practically, compared to \stcmd{cqiv}, the main advantages of \stcmd{sivqr} are handling multiple and/or non-continuous endogenous regressors, as well as allowing simultaneity and reverse causality. For further comparison of the IVQR and triangular models, see section 9.2.5 of ChernozhukovHansenWuthrich2017.
Unfortunately, all these commands (including \stcmd{sivqr}) require “strong” instruments for valid standard errors. In the future, it would be valuable to have \Stata\ implementations of IVQR inference methods robust to weak instruments, such as those of ChernozhukovHansen2008, ChernozhukovEtAl2009, ChernozhukovHansenWuthrich2017, and references therein. See also KaplanLiu2021 for a “$k$-class” IVQR estimator that can reduce estimation bias from weak instruments, as well as simulation evidence about the usefulness of the conventional “first-stage” \fstat\ for assessing instrument strength in IVQR.
Section (ref) discusses the methodology at a relatively intuitive level. Section (ref) describes syntax and usage of \stcmd{sivqr}. Section (ref) provides examples that can be replicated with the provided do-file. Section (ref) shows some of the theoretical foundations before concluding. Abbreviations are used for instrumental variables (IV), quantile regression (QR), IV quantile regression (IVQR), and two-stage least squares (2SLS).
This section discusses methodology at a relatively non-technical level (compared to section (ref)).
First, consider interpretation of the parameters being estimated. For simplicity, imagine a single regressor $x$ and outcome variable $y$, both scalars. Scalars $u$ and $v$ represent unobserved variables.
Usually there is a structural model like $y=\beta_0+\beta_1x+v$. Here, $\beta_0$ and $\beta_1$ are unknown constants and $v$ is everything besides $x$ that causally determines $y$. Here, $\beta_1$ has a causal interpretation as some effect of $x$ on $y$, but in reality we rarely believe such an effect is the same for all individuals (or firms, or schools, etc.). There are two approaches: either let differences in effects go into $v$ and interpret $\beta_1$ as some sort of average, or try to learn about the effect heterogeneity.
Consider a structural model that allows individuals to each have their own intercept and slope. Since the parameters are now individual-specific, they are not constants like $\beta_0$ and $\beta_1$ were, but rather random variables in a “random coefficients” model. The additive error term $v$ from before is simply absorbed into the random intercept, so the structural model is $y=b_0+b_1x$. That is, each individual has their own $(y,x,b_0,b_1)$, but only $(y,x)$ is observable.
Now imagine the random coefficients can each be written as deterministic functions of a scalar unobservable $u$: $b_0=\beta_0(u)$ and $b_1=\beta_1(u)$. Because $\beta_0(\cdot)$ and $\beta_1(\cdot)$ are unrestricted, the distribution of $u$ can be normalized to uniform over the unit interval $[0,1]$. The functions $\beta_0(\cdot)$ and $\beta_1(\cdot)$ are unknown but deterministic; evaluated at a fixed $0<\tau<1$, $\beta_0(\tau)$ and $\beta_1(\tau)$ are unknown constants, just like $\beta_0$ and $\beta_1$ were before. The differences across individuals are driven by $u$. Each individual has their own $(y,x,u)$, with $y=\beta_0(u)+\beta_1(u)x$, whereas functions $\beta_0(\cdot)$ and $\beta_1(\cdot)$ are not specific to any individual.
A special case of this random coefficient model is the usual structural model with constant (non-random) coefficients. Let $\beta_1(u)=\beta_1$, a constant that does not depend on $u$. Define $v\equiv\beta_0(u)-\beta_0$, so $\beta_0(u)=\beta_0+v$. The function $\beta_0(\cdot)$ can be interpreted as the inverse CDF (quantile function) of $v$, shifted by $\beta_0$; for example, with $\Phi(\cdot)$ as the standard normal CDF, $v=\Phi^{-1}(u)$ has a standard normal distribution. Then \[ y=\beta_0(u)+\beta_1(u)x=\beta_0+v+\beta_1x = \beta_0+\beta_1x+v. \]
Some additional restrictions are required in order to learn about the structural model. Imagine further that given $x$, $\beta_0(u)+\beta_1(u)x$ is increasing in $u$. This is known as a “monotonicity” assumption. It also explains why $u$ is often called the “rank variable”: it describes how somebody's $y$ would rank in the population if everyone were forced to have the same $x$. For example, somebody with $u=0.5$ would have median $y$, and somebody with $u=0.9$ would have $90$th percentile $y$. If everybody keeps the same “rank” regardless of the $x$ value, then “rank invariance” holds. A weaker assumption called “rank similarity” allows the ranking to differ across $x$ as long as the differences are exogenous.
Even if $x$ is endogenous, IVQR can estimate $\beta_0(\tau)$ and $\beta_1(\tau)$ for any $0<\tau<1$ if an instrument $z$ is available that is related to $x$ but independent of $u$ ChernozhukovHansen2005. The interpretation of these parameters depends partly on the rank assumption. If rank invariance holds, then $\beta_0(\tau)+\beta_1(\tau)x_0$ is the $y$ value that somebody with rank $u=\tau$ would have if we assign them to have value $x=x_0$. Even with the weaker rank similarity assumption, this is the $\tau$-quantile structural function of ImbensNewey2009: given any $x=x_0$, it provides the $\tau$-quantile of $\beta_0(u)+\beta_1(u)x_0$ over the unconditional population distribution of $u$ (uniform over $[0,1]$), which is $\beta_0(\tau)+\beta_1(\tau)x_0$ due to monotonicity. Similarly, $\beta_1(\tau)$ can be interpreted as a $\tau$-quantile treatment effect, capturing the difference in the $\tau$-quantile of $y$ between the counterfactual “untreated” distribution for which everyone has $x=x_0$ and the counterfactual “treated” distribution for which everyone has $x=x_0+1$.
\Citet{ChernozhukovHansen2005} show how to derive moment conditions (or “estimating equations”) to characterize the parameters $\beta_0(\tau)$ and $\beta_1(\tau)$ given a valid instrument and the assumptions discussed above. For comparison, with $v\equiv y-\beta_0-\beta_1x$ the standard IV moment conditions are
The IVQR moment conditions are
where $\operatorname{\mathbf{1}}\{\cdot\}$ is the indicator function defined as $\operatorname{\mathbf{1}}\{\mathcal{A}\}=1$ if $\mathcal{A}$ is true and otherwise $\operatorname{\mathbf{1}}\{\mathcal{A}\}=0$. This is implied by a conditional quantile restriction on $y-\beta_0(\tau)-\beta_1(\tau)x$ (given $z$), mirroring how the standard IV moments are implied by a conditional mean restriction on $y-\beta_0-\beta_1x$ (given $z$).
Smoothing solves the computational difficulties inherent in (ref). Unlike the standard IV moment conditions, the IVQR moment conditions cannot be solved explicitly for the parameters, nor are the sample moment conditions smooth (differentiable) functions of the parameters. This computational challenge is addressed by “smoothing” the indicator function: replacing it with a function that decreases continuously from $1$ to $0$ rather than discontinuously jumping from $1$ to $0$. The smoothed function differs from the indicator function only near the jump point; they are identical outside of a window whose width is controlled by a bandwidth (which asymptotically goes to zero). This smoothing allows the sample system of equations defining the estimator to be solved by standard numerical methods like those available in Mata. As a bonus, smoothing also improves the statistical properties of the estimator in theory and simulations; see sections 5 and 7 of KaplanSun2017.
The \stcmd{sivqr} command estimates the coefficients in an instrumental variables quantile regression (IVQR) model, and it provides standard errors.
Syntax, options, and stored results are now shown. Prefix \stcmd{by} is allowed (see \dref{by}). After running, postestimation commands like \rref{test} and \rref{predict} can be used as usual.
\hangpara As in \stcmd{ivregress}: \varlist$_1$ has exogenous regressors (or control variables), \varlist$_2$ has endogenous regressors, and \varlist$_\textup{iv}$ has excluded instruments (exogenous variables that instrument for \varlist$_2$).
\hangpara \stcmd{pweight}s, \stcmd{iweight}s, and \stcmd{fweight}s are allowed; see \uref{11.1.6 weight}.
\hangpara quantile(\num) specifies the quantile level as in \rref{qreg}: a real number strictly between $0$ and $1$, or alternatively a number between $1$ and $100$ interpreted as a percentile. For example, \stcmd{quantile(0.5)} specifies the median, as does \stcmd{quantile(50)}. To prevent inadvertent mistakes, there is no default; this must be specified explicitly.
\hangpara bandwidth(\num) manually specifies the desired smoothing bandwidth. Leaving it unspecified invokes a plug-in bandwidth based on KaplanSun2017 as the default. If the desired bandwidth is too small to find a numerical solution, then it is increased until a solution is found. For example, specifying \stcmd{bandwidth(0)} requests the smallest bandwidth that is computationally feasible. Generally, the default plug-in bandwidth is recommended; compare to an arbitrarily chosen bandwidth, the plug-in usually produces a more accurate estimator (lower mean squared error). However, for a final analysis, it is prudent to also try a smaller bandwidth like \stcmd{bandwidth(0)} as a sensitivity check (results should be qualitatively similar, though not identical). A smaller bandwidth tends to yield an estimator with larger variance but smaller (closer to zero) bias.
\hangpara level(\num) specifies the confidence level, as a percentage, for confidence intervals. The default is \stcmd{level(95)} or as set by \stcmd{set level}; see \uref{20.7 Specifying the width of confidence intervals}.
\hangpara reps(\num) specifies the number of bootstrap replications for estimating the variance--covariance matrix and standard errors. With the default \stcmd{reps(0)}, heteroskedasticity-robust analytic standard errors are reported. Otherwise, the Bayesian bootstrap of Rubin1981 is used; it is a valid frequentist bootstrap that also has a nonparametric Bayesian interpretation (see section (ref)).
\hangpara logiterations prints each iteration of the numerical solver; see \mrefe{solvenl(\,)}, in particular \stcmd{solvenl_init_iter_log(\,)}. The default is not to print such information, which usually only helps debugging and troubleshooting.
\hangpara noconstant omits the intercept term that otherwise is included automatically.
\hangpara seed(\num) sets the random-number seed, to make results replicable; default is \stcmd{seed(112358)}. The current seed is restored at the end of execution.
\hangpara nodots suppresses display of the replication dots (see \rref{bootstrap}).
\hangpara initial( matname) sets the initial coefficient values for the numerical search, using the values stored in the matrix (row vector) named matname. If not specified, then \stcmd{qreg} is used to generate initial values. This option is particularly helpful when estimating the same model across many quantile levels; examples are provided in the sivqr_examples_init.do file.
\stcmd{sivqr} stores the following in e():
Almost all of the above are the same as for the familiar commands \stcmd{ivregress} and/or \stcmd{qreg} (or \stcmd{bsqreg}). The only exceptions are for the smoothing bandwidth. By default, the plug-in bandwidth is computed and returned in \stcmd{e(bwidth_req)}. More specifically, three plug-in bandwidths are computed (see section (ref)), of which \stcmd{e(bwidth_req)} is the minimum and \stcmd{e(bwidth_max)} is the maximum. If the numerical solver cannot find the solution with \stcmd{e(bwidth_req)} because it is too close to zero, then the bandwidth is increased until the numerical solver finds the solution. This feasible bandwidth is returned in \stcmd{e(bwidth)}. There is nothing wrong with these values being different; for example, specifying \stcmd{bandwidth(0)} makes \stcmd{e(bwidth_req)} zero and \stcmd{e(bwidth)} the smallest possible bandwidth for which the numerical solver finds a solution.
If the default plug-in bandwidth is used and \stcmd{e(bwidth)} is greater than \stcmd{e(bwidth_max)}, then there may be deeper problems such as weak instruments. Although there is currently no weak instrument test specific to IVQR, it would be helpful to run a weak instrument test for 2SLS in this case, like by running \stcmd{ivregress} and then \stcmd{estat firststage}. If \stcmd{e(bwidth)} is above \stcmd{e(bwidth_req)} but less than \stcmd{e(bwidth_max)}, then it may be simply that the plug-in bandwidth was chosen too conservatively (see section (ref)), but it may still be worth running \stcmd{ivregress} and \stcmd{estat firststage}.
\stcmd{sivqr} stores the following in r():
These are the standard values stored by \stcmd{ereturn display}; see \pref{ereturn}.
All results (besides runtimes) can be replicated with the provided sivqr_examples.do file. I ran them in \Stata/MP 17.0 with 2 cores Stata17 on Windows 10, on a few-years-old standard-issue Dell desktop computer (Intel i5-8500 3GHz processor, 8GB RAM). Besides using version 1.1.0 (22jan2022) of \stcmd{sivqr}, I used version 2.6 (29sep2021) of \stcmd{ivqreg2} from SSC Stata.ivqreg2, and version 1.0.0 (28jul2010) of \stcmd{ivqreg} Stata.ivqreg. The most recent version 1.0.1 (02dec2010) from Stata.ivqregDEC unfortunately gives an error when trying to use an overidentified model (as in section (ref)), “st_numscalar(): 3204 matrix found where scalar required.” Other changes in version 1.0.1 are both positive and negative; modifying the section (ref) model to only have \stcmd{union} as an excluded instrument, then version 1.0.1 reports standard errors whereas version 1.0.0 does not (instead reporting “Warning: variance matrix is nonsymmetric or highly singular”), although the version 1.0.1 runtime is over ten times longer.
The primary focus is comparison with \stcmd{ivqreg} and \stcmd{ivqreg2}, in terms of estimates, speed, and capabilities, including:
As noted earlier, \stcmd{ivqreg} does not allow multiple endogenous regressors, which precludes even a quadratic function of the endogenous regressor, and the factorial operator \stcmd{\#\#} is supported only by \stcmd{sivqr}.
The following revisits Example 4 in \rref{ivregress}. Briefly, it hopes to estimate a structural wage model, treating job tenure as endogenous; see \rref{ivregress} for details. I focus on the coefficient on job tenure.
For comparison, \stcmd{ivregress 2sls} is run first. This provides a helpful “sanity check”: even if there are differences across quantile levels, the IVQR coefficient estimates and signs should still be similar overall, especially around the median.
The following shows results from \stcmd{sivqr} at the median level. To cluster standard errors at the individual level, the \stcmd{bootstrap} prefix is used after using \stcmd{xtset} appropriately. \footnote{For example, see \url{https://www.stata.com/support/faqs/statistics/bootstrap-with-panel-data/}} With \stcmd{bootstrap}, it would only waste time (without changing the reported estimates or standard errors) to specify the \stcmd{reps} option (with a non-zero number) for \stcmd{sivqr}. At the median, the \stcmd{sivqr} coefficient estimates (other than the intercept) as well as the corresponding cluster-robust standard errors are similar to those from \stcmd{ivregress 2sls}.
The following code produces the results in tables (ref) and (ref), besides the cluster bootstrap standard errors computed above (that are repeated for other quantiles). Compared to the provided do-file, below I omit various \stcmd{timer} and matrix storage commands (and output). I run \stcmd{sivqr} with both the default plug-in bandwidth as well as \stcmd{bandwidth(0)}, which tries to find the smallest computationally feasible bandwidth. The squared age variable is generated manually because \stcmd{ivqreg} and \stcmd{ivqreg2} do not support the factorial operator \stcmd{\#\#}.
\sisetup{detect-mode, tight-spacing = true, group-digits = false, input-signs = , input-symbols = , input-open-uncertainty = , input-close-uncertainty = , table-align-text-pre = false, round-mode = places, round-precision = 4, table-space-text-pre = (, table-space-text-post = ), table-number-alignment=center, table-format=1.4 }
Table (ref) shows the different \stcmd{tenure} coefficient estimates at different quantiles, with the 2SLS results for comparison. To help gauge the magnitude of the differences, table (ref) also shows cluster-robust standard errors from 2SLS and from using \stcmd{bootstrap} with \stcmd{sivqr}. (None of the commands produces cluster-robust standard errors by itself, so the default standard errors are not reported; they are generally too small here, as expected.)
At the median, all estimates are similar. All are well within one standard error of the 2SLS estimate. (The 2SLS estimate can differ greatly from median IVQR estimates, just as the sample mean and median can differ greatly, but they are often similar.) The cluster-bootstrap \stcmd{sivqr} standard error is also similar to the 2SLS cluster-robust standard error.
Away from the median, \stcmd{sivqr} (either bandwidth) and \stcmd{ivqreg} remain similar to each other, but other differences arise. First, the IVQR estimates are all increasing with the quantile level, differing significantly from the 2SLS estimate both economically and statistically. Compared to the $0.25$ quantile estimates, most of the $0.75$ quantile estimates are nearly twice as large, suggesting substantial heterogeneity. This could be interpreted as individuals with higher levels of “ability” (and other unobserved wage determinants) having higher returns to job tenure.
Second, \stcmd{ivqreg2} differs from the other IVQR estimators. This is especially seen at the $0.75$ quantile: the \stcmd{ivqreg2} estimate of the return to tenure is below the other three estimates by around $0.03$ (three percentage points), which is around a three standard error difference. It is possible that the location--scale model of \stcmd{ivqreg2} is correctly specified and that this difference is only due to larger estimation error of \stcmd{sivqr} and \stcmd{ivqreg}, but given the consistency of both \stcmd{sivqr} and \stcmd{ivqreg} along with the difference being three standard errors, it seems likely that at least part of the difference is explained by misspecification of the \stcmd{ivqreg2} location--scale model.
As a side note, these \stcmd{ivqreg2} results were using the regressor order that gave the most economically plausible estimates. Keeping everything the same but with \stcmd{tenure} last in the regressor list, the estimates of the \stcmd{tenure} coefficient become significantly negative, and two warning messages appear, “WARNING: .03758389% of the fitted values of the scale function are not positive” and “Warning: variance matrix is nonsymmetric or highly singular” (with missing standard errors).
Table (ref) shows the runtime of the different commands at different quantiles. The fastest is \stcmd{sivqr} with the default plug-in bandwidth and analytic standard errors, taking around $25$ seconds to compute at all three quantiles. Using \stcmd{bandwidth(0)} is significantly slower (a few minutes total, instead of seconds) because it requires many more calls to the numerical solver, trying to find the very smallest usable bandwidth. The \stcmd{ivqreg} time is yet slower than \stcmd{bandwidth(0)}, also multiple minutes total. Naturally, running a cluster bootstrap with $100$ replications takes even longer, even using the fast \stcmd{sivqr} with plug-in bandwidth, taking around $40$ minutes total. Even without bootstrap, \stcmd{ivqreg2} takes over $20$ minutes, an order of magnitude more than the $25$ seconds taken by \stcmd{sivqr}.
The provided do-file illustrates some additional important considerations about computation time. Most of the \stcmd{sivqr} computation time is from calling the \stcmd{solvenl} numerical solver. Naturally, some numerical problems are more difficult to solve quickly. In particular, difficulty is higher with less smoothing or when there is less data near the quantile of interest (roughly speaking), which most often happens farther into the tails at quantile levels near zero or one.
Table (ref) illustrates this by running \stcmd{sivqr} over quantiles $0.10,0.15,\ldots,0.90$, using either the plug-in bandwidth (always between $0.05$ and $0.08$) or \stcmd{bandwidth(0)} (always between $0.00001$ and $0.00012$, except at the $0.9$ quantile). First, although there are other factors, runtime is generally longer in the tails: with the plug-in bandwidth, runtime is at least $8$ seconds at quantile levels $\{0.10,0.80,0.85,0.90\}$ but under $8$ seconds in between. The pattern is less clear for \stcmd{bandwidth(0)} because it depends not only on the runtime per \stcmd{solvenl} call, but also on the number of such calls. Second, the runtime is longer with the smaller bandwidth at $16$ out of $17$ quantile levels.
Also, even with the same original dataset, some bootstrap datasets are naturally more difficult than others, so bootstrap runtime can depend on the random seed. With \stcmd{quantile(0.5)}, out of $100$ runs with different seeds, runtimes for the built-in \stcmd{sivqr} bootstrap with $20$ replications ranged from $9$ seconds to $30$ seconds, with lower and upper quartiles $11$ seconds and $16$ seconds. Among five runs with the \stcmd{bootstrap} prefix (also $20$ replications), runtimes ranged from $79$ to $90$ seconds. (I conjecture this is much slower because it resamples a new dataset rather than merely reweighting the original dataset like the Bayesian bootstrap of \stcmd{sivqr}, but I have not broken down the runtime by function to investigate further.) However, it cannot be known beforehand which seeds are “faster” for a particular dataset, so using the default seed is recommended, which also avoids any appearance of trying to manipulate results.
This section contains additional theoretical details. After introducing notation in section (ref), section (ref) concerns the difficulty of estimation and the solution of KaplanSun2017, while section (ref) has some of the identification arguments from ChernozhukovHansen2005 that characterize causal parameters as solutions to moment conditions (estimating equations). Section (ref) describes various bandwidths. Section (ref) has brief notes on the Bayesian bootstrap used to compute standard errors.
These details may provide a deeper understanding for some readers, but they may also be skipped without hindering successful application of \stcmd{sivqr} in practice.
Notationally, let $y$ be the scalar outcome variable, $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ a column vector of regressors, and $u$ a scalar unobserved variable. Let $\boldsymbol{\mathbf{\MakeLowercase{d}}}$ be a subset of $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ containing all the endogenous regressors. Let $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ be a column vector of all exogenous variables: both exogenous regressors in $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ as well as the excluded instruments (instrumental variables). The scalar random variable $u$ is unobserved. Following \sj convention, vectors are lowercase bold (like $\boldsymbol{\mathbf{\MakeLowercase{x}}}$) and scalars are lowercase plain (like $y$) since uppercase is reserved for matrices, but with some luck, random variables can be distinguished from their possible non-random realizations.
Additional notation is introduced as needed below.
First, as a helpful reference point for intuition, recall the standard (“mean”) instrumental variables (IV) model and estimator. The goal is to estimate the non-random parameter vector $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$ in the structural model
If the instruments $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ are independent of the structural error $u$, then
where the $\implies$ follows from the law of iterated expectations. If the dimensions of $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ and $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ are the same (“exact identification”), then this can be solved explicitly for $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$:
and (with i.i.d.\ data) the expectations are replaced by sample averages to yield the familiar IV estimator.
Superficially, imagine replacing the conditional mean restriction $\operatorname{E}(y\mid\boldsymbol{\mathbf{\MakeLowercase{z}}})$ in (ref) with a conditional $\tau$-quantile restriction for some $0<\tau<1$:
That is, conditional on any value of $\boldsymbol{\mathbf{\MakeLowercase{z}}}$, the conditional distribution of $u$ has its $\tau$-quantile equal to zero. By definition, the $\tau$-quantile of a distribution is the value with $\tau$ probability below that value, and the same is true conditional on $\boldsymbol{\mathbf{\MakeLowercase{z}}}$. Thus, (ref) is equivalent to
Let $\operatorname{\mathbf{1}}\{\cdot\}$ be the indicator function such that $\operatorname{\mathbf{1}}\{\mathcal{A}\}=1$ if $\mathcal{A}$ is true and otherwise $\operatorname{\mathbf{1}}\{\mathcal{A}\}=0$. Rewriting $\Pr(\cdot)$ as $\operatorname{E}[\operatorname{\mathbf{1}}\{\cdot\}]$ and applying the law of iterated expectations as in (ref),
Despite some similarities, computing an IVQR estimator based on (ref) is much more difficult than computing the mean IV estimator based on (ref). With mean IV, it was possible to solve for $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$ and replace the mean (expectation) with the sample average. Here, the $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$ is stuck inside the indicator function $\operatorname{\mathbf{1}}\{\cdot\}$, so it cannot be solved for explicitly. Further, it may be impossible to solve the equation exactly after substituting in the sample average.
These difficulties are both addressed by “smoothing” the indicator function in (ref). Replacing $\operatorname{\mathbf{1}}\{\cdot\}$ with a very similar but continuously differentiable function allows the sample system of equations to be solved quickly by standard numerical methods. Specifically, replacing $\operatorname{\mathbf{1}}\{v\le0\}$, $\tilde{I}(v)$ is a function of $v$ that smoothly decreases from $1$ to $0$ over $-1\le v\le1$ instead of decreasing discontinuously from $1$ to $0$ at $v=0$. \footnote{This notation is simpler and equivalent to (but different than) that in KaplanSun2017.} Adding a bandwidth $h$ allows $\tilde{I}(v/h)$ to decrease more steeply over $-h\le v\le h$. The estimator $\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}$ solves the “smoothed estimating equations” (or “smoothed moment conditions”)
As a bonus, this smoothing also improves estimation precision, as shown both theoretically and in simulations by KaplanSun2017. Weights $w_i$ can be inserted into (ref) as needed:
Specifically, \stcmd{sivqr} uses a piecewise linear $\tilde{I}(\cdot)$ that connects the smoothed IVQR estimator to other estimators. Instead of jumping down from $1$ to $0$ discontinuously, $\tilde{I}(\cdot)$ has $\tilde{I}(v)=1$ for $v\le-1$ and $\tilde{I}(v)=0$ for $v\ge1$, but transitions linearly with $\tilde{I}(v)=(1-v)/2$ for $-1<v<1$. \Citet[\S2.2, p.\ 111]{KaplanSun2017} show how this choice produces the Winsorized mean estimator of the type in Huber1964 in the special case with $\tau=0.5$ and an intercept-only model ($\boldsymbol{\mathbf{\MakeLowercase{x}}}=\boldsymbol{\mathbf{\MakeLowercase{z}}}=1$). They also show (in Section 2.2) how using a very large amount of smoothing ($h$) turns the smoothed IVQR estimator into the usual (mean) IV estimator, with an adjusted intercept. Intuitively, if $h$ is large enough that every observation is smoothed, then $\tilde{I}(\cdot)$ is a linear function of the residuals $y_i-\boldsymbol{\mathbf{\MakeLowercase{x}}}_i'\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}$, just as in the mean IV moment conditions; if $\tau\ne0.5$, then the intercept is different, but the slope estimates are identical to the IV slope estimates for any $\tau$. Of course, in practice $h$ is small, but this shows that the worst-case effect of choosing $h$ too large is that you simply get the usual IV slope estimates.
If there are more instruments than parameters (overidentification), then a different $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ is used in (ref). Specifically, it is replaced by the linear projection of the regressor vector $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ onto the instruments $\boldsymbol{\mathbf{\MakeLowercase{z}}}$, which has the same dimension as $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ and thus $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$. This linear projection is motivated by the two-stage least squares estimator, which can also be written as the mean IV estimator when replacing $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ with the linear projection of $\boldsymbol{\mathbf{\MakeLowercase{x}}}$ onto $\boldsymbol{\mathbf{\MakeLowercase{z}}}$. Although theoretically more efficient estimators may exist, this estimator remains consistent, reliable, and fast to compute regardless of the degree of overidentification.
It remains to motivate the conditional quantile restriction in (ref) from a causal model, as originally done by ChernozhukovHansen2005. \footnote{See also Chapter 6 of Kaplan2020distnp for an introductory discussion.} Instead of assuming every individual (or firm, or school, etc.)\ has the same coefficients $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}$, imagine each individual has their own coefficient vector $\boldsymbol{\mathbf{\MakeLowercase{b}}}$. This is a “random coefficient” model, meaning $\boldsymbol{\mathbf{\MakeLowercase{b}}}$ can differ by individual the same way that $y$, $\boldsymbol{\mathbf{\MakeLowercase{x}}}$, and $\boldsymbol{\mathbf{\MakeLowercase{z}}}$ do:
An additional error term would be redundant; for example, if $y=b_0+b_1x+v$, then the random intercept can simply absorb $v$ to become $b_0+v$. To make the model more tractable empirically, imagine $\boldsymbol{\mathbf{\MakeLowercase{b}}}$ can be written as a deterministic (but unknown) vector-valued function $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\cdot)$ applied to a scalar unobserved $u$:
Assuming $u$ is continuous, it can be normalized to have a $\mathrm{Unif}(0,1)$ distribution (uniformly distributed between $0$ and $1$): any transformation is simply absorbed into $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\cdot)$. Assume a “monotonicity” property that given any value of $\boldsymbol{\mathbf{\MakeLowercase{x}}}$, $\boldsymbol{\mathbf{\MakeLowercase{x}}}'\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(u)$ is an increasing function of $u$, and like before assume the instruments are independent of $u$. Then,
That is, for any $u=\tau$, $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\tau)$ is identified by the conditional quantile restriction from (ref), which can be estimated by the \stcmd{sivqr} command using (ref). Given $\tau$, the function $q_\tau(\boldsymbol{\mathbf{\MakeLowercase{x}}})=\boldsymbol{\mathbf{\MakeLowercase{x}}}'\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\tau)$ is also the $\tau$-quantile structural function introduced by ImbensNewey2009.
This section only applies to the plug-in bandwidth. If instead the \stcmd{bandwidth} option is manually specified, then none of the following applies.
\Citet[prop.\ 2, p.\ 117]{KaplanSun2017} provide a theoretical “optimal” smoothing bandwidth. Specifically, it minimizes the asymptotic mean squared error of the smoothed estimating equations (moment conditions) themselves. It also minimizes the asymptotic mean squared error of a particular linear combination of the estimated coefficients, although it does not do so for every possible linear combination; see section 5 of KaplanSun2017.
A plug-in (data-dependent) version of the optimal bandwidth from KaplanSun2017 is implemented in \stcmd{sivqr} as follows. First, as in their proposition 2, the bandwidth is simplified by assuming $v\equiv y-\boldsymbol{\mathbf{\MakeLowercase{x}}}'\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\tau)$ is independent of the full vector of instruments. Second, the smoothed indicator function is piecewise linear, as mentioned in section (ref). In the notation of KaplanSun2017,
Third, the smallest amount of underlying smoothness in the data-generating process is assumed, $r=2$ KaplanSun2017.
With the chosen $G(\cdot)$ in (ref) and $r=2$, simplifying the optimal bandwidth $h^*$ from propositon 2 of KaplanSun2017 with $v\protect\mathpalette{\protect\independenT}{\perp}\boldsymbol{\mathbf{\MakeLowercase{z}}}$,
where $f_v(\cdot)$ is the probability density function of $v$ and $f_v''(\cdot)$ its second derivative, and $d$ is the dimension (length) of $\boldsymbol{\mathbf{\MakeLowercase{\beta}}}(\tau)$.
The \stcmd{sivqr} plug-in bandwidth is the minimum of a few different values. By taking the minimum (instead of average or maximum), mistakes tend to be in the direction of undersmoothing compared to the “optimal” smoothing that minimizes mean squared error. This results in lower bias (but higher variance) than is optimal, because smoothing reduces variance at the cost of increased bias. The minimum is taken among a nonparametrically estimated version of $h^*$, a Gaussian reference version, and the Silverman rule of thumb bandwidth Silverman1986 for estimating the density of $v$, as suggested by FernandesGuerreHorta2021 for non-IV smoothed QR. The former two are detailed further below.
After finally using the plug-in bandwidth to compute the smoothed IVQR estimator, the residuals $\sthat{v}_i$ are then recomputed from the new $\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}(\tau)$, leading to an updated plug-in bandwidth, and the resulting IVQR estimate is reported by \stcmd{sivqr}.
One approach is to replace $f_v(0)$ and $f_v'(0)$ in $h^*$ by nonparametric kernel estimates. Since $v$ is unobserved, an initial estimate $\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}(\tau)$ must be used to construct residuals
For $f(0)$ (dropping the $v$ subscript for now), given the $\sthat{v}_i$, the usual kernel density estimator is
where $s$ is the kernel bandwidth (not to be confused with the IVQR smoothing bandwidth $h$) and the kernel function $K(\cdot)$ is chosen to be the Gaussian kernel. For $s$, a pointwise version of Silverman's rule of thumb Silverman1986 is used because interest is in $f(0)$, not the full function $f(\cdot)$, so it is better to minimize the pointwise mean squared error than the integrated mean squared error. The (asymptotic) pointwise mean squared error and optimal bandwidth can be found in DasGupta2008, for example:
Following the Gaussian reference approach of Silverman1986, assume $f(\cdot)$ is the density of a $\mathrm{N}(\mu,\sigma^2)$ distribution; that is, $f(x)=\phi((x-\mu)/\sigma)/\sigma$, and the corresponding cumulative distribution function is $F(x)=\Phi((x-\mu)/\sigma)$, where $\phi(\cdot)$ and $\Phi(\cdot)$ are respectively the density and distribution functions of the standard normal distribution. Given the restriction that the $\tau$-quantile of $v$ is zero,
Thus, the density at zero is
Similarly, assuming normality allows us to express $\sthat{f}''(0)$ in terms of $\sthat\sigma$. The second derivative of the $\mathrm{N}(\mu,\sigma^2)$ density evaluated at zero is \footnote{D[PDF[NormalDistribution[\textbackslash[Mu], \textbackslash[Sigma]], x], \{x,2\}] at \url{http://www.wolframalpha.com} for example; then plug in $x=0$ and $-\mu/\sigma=\Phi^{-1}(\tau)$.}
Thus,
Also, for the Gaussian kernel,\footnote{Unfortunately, the value in table 32.1 in DasGupta2008 is incorrectly stated as $1/\sqrt{2\pi}$.}
Finally, plugging everything into the formula for $s^*$,
Following Silverman1986, to get a feasible bandwidth $\sthat{s}$, the unknown $\sigma$ is replaced by either the sample standard deviation of the $\sthat{v}_i$ or their sample interquartile range divided by $1.349$ (the standard normal interquartile range), whichever is smaller. This $\sthat{s}$ is then used to compute $\sthat{f}_v(0)$ using (ref).
For the density derivative $f_v'(0)$ (here simply $f'(0)$), the approach is the same: use a nonparametric kernel estimator with a Silverman-type bandwidth. The usual estimator is
where now $b$ is the bandwidth (to avoid confusion with $h$ and $s$). The asymptotic mean squared error for this estimator at the point zero is from WandJones1994, \footnote{Unfortunately, their (2.34) in section 2.12 has a typo, but the correct version of (2.34) can be derived from the bias and variance expressions in their exercise 2.6(a) on pages 52--53.}
where $\int_\mathbb{R} K(v)v^2\mathop{}\!\mathrm{d}{v}=1$ as in (ref). Solving the first-order condition (setting the derivative with respect to $b$ equal to zero), this leads to the optimal bandwidth
Plugging in the standard normal density for $K(\cdot)$,
Plugging in $f(0)$ from (ref) and \footnote{D[PDF[NormalDistribution[\textbackslash[Mu], \textbackslash[Sigma]], x], \{x,3\}] at \url{http://www.wolframalpha.com} for example; then plug in $x=0$ and $-\mu/\sigma=\Phi^{-1}(\tau)$.}
yields
Again, to get a feasible bandwidth $\sthat{b}$, $\sigma$ is replaced by the sample standard deviation or interquartile range divided by $1.349$, whichever is smaller. This $\sthat{b}$ is then used to compute $\sthat{f}_v'(0)$ using (ref).
Finally, $\sthat{f}_v(0)$ and $\sthat{f}_v'(0)$ are plugged into the expression for $h^*$ to get the feasible plug-in bandwidth $\sthat{h}$.
Here, the Gaussian reference approach is used directly for $h^*$, to simplify the density and density derivative of $v$.
From (ref), assuming normality of $v$ yields $f(0)=\sigma^{-1}\phi(\Phi^{-1}(\tau))$.
Assuming normality, the first derivative of the density of $v$ evaluated at zero is \footnote{D[PDF[NormalDistribution[\textbackslash[Mu], \textbackslash[Sigma]], x], \{x,1\}] at \url{http://www.wolframalpha.com} for example; then plug in $x=0$ and $-\mu/\sigma=\Phi^{-1}(\tau)$.}
Plugging this into the expression for $h^*$ yields
and as usual $\sigma$ is replaced by the smaller of the sample standard deviation of the $\sthat{v}_i$ or the sample interquartile range of the $\sthat{v}_i$ divided by $1.349$.
By default or with \stcmd{reps(0)}, analytic standard errors are computed as follows; otherwise, bootstrap is used as in section (ref). The formulas here are special cases of (3.11) of ChernozhukovHansen2006; (18), (23), and Theorem 5 of deCastroGalvaoKaplanLiu2019; or Theorem 7 of KaplanSun2017. For notational simplicity, I leave the dependence on $\tau$ implicit. Define an additive error term
The first-order asymptotic distribution of the \stcmd{sivqr} estimator is
where $f_{\epsilon|Z,X}(0\mid\boldsymbol{\mathbf{\MakeLowercase{z}}},\boldsymbol{\mathbf{\MakeLowercase{x}}})$ is the conditional PDF of $\epsilon$ given $(\boldsymbol{\mathbf{\MakeLowercase{z}}},\boldsymbol{\mathbf{\MakeLowercase{x}}})$ evaluated at $\epsilon=0$. The expression $\boldsymbol{\mathbf{\MakeUppercase{J}}}^{-1}\boldsymbol{\mathbf{\MakeUppercase{S}}}(\boldsymbol{\mathbf{\MakeUppercase{J}}}')^{-1}$ is mathematically equivalent but requires inverting the asymmetric matrix $\boldsymbol{\mathbf{\MakeUppercase{J}}}$ by itself. The estimated asymptotic covariance matrix is
using a conventional kernel density estimator for $\sthat{\boldsymbol{\mathbf{\MakeUppercase{J}}}}$ with a Gaussian kernel $\phi(\cdot)$ and bandwidth $h$ being Silverman's rule of thumb Silverman1986, specifically $1.06n^{-1/5}$ times the minimum of either the sample standard deviation of the $\boldsymbol{\mathbf{\MakeLowercase{y}}}_i-\boldsymbol{\mathbf{\MakeLowercase{x}}}_i'\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}$ or their sample interquartile range divided by $1.349$. Alternatively, the bandwidth in (A.8) of Liu2019 would also work well. Note $\sthat{\boldsymbol{\mathbf{\MakeUppercase{J}}}}$ is a Powell1986 type estimator, similar to (3.14) of ChernozhukovHansen2006 but here with a Gaussian instead of uniform kernel. The reported standard errors are then the square roots of the diagonals of the matrix $(\sthat{\boldsymbol{\mathbf{\MakeUppercase{J}}}}{}'\sthat{\boldsymbol{\mathbf{\MakeUppercase{S}}}}{}^{-1}\sthat{\boldsymbol{\mathbf{\MakeUppercase{J}}}})^{-1}/n$.
When \stcmd{reps(\num)} specifies multiple bootstrap replications, standard errors are computed by Bayesian bootstrap Rubin1981, assuming i.i.d.\ sampling. (With non-i.i.d.\ sampling, \stcmd{bootstrap} can be used; see \rref{bootstrap}.) This is a particular type of frequentist exchangeable bootstrap vanderVaartWellner1996 that also has a nonparametric Bayesian interpretation ChamberlainImbens2003. In each replication, weights are set as $w_i=\xi_i/\stbar{\xi}$ where the $\xi_i$ are i.i.d.\ standard exponential random variables and $\stbar{\xi}$ is their average; equivalently, the vector $(w_1,\ldots,w_n)/n$ follows a Dirichlet distribution with every parameter equal to one. Using these $w_i$, the estimator is computed by solving (ref). This is repeated many times; the reported standard error is the standard deviation of the different $\sthat{\boldsymbol{\mathbf{\MakeLowercase{\beta}}}}$.
The \stcmd{sivqr} command offers smoothed estimation of the instrumental variables quantile regression model of ChernozhukovHansen2005. The smoothing improves computation and accuracy, helping \stcmd{sivqr} to overcome several limitations of \stcmd{ivqreg} and \stcmd{ivqreg2}, while complementing commands for alternative quantile models with endogeneity (\stcmd{cqiv} and \stcmd{ivqte}). The \stcmd{sivqr} command can help \Stata\ users reliably estimate heterogeneous effects across a variety of settings and datasets.
There would be no \stcmd{sivqr} \Stata\ command without the work of Yixiao Sun on KaplanSun2017. Thanks to Xin Liu for research assistance. Thanks to Joao Santos Silva for suggesting the \stcmd{initial()} option. Thanks to Di Liu for several helpful suggestions and encouragement, including \stcmd{e(bwidth_max)} and analytic standard errors. Many thanks to a particularly thorough anonymous reviewer who found my mistakes and helped improve both the content and presentation tremendously.
\ifnum 35=1 \else \fi