IPARCOS-UCM-23-019

# Operational realization of quantum vacuum ambiguities

Álvaro Álvarez-Domínguez<EMAIL_ADDRESS>Departamento de Física Teórica and
IPARCOS, Universidad Complutense de Madrid, Plaza de las Ciencias 1, 28040
Madrid, Spain  José A.R. Cembranos<EMAIL_ADDRESS>Departamento de Física
Teórica and IPARCOS, Universidad Complutense de Madrid, Plaza de las Ciencias
1, 28040 Madrid, Spain  Luis J. Garay<EMAIL_ADDRESS>Departamento de
Física Teórica and IPARCOS, Universidad Complutense de Madrid, Plaza de las
Ciencias 1, 28040 Madrid, Spain  Mercedes Martín-Benito
<EMAIL_ADDRESS>Departamento de Física Teórica and IPARCOS, Universidad
Complutense de Madrid, Plaza de las Ciencias 1, 28040 Madrid, Spain  Álvaro
Parra-López<EMAIL_ADDRESS>Departamento de Física Teórica and IPARCOS,
Universidad Complutense de Madrid, Plaza de las Ciencias 1, 28040 Madrid,
Spain  Jose M. Sánchez Velázquez<EMAIL_ADDRESS>Instituto de
Física Teórica UAM/CSIC, c/ Nicolás Cabrera 13-15, Cantoblanco, 28049 Madrid,
Spain

###### Abstract

We provide a reinterpretation of the quantum vacuum ambiguities that one
encounters when studying particle creation phenomena due to an external and
time-dependent agent. We propose a measurement-motivated understanding: Each
way of measuring the number of created particles selects a particular vacuum.
This point of view gives a clear and physical meaning to the time evolution of
the number of particles produced by the agent as the counts in a specific
detector, and at the same time relates commonly used quantization
prescriptions to particular measurement setups.

## Introduction.–

When quantizing a matter field in the presence of a classical, external, and
time-dependent agent, there arises the phenomenon of particle production.
Examples of this are gravitational particle creation caused by a time-
dependent geometry Parker (1969); Ford (1987); Weinberg (2008), or pair
production in the Schwinger effect due to an external electric field Sauter
(1931); Schwinger (1951).

These phenomena are usually studied within the framework of quantum field
theory in curved spacetimes Birrell and Davies (1982); Mukhanov and Winitzki
(2007); Calzetta and Hu (2008); Parker and Toms (2009). In the canonical
quantization of the matter fields, one typically requires the resulting
quantum theory to preserve the symmetries of the classical theory. In the
well-known case of flat spacetime devoid of any external background field,
this criterion turns out to select a unique quantization respecting Poincaré
symmetry and the corresponding so-called Minkowski vacuum. However, when we
introduce a time-dependent, external agent such as an electric field, time-
translational invariance is broken. Then, the symmetry group of the classical
theory is not large enough to completely determine a preferred quantum theory.
Indeed, one encounters ambiguities in the choice of annihilation and creation
operators, which lead to different quantum theories, each with their
respective vacuum and particle notions. Even in vacuum, the action of the
external field results in a non-vanishing number $\mathcal{N}(\tau)$ of
particles produced at time $\tau$. This can be understood from the fact that
the notion of particle has changed from the initial time to the time $\tau$.
The precise number of particles, however, strongly depends on the specific
vacuum states that we are comparing. This clearly poses questions about the
physical interpretation of $\mathcal{N}(\tau)$, and the discussion in the
literature is still open Ilderton (2022); Álvarez-Domínguez  _et al._ (2022);
Dabrowski and Dunne (2014, 2016); Yamada (2021); Domcke  _et al._ (2022).

Recent works have experimentally implemented gravitational particle production
in black-hole Steinhauer (2016); Muñoz de Nova  _et al._ (2019) and
cosmological Eckel  _et al._ (2018); Wittemer  _et al._ (2019); Steinhauer
_et al._ (2022); Viermann  _et al._ (2022) analogue systems, where by means of
two-point correlation functions of the density contrast, the number of
produced particles after the expansion was measured. Motivated by the
experimental accessibility of this quantity, we provide a way of understanding
the physical meaning of the possible definitions of $\mathcal{N}(\tau)$ in
terms of the number of particles measured well after the time $\tau$ at which
the external agent has been switched off.

We consider a simple setup in which an electric field is switched on smoothly
from zero, so that there is no ambiguity in the choice of initial vacuum. In
order to measure the actual number of particles at a certain time $\tau_{1}$,
one would need to instantaneously disconnect the external agent, here the
background electric field, and measure afterwards. However, instantaneous
processes are unfeasible, and thus we cannot have experimental access to that
magnitude. Instead, one possibility is to start switching the electric field
off smoothly at that time, wait some time until the electric field is
completely switched off, and finally measure. We denote this outcome as
$\mathcal{N}^{\text{exp}}_{\tau_{1}}$. In order to measure the number of
particles at a later time $\tau_{2}$, we would need to repeat the experiment
switching the external field off at that new instant. In this way, we would
obtain a set of measurement results
$\mathcal{N}^{\text{exp}}_{\tau_{1}},\mathcal{N}^{\text{exp}}_{\tau_{2}},\cdots$,
which tells us what is the number of particles measured in our experiment if
we start to switch off the electric field at $\tau_{1},\tau_{2},\cdots$.
Nevertheless, as we remarked before, this procedure and the results of the
measurement will depend on how we switch the electric field off, which might
be conditioned by the particular characteristics of the detector that we are
using.

Here we propose to relate the different ways in which we can switch the
electric field off and measure the number of particles on the one hand, with
the theoretical ambiguities in the choice of the quantum vacuum on the other.
For each measurement setup, leading to a family of results
$\\{\mathcal{N}_{\tau_{i}}^{\text{exp}}\\}$, we can find among all possible
quantizations a notion of vacuum such that
$\mathcal{N}(\tau_{i})=\mathcal{N}^{\text{exp}}_{\tau_{i}}$. The meaning of
$\mathcal{N}(\tau)$ becomes clear in this case: It is the resulting number of
particles that would be measured, following our particular measurement
process, if we switched off our experiment at time $\tau$. Thus, canonical
quantum ambiguities are inherently physical in the sense that they are
intimately related to the infinitely many different ways of measuring.

## Preliminaries.—

As a working case, let us consider a charged scalar field $\phi(t,x)$ in
(1+1)-dimensional Minkowski spacetime in the presence of a spatially-
homogeneous time-dependent classical electric field (although our analysis can
be extrapolated to higher dimensions and to similar particle creation
scenarios due to an external time-dependent agent). Its dynamics is determined
by the Klein-Gordon equation

$\left[\left(\partial_{\mu}+iqA_{\mu}\right)\left(\partial^{\mu}+iqA^{\mu}\right)+m^{2}\right]\phi(t,x)=0,$
(1)

where $m$ and $q$ are the mass and the charge of the field, respectively, and
$A_{\mu}$ is the external electric potential. We choose the temporal gauge,
namely $A_{\mu}=(0,A(t))$, and thus $E(t)=-\dot{A}(t)$, so that we explicitly
preserve the spatial homogeneity of the electric field $E(t)$ in the
quantization of the scalar field.

In order to construct a quantum field operator, we expand the matter field in
Fourier time-dependent modes $\phi_{k}$ which verify decoupled harmonic
oscillator equations

$\ddot{\phi}_{k}(t)+\Omega_{k}(t)^{2}\phi_{k}(t)=0$ (2)

with time-dependent frequencies

$\Omega_{k}(t)^{2}=k^{2}+2qA(t)k+q^{2}A(t)^{2}+m^{2}.$ (3)

Both the real and imaginary parts of the Fourier modes $\phi_{k}$ satisfy the
same equation. So, from now on and for simplicity in the notation, the symbol
$\phi_{k}$ will be a real variable that represents each of them. The final
particle number will simply require an extra factor of 2 to account for this
duplicity. The expansion of each (real) mode in terms of a basis of solutions
$\varphi_{k}$ of Eq. (2),

${\phi}_{k}(t)={a}_{k}\varphi_{k}(t)+{a}^{*}_{k}\varphi_{k}^{*}(t),$ (4)

determines the Fock quantization associated with this choice of basis,
provided that the expansion coefficients are promoted to annihilation and
creation operators $\hat{a}_{k}$ and $\hat{a}^{\dagger}_{k}$, respectively.
The corresponding quantum vacuum $\ket{0}$ satisfies, by definition,
$\hat{a}_{k}\ket{0}=0$ for all $k$. Note that modes with different wavenumbers
$k$ are decoupled from one another. Thus, in the following we will drop the
index $k$, recovering it explicitly whenever is necessary.

We want to account for realistic (non-instantaneous) switch-ons and switch-
offs of the electric field (see Fig. 1). We start the experiment switching the
electric field on at a time $t_{\text{on}}$, and after a time
$\delta_{\text{on}}$ the electric field smoothly reaches the constant value
$E_{0}$. Then, in order to measure the number of created particles at a given
time $\tau$, we start switching it off at that time, and after a lapse
$\delta$, the electric field vanishes. This process is characterized by the
properties of the experimental setup. The finite duration of the on and off
processes can be described by a $C^{\infty}$ step function
$\Theta_{\sigma}(t)$ of width $\sigma$ (i.e., $\delta_{\text{on}}$ and
$\delta$, respectively, in our case at hand). The step Heaviside function
corresponds to $\sigma=0$. The specific form of $\Theta_{\sigma}(t)$ does not
qualitatively affect the results. In our numerical computation, we use a
$C^{\infty}$ regularization which interpolates between 0 and 1 in the interval
$(-\sigma/2,\sigma/2)$ by means of the function
$\bm{(}1+\tanh\\{\cot[\pi(1/2-t/\sigma)]\\}\bm{)}/2$ and is constant outside
said interval. In all figures we fix $m=1$ and $qE_{0}=1$, so that the
electric field reaches the critical Schwinger limit $m^{2}/q$ Schwinger
(1951); Yakimenko  _et al._ (2019). For lower strengths, the probability of
producing a particle pair is negligible. In addition, we set $t_{\text{on}}=0$
and times are given in units of the switch-on duration $\delta_{\text{on}}$,
while we parametrize different switch-offs by varying $\delta$.

Figure 1: Time evolution of the electric field (solid line) with different
switch-off profiles (dashed/dotted lines) corresponding to different values of
$\delta$, starting at $\tau$.

## Measured particle number.—

Given a particular experimental setting, we can compute the asymptotic number
of created particles $N^{\text{exp}}_{\tau}$ in the mode $k$ that would be
measured by our detector when we start the switch-off at the time $\tau$.
Initially, when the electric field is not on yet, the matter field is in the
Minkowski vacuum. This state is determined by the solution
$\varphi^{\text{in}}$ to the mode equation (2) which behaves as a positive-
frequency plane wave before $t_{\text{on}}$. This is the only solution
compatible with Poincaré symmetry in the asymptotic past. Analogously, there
is another Minkowski out-vacuum, associated with $\varphi^{\text{out}}$, which
is a positive-frequency plane wave at times after the electric field is
completely switched off, namely, at $t_{\text{off}}=\tau+\delta$. However,
because of the presence of the electric field, these two vacua are different.
Indeed, $N_{\tau}^{\text{exp}}$ measures how excited is the in-vacuum with
respect to the out-vacuum. To obtain this quantity, we need to compare both
solutions to Eq. (2) at the same time. Explicitly, we evolve the in-solution
from the initial time and compare it with $\varphi^{\text{out}}$, both
evaluated at $t_{\text{off}}$. We can then calculate the asymptotic number of
created particles in the mode $k$ according to the Bogoliubov formalism as

$N_{\tau}^{\text{exp}}=\left|\varphi^{\text{in}}(t_{\text{off}})\dot{\varphi}^{\text{out}}(t_{\text{off}})-\dot{\varphi}^{\text{in}}(t_{\text{off}})\varphi^{\text{out}}(t_{\text{off}})\right|^{2}.$
(5)

Figure 2: Spectra of the asymptotic number of created particles
$N_{\tau}^{\text{exp}}$ for different switch-offs, starting at times
$\tau=2.5$ (left) and $\tau=3$ (right) for different switch-off durations
$\delta$ (dashed lines). We also represent the computed number of created
particles $N^{\text{ad}}(\tau)$ in the zeroth-order adiabatic vacuum at those
times (solid line).

It is important to remark that for a different measurement process starting at
the same $\tau$, here characterized by another duration of the switch-off
$\delta$, the outcome $N_{\tau}^{\text{exp}}$ will change. In Fig. 2, we show
the spectra of asymptotically produced particles $N_{\tau}^{\text{exp}}$ at
$\tau$ for different switch-off durations $\delta$. Observe that each spectrum
oscillates in $k$, and increasing the value of $\delta$ shifts the phase in a
way such that fluctuations are displaced towards the infrared. Note that these
oscillations in the spectral distribution have already been observed in analog
experiments Steinhauer  _et al._ (2022); Viermann  _et al._ (2022).

## Theoretical particle number.—

On the other hand, one may be interested in theoretically computing the number
of particles that have been created from an initial time $t_{\text{on}}$ to
some time $\tau$, without relying on any measurement. This question is,
however, much more subtle. One would need to choose a particular solution
$\varphi^{\tau}$ to Eq. (2) by imposing initial conditions at time $\tau$,
i.e., $(\varphi^{\tau}(\tau),\dot{\varphi}^{\tau}(\tau))$. Contrarily to what
happens for the out-solution $\varphi^{\text{out}}$, there is an ambiguity in
the selection of $\varphi^{\tau}$. Indeed, when the electric field is on, the
frequency (3) is not constant and the physical criterion of preservation of
the classical symmetries in the quantum theory is not strong enough to fix a
unique vacuum. The spectral number of created particles between
$t_{\text{on}}$ and $\tau$ is defined by

$N(\tau)=\left|\varphi^{\text{in}}(\tau)\dot{\varphi}^{\tau}(\tau)-\dot{\varphi}^{\text{in}}(\tau)\varphi^{\tau}(\tau)\right|^{2}.$
(6)

In contrast with the measured value $N_{\tau}^{\text{exp}}$, which is unique
given a particular switch-off of the electric field starting at a certain time
$\tau$, $N(\tau)$ strongly depends on the choice of vacuum at that time. This
is precisely the reason why the interpretation of this magnitude is not clear
in the literature yet Ilderton (2022); Álvarez-Domínguez  _et al._ (2022);
Dabrowski and Dunne (2014, 2016); Yamada (2021).

As a well-known example, the computed number of created particles
$N^{\text{ad}}(\tau)$ for the zeroth-order adiabatic initial conditions at
$\tau$ Birrell and Davies (1982) is depicted in Fig. 2 by the solid line. Note
that this curve approaches that of $N_{\tau}^{\text{exp}}$ for large values of
$\delta$. We will comment on this later on.

## Relation between measured and theoretical particle numbers.—

Each measurement procedure selects a particular vacuum for which the
theoretical number of particles $N(\tau)$ has a well-defined physical meaning.
Indeed, among all possibilities for choosing a particular mode
$\varphi^{\tau}$ in Eq. (6), there is one for which the number of predicted
particles at time $\tau$ coincides with the outcome $N^{\text{exp}}_{\tau}$
that a particular measurement device would yield. More explicitly, for each
time $\tau$, we choose the mode $\varphi^{\tau}$ as the out-vacuum
$\varphi^{\text{out}}$ associated with the switch-off starting at that time.
In fact, replacing $t_{\text{off}}$ by $\tau$ in Eq. (5) does not change the
resulting $N^{\text{exp}}_{\tau}$, since this magnitude does not depend on the
instant at which it is evaluated. Consequently, the choice
$\varphi^{\tau}=\varphi^{\text{out}}$ makes Eq. (6) equal to Eq. (5). The set
of modes $\\{\varphi^{\text{exp},\tau}\\}$ and corresponding vacua
$\\{\ket{0}^{\text{exp},\tau}\\}$ defined in this way, one for each $\tau$,
allows one to construct $N^{\text{exp}}(\tau)$, viewed as function of the time
at which we start switching the electric field off. At each time $\tau$, the
in-vacuum is an excited state with respect to the vacuum
$\ket{0}^{\text{exp},\tau}$: Its excitations correspond precisely to the
quanta which would be measured by our detector.

This prescription defines a family of physical vacua: Those for which there
exists a switch-off giving the same particle number as the one predicted by
the vacuum. Furthermore, all these vacua unitarily implement the dynamics, as
they are associated with a finite number of particles by construction, thus
verifying Shale’s theorem Shale (1962).

In our simple setup, the measurement device is characterized by $\delta$, and
for each of its values we have a different set of modes
$\\{\varphi^{\text{exp},\tau}\\}$ and hence different notions of
$N^{\text{exp}}(\tau)$. This is illustrated in Fig. 3, where one can see the
time evolution of $N^{\text{exp}}(\tau)$ for a certain infrared mode, for
different durations of the switch-off $\delta$. For each time $\tau$, we
compute the asymptotic number of particles $N_{\tau}^{\text{exp}}$ when we
start switching the electric field off at $\tau$. The observed oscillations in
$\tau$ were already present in references Kluger  _et al._ (1998); Schmidt
_et al._ (1998); Dabrowski and Dunne (2014, 2016), but it is now when we can
provide them with a full physical meaning, as they follow from a measurement-
based notion of particle. Moreover, recent works try to implement experimental
setups that make use of this behavior to enhance particle production (see the
recent study Aleksandrov  _et al._ or other references on the dynamically
assisted Schwinger effect Schützhold  _et al._ (2008); Bulanov  _et al._
(2010)).

The amplitudes of the fluctuations are smaller as we increase the value of
$\delta$. In agreement with Fig. 2, the evolution of $N^{\text{exp}}(\tau)$
for long switch-offs is well approximated by the prediction obtained with the
zeroth-order adiabatic prescription at time $\tau$, to which we dedicate the
next section.

Figure 3: Evolution of the number of created particles $N^{\text{exp}}(\tau)$
with $k=0.5$ for different switch-off durations $\delta$. The solid line
corresponds to the zeroth-order adiabatic prescription.

As an aside, for each measurement process we could have reassigned the
asymptotic outcome $N_{\tau}^{\text{exp}}$ to any other time different to the
time at which the switch-off starts (e.g., results in Ilderton (2022) are
obtained by taking $\tau+\delta/2$ as the reference time instead). However,
note that this would lead to a simple relabelling of the $\tau$-axis in Fig.
3, shifting each curve proportionally to its value of $\delta$.

Finally, one may wonder whether it is possible to find, given a choice of
vacuum, a particular switch-off starting at $\tau$ such that the measurement
of the associated asymptotic number of particles $N_{\tau}^{\text{exp}}$
coincides with the value of $N(\tau)$, computed using said vacuum
prescription. This requirement does not unequivocally determine the time
evolution of the in-mode $\varphi^{\text{in}}$ after time $\tau$. Therefore,
each function $\varphi^{\text{in}}$ compatible with the previous condition
would lead to a different mode equation. However, it is definitely non-trivial
that one can find a time-dependent frequency $\Omega(t)$ of the form of Eq.
(3) fulfilling this requirement for all values of $k$.

## Zeroth-order adiabatic vacuum.—

Suppose the unfeasible situation in which the electric field is
instantaneously switched off at a particular time $\tau$, i.e., $\delta=0$.
This can be easily implemented in the electric potential with a continuous but
non-differentiable step function at $\tau$. For the out-solution to be smooth,
its form at $\tau$ has to be that of a plane wave with frequency
$\Omega(\tau)$ as given in Eq. (3). Then, this out-solution coincides with the
zeroth-order adiabatic vacuum at time $\tau$. This is precisely the case
analyzed in Ilderton (2022), and the particular number of particles
$N^{\text{ad}}(\tau)$ calculated using this well-known adiabatic vacuum is
indeed the one used in many references, especially those studying the Quantum
Vlasov Equation Kluger  _et al._ (1998); Schmidt  _et al._ (1998); Aleksandrov
_et al._ (2022); Anderson and Mottola (2014); Ruffini  _et al._ (2010);
Roberts and Schmidt (2000); Dunne (2009); Dumlu and Dunne (2011); Hebenstreit
_et al._ (2009).

Figure 4: The solid lines represent the total measured number of created
particles as a function of the switch-off duration $\delta$, for two different
starting times $\tau=2.5$ (blue) and $\tau=3$ (red). The dashed lines are the
integrated number of particles at each $\tau$, obtained by imposing the
zeroth-order adiabatic prescription.

We infer from Figs. 2 and 3 that the prediction of the adiabatic vacuum is the
opposite of that of a very fast but non-instantaneous switch-off (in contrast
with Ilderton (2022)). Furthermore, in Fig. 4 we consider smooth switch-offs
of various durations $\delta$ and compute the integrated number of created
particles $\mathcal{N}^{\text{exp}}(\tau)=\int dk\ N_{k}^{\text{exp}}(\tau)$
which would be measured in the experiment, at two different times. We identify
two asymptotic regions separated by a transition regime for switch-offs with a
length of the order of the switch-on duration $\delta_{\text{on}}=1$. For
larger values of $\delta$, $\mathcal{N}^{\text{exp}}(\tau)$ approaches
$\mathcal{N}^{\text{ad}}(\tau)=\int dk\ N^{\text{ad}}_{k}(\tau)$, i.e., the
zeroth-order adiabatic result—the limit is not exact because of the difference
at sufficiently infrared scales, as it was shown in Figs. 2 and 3. For fast
switch-offs the measured particle number significantly deviates from the
adiabatic prediction depending on $\tau$, which affects the distribution of
maxima and minima of the particle spectrum (see Fig. 2).

To summarize, the measurement process provides a physical meaning to the
quanta predicted by the ubiquitous zeroth-order adiabatic vacuum: They
correspond to those measured by a sufficiently long switch-off.

## Conclusions.—

Quantum vacuum ambiguities are inherent to quantum field theory in the
presence of an external, time-dependent agent. In this work, we show that
knowing the particularities of how we measure the particle number allows us to
identify a particular quantum vacuum with clear physical meaning: Its
associated notion of particle is that which would be measured by our detector
in a potential experiment.

As an example, we show how the zeroth-order adiabatic vacuum prescription
approximates the result corresponding to a slow switch-off. Definitely, it
does not correspond to the limit when the duration of the switch-off goes to
zero.

In conclusion, we select a family of vacua—those related to a realistic
switch-off of the external agent—which are physical in the sense that they
accommodate information about real outcomes.

###### Acknowledgements.

## Acknowledgments.—

This work is partially supported by the MICINN (Spain) projects
PID2019-107394GB-I00/AEI/10.13039/501100011033 (AEI/FEDER, UE),
PID2019-108655GB-I00 (AEI) and PID2020-118159GB-C44, the EU STRONG-2020
project [grant no. 824093], COST (European Cooperation in Science and
Technology) Actions CA21106 and CA21136, and STMS Grant from COST Action
CA16108. Á.P.-L. acknowledges financial support from the MIU (Ministerio de
Universidades, Spain) fellowship FPU20/05603. JMSV acknowledges the support of
the Spanish Agencia Estatal de Investigación through the grant “IFT Centro de
Excelencia Severo Ochoa CEX2020-001007-S".

## References

  * Parker (1969) L. Parker, Phys. Rev. 183, 1057 (1969).
  * Ford (1987) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * Weinberg (2008) S. Weinberg,  _Cosmology_ (Oxford University Press, Oxford, 2008).
  * Sauter (1931) F. Sauter, Z. Phys. 69, 742 (1931).
  * Schwinger (1951) J. S. Schwinger, Phys. Rev. 82, 664 (1951).
  * Birrell and Davies (1982) N. D. Birrell and P. C. W. Davies,  _Quantum Fields in Curved Space_ (Cambridge University Press, Cambridge, 1982).
  * Mukhanov and Winitzki (2007) V. Mukhanov and S. Winitzki,  _Introduction to Quantum Effects in Gravity_ (Cambridge University Press, Cambridge, 2007).
  * Calzetta and Hu (2008) E. A. Calzetta and B.-L. B. Hu,  _Nonequilibrium Quantum Field Theory_ (Cambridge University Press, Cambridge, 2008).
  * Parker and Toms (2009) L. E. Parker and D. Toms,  _Quantum Field Theory in Curved Spacetime: Quantized Field and Gravity_ (Cambridge University Press, Cambridge, 2009).
  * Ilderton (2022) A. Ilderton, Phys. Rev. D 105, 016021 (2022).
  * Álvarez-Domínguez  _et al._ (2022) A. Álvarez-Domínguez, L. J. Garay, and M. Martín-Benito, Phys. Rev. D 105, 125012 (2022).
  * Dabrowski and Dunne (2014) R. Dabrowski and G. V. Dunne, Phys. Rev. D 90, 025021 (2014).
  * Dabrowski and Dunne (2016) R. Dabrowski and G. V. Dunne, Phys. Rev. D 94, 065005 (2016).
  * Yamada (2021) Y. Yamada, JCAP 2021, 009 (2021) (2021).
  * Domcke  _et al._ (2022) V. Domcke, Y. Ema, and K. Mukaida, JHEP 2022, 33 (2022).
  * Steinhauer (2016) J. Steinhauer, Nat. Phys. 12, 959 (2016).
  * Muñoz de Nova  _et al._ (2019) J. R. Muñoz de Nova, K. Golubkov, V. I. Kolobov, and J. Steinhauer, Nature 569, 688 (2019).
  * Eckel  _et al._ (2018) S. Eckel, A. Kumar, T. Jacobson, I. B. Spielman, and G. K. Campbell, Phys. Rev. X 8, 021021 (2018).
  * Wittemer  _et al._ (2019) M. Wittemer, F. Hakelberg, P. Kiefer, J.-P. Schröder, C. Fey, R. Schützhold, U. Warring, and T. Schaetz, Phys. Rev. Lett. 123, 180502 (2019).
  * Steinhauer  _et al._ (2022) J. Steinhauer, M. Abuzarli, T. Aladjidi, T. Bienaimé, C. Piekarski, W. Liu, E. Giacobino, A. Bramati, and Q. Glorieux, Nat. Commun. 13, 2890 (2022).
  * Viermann  _et al._ (2022) C. Viermann _et al._ , Nature 611, 260 (2022).
  * Yakimenko  _et al._ (2019) V. Yakimenko _et al._ , Phys. Rev. Lett. 122, 190404 (2019).
  * Shale (1962) D. Shale, Trans. Am. Math. Soc. 103, 149 (1962).
  * Kluger  _et al._ (1998) Y. Kluger, E. Mottola, and J. M. Eisenberg, Phys. Rev. D 58, 12 (1998).
  * Schmidt  _et al._ (1998) S. M. Schmidt, D. Blaschke, G. Ropke, S. A. Smolyansky, A. V. Prozorkevich, and V. D. Toneev, Int. J. Mod. Phys. E 7, 709 (1998).
  * (26) I. A. Aleksandrov, D. G. Sevostyanov, and V. M. Shabaev, “Schwinger particle production: rapid switch off of the external field versus dynamical assistance,” arXiv:2210.15626 .
  * Schützhold  _et al._ (2008) R. Schützhold, H. Gies, and G. Dunne, Phys. Rev. Lett. 101, 130404 (2008).
  * Bulanov  _et al._ (2010) S. S. Bulanov, V. D. Mur, N. B. Narozhny, J. Nees, and V. S. Popov, Phys. Rev. Lett. 104, 220404 (2010).
  * Aleksandrov  _et al._ (2022) I. A. Aleksandrov, D. G. Sevostyanov, and V. M. Shabaev, Symmetry 14, 2444 (2022).
  * Anderson and Mottola (2014) P. R. Anderson and E. Mottola, Phys. Rev. D 89, 104038 (2014).
  * Ruffini  _et al._ (2010) R. Ruffini, G. Vereshchagin, and S.-S. Xue, Phys. Rep. 487, 1 (2010).
  * Roberts and Schmidt (2000) C. D. Roberts and S. M. Schmidt, Prog. Part. Nucl. Phys 45, 1 (2000).
  * Dunne (2009) G. V. Dunne, Eur. Phys. J. D 55, 327 (2009).
  * Dumlu and Dunne (2011) C. K. Dumlu and G. V. Dunne, Phys. Rev. D 83, 065028 (2011).
  * Hebenstreit  _et al._ (2009) F. Hebenstreit, R. Alkofer, G. V. Dunne, and H. Gies, Phys. Rev. Lett. 102, 150404 (2009).


David M. Kaplan
Department of Economics
University of Missouri
Columbia, MO
<EMAIL_ADDRESS>
# Smoothed instrumental variables quantile regression

D. M. Kaplan

###### Abstract

In this article, I introduce the sivqr command, which estimates the
coefficients of the instrumental variables (IV) quantile regression model
introduced by Chernozhukov and Hansen (2005). The sivqr command offers several
advantages over the existing ivqreg and ivqreg2 commands for estimating this
IV quantile regression model, which complements the alternative “triangular
model” behind cqiv and the “local quantile treatment effect” model of ivqte.
Computationally, sivqr implements the smoothed estimator of Kaplan and Sun
(2017), who show that smoothing improves both computation time and statistical
accuracy. Standard errors are computed analytically or by Bayesian bootstrap;
for non-i.i.d. sampling, sivqr is compatible with bootstrap. I discuss syntax
and the underlying methodology, and I compare sivqr with other commands in an
example.

###### keywords:

st0001, sivqr, endogeneity, instrumental variables, quantile regression

##  1 Introduction

The instrumental variables quantile regression (IVQR) model of Chernozhukov
and Hansen (2005) 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 ivqreg (Kwak, 2010a) and ivqreg2 (Machado and Santos Silva,
2018) have several limitations. First, ivqreg (implementing the estimator of
Chernozhukov and Hansen, 2006) 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,
ivqreg2 (implementing the estimator of Machado and Santos Silva, 2019) imposes
a location–scale model that can help if correct but that can result in
inconsistency if misspecified. Machado and Santos Silva (2019, p. 159) 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, sivqr and ivqreg2 could be
combined using a “model averaging” approach like in Cheng et al. (2019) or Liu
(2019); meanwhile, simply running both may be insightful. Third, the existing
commands’ computation time can be long. Fourth, ivqreg and ivqreg2 share some
minor inconveniences, like not supporting syntax for factor variables and
interaction terms. Fifth, 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 xvar*), but it introduces an error when
there is overidentification (more excluded instruments than endogenous
regressors).

To address these limitations, I introduce the sivqr command. The 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 4.1 example is around $25$ seconds for sivqr, as opposed to a
few minutes (ivqreg) or over $20$ minutes (ivqreg2). Although these
differences become trivial on small datasets where everything computes
quickly, they become even more important when using bootstrap for cluster-
robust standard errors, which requires computing the estimator many times.

The sivqr command implements the smoothed IVQR estimator of Kaplan and Sun
(2017), using familiar syntax. Syntax is similar to the ivregress command (see
[R] ivregress) but additionally specifying a quantile level. Likewise, syntax
is similar to the qreg command (see [R] qreg) but additionally specifying
instruments for the endogenous regressors. With non-i.i.d. sampling, the
bootstrap prefix can be used for standard errors; see [R] bootstrap and the
example in section 4.

Other Stata commands that address endogeneity in quantile models are based on
different models than that of Chernozhukov and Hansen (2005). These
alternative models have neither “stronger” nor “weaker” assumptions; they may
be more plausible in certain settings and less plausible in others. The ivqte
command (Frölich and Melly, 2010) estimates the unconditional (default) or
conditional (option aai) “local quantile treatment effect,” analogous to the
local average treatment effect model of Imbens and Angrist (1994). This
conditional (aai) estimand and the IVQR estimand are different but closely
related, as explained in section 10.5 of Melly and Wüthrich (2017) and
detailed in Wüthrich (2020), both of whom compare the models more generally,
too. See also Frölich and Melly (2013) for details. Practically, compared to
ivqte with the aai option, the main advantage of sivqr is the ability to
handle more than a single, binary endogenous regressor. The cqiv command
(Chernozhukov et al., 2019, 2012) uses a control function estimator based on a
triangular model (Lee, 2007), also allowing for censoring (Chernozhukov et
al., 2015). The estimand is the same as for IVQR (see section 2.1).
Practically, compared to cqiv, the main advantages of 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 Chernozhukov et al. (2017).

Unfortunately, all these commands (including 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 Chernozhukov and Hansen (2008), Chernozhukov et
al. (2009), Chernozhukov et al. (2017, §9.3.3), and references therein. See
also Kaplan and Liu (2021) 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” $F$ statistic for assessing
instrument strength in IVQR.

Section 2 discusses the methodology at a relatively intuitive level. Section 3
describes syntax and usage of sivqr. Section 4 provides examples that can be
replicated with the provided do-file. Section 5 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).

##  2 A gentle introduction to methodology

This section discusses methodology at a relatively non-technical level
(compared to section 5).

###  2.1 Parameter interpretation

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_{1}x+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_{1}x$. 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_{1}x=\beta_{0}+\beta_{1}x+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$ (Chernozhukov and Hansen, 2005). 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 Imbens and Newey (2009, §3.1): 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$.

###  2.2 Estimation

Chernozhukov and Hansen (2005) 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_{1}x$ the standard IV
moment conditions are

$\begin{split}0&=\operatorname{E}(v)=\operatorname{E}(y-\beta_{0}-\beta_{1}x)\\\
0&=\operatorname{E}(zv)=\operatorname{E}[z(y-\beta_{0}-\beta_{1}x)]\end{split}$

The IVQR moment conditions are

$\begin{split}0&=\operatorname{E}[\operatorname{\mathbf{1}}\\{y-\beta_{0}(\tau)-\beta_{1}(\tau)x\leq
0\\}-\tau]\\\
0&=\operatorname{E}[z(\operatorname{\mathbf{1}}\\{y-\beta_{0}(\tau)-\beta_{1}(\tau)x\leq
0\\}-\tau)]\end{split}$ (1)

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_{1}x$ (given $z$).

Smoothing solves the computational difficulties inherent in (1). 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 Kaplan and Sun (2017).

##  3 The sivqr command

The 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 by is allowed (see
[D] by). After running, postestimation commands like [R] test and [R] predict
can be used as usual.

###  3.1 Syntax

sivqr depvar $\bigl{[}$ varlist1 $\bigr{]}$ (varlist2 =
varlist${}_{\textup{iv}}$) $\bigl{[}$ if $\bigr{]}$ $\bigl{[}$ in $\bigr{]}$
$\bigl{[}$ weight $\bigr{]}$ , quantile(#) $\bigl{[}$ bandwidth(#) level(#)
reps(#) logiterations noconstant seed(#) nodots initial(matname) $\bigr{]}$

As in ivregress: varlist1 has exogenous regressors (or control variables),
varlist2 has endogenous regressors, and varlist${}_{\textup{iv}}$ has excluded
instruments (exogenous variables that instrument for varlist2).

pweights, iweights, and fweights are allowed; see [U] 11.1.6 weight.

###  3.2 Options

quantile(#) specifies the quantile level as in [R] qreg: a real number
strictly between $0$ and $1$, or alternatively a number between $1$ and $100$
interpreted as a percentile. For example, quantile(0.5) specifies the median,
as does quantile(50). To prevent inadvertent mistakes, there is no default;
this must be specified explicitly.

bandwidth(#) manually specifies the desired smoothing bandwidth. Leaving it
unspecified invokes a plug-in bandwidth based on Kaplan and Sun (2017) 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
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 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.

level(#) specifies the confidence level, as a percentage, for confidence
intervals. The default is level(95) or as set by set level; see [U] 20.7
Specifying the width of confidence intervals.

reps(#) specifies the number of bootstrap replications for estimating the
variance–covariance matrix and standard errors. With the default reps(0),
heteroskedasticity-robust analytic standard errors are reported. Otherwise,
the Bayesian bootstrap of Rubin (1981) is used; it is a valid frequentist
bootstrap that also has a nonparametric Bayesian interpretation (see section
5.6).

logiterations prints each iteration of the numerical solver; see [M-5]
solvenl( ), in particular solvenl_init_iter_log( ). The default is not to
print such information, which usually only helps debugging and
troubleshooting.

noconstant omits the intercept term that otherwise is included automatically.

seed(#) sets the random-number seed, to make results replicable; default is
seed(112358). The current seed is restored at the end of execution.

nodots suppresses display of the replication dots (see [R] bootstrap).

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 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.

###  3.3 Stored results

sivqr stores the following in e():

Scalars |  |  |
---|---|---|---
e(N) |  number of observations |  e(reps) |  number of bootstrap replications
e(bwidth) |  smoothing bandwidth used |  e(q) |  quantile level requested
e(bwidth_req) |  smoothing bandwidth requested requested (or plug-in value) |  e(bwidth_max) |  maximum plug-in bandwidth

Macros |  |  |
---|---|---|---
e(cmd) |  sivqr |  e(constant) |  noconstant if specified
e(instd) |  instrumented variable(s) |  e(insts) |  instrument(s)
e(depvar) |  name of dependent variable |  e(exogr) |  exogenous regressors
e(wtype) |  weight type |  e(wexp) |  weight expression
e(properties) |  b V |  e(title) |  title in estimation output
e(vcetype) |  Robust or Bootstrap |  |

Matrices |  |  |
---|---|---|---
e(b) |  estimated coefficient vector |  e(V) |  estimated variance–covariance matrix of the estimator

Functions |  |  |
---|---|---|---
e(sample) |  marks estimation sample |  |

Almost all of the above are the same as for the familiar commands ivregress
and/or qreg (or bsqreg). The only exceptions are for the smoothing bandwidth.
By default, the plug-in bandwidth is computed and returned in e(bwidth_req).
More specifically, three plug-in bandwidths are computed (see section 5.4), of
which e(bwidth_req) is the minimum and e(bwidth_max) is the maximum. If the
numerical solver cannot find the solution with 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 e(bwidth). There is
nothing wrong with these values being different; for example, specifying
bandwidth(0) makes e(bwidth_req) zero and e(bwidth) the smallest possible
bandwidth for which the numerical solver finds a solution.

If the default plug-in bandwidth is used and e(bwidth) is greater than
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 ivregress and then estat firststage. If e(bwidth) is above
e(bwidth_req) but less than e(bwidth_max), then it may be simply that the
plug-in bandwidth was chosen too conservatively (see section 5.4), but it may
still be worth running ivregress and estat firststage.

sivqr stores the following in r():

Scalar |  |  |
---|---|---|---
r(level) |  confidence level |  |
Matrix |  |  |
r(table) |  table of results |  |

These are the standard values stored by ereturn display; see [P] ereturn.

##  4 Example: Example 4 from ivregress

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 (StataCorp,
2021) 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 sivqr, I used version 2.6 (29sep2021) of ivqreg2 from SSC
(Machado and Santos Silva, 2018), and version 1.0.0 (28jul2010) of ivqreg
(Kwak, 2010a). The most recent version 1.0.1 (02dec2010) from Kwak (2010b)
unfortunately gives an error when trying to use an overidentified model (as in
section 4.1), “st_numscalar(): 3204 matrix found where scalar required.” Other
changes in version 1.0.1 are both positive and negative; modifying the section
4.1 model to only have 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 ivqreg and ivqreg2, in terms of
estimates, speed, and capabilities, including:

  1. 1.

sivqr with a small bandwidth produces estimates very similar to ivqreg;

  2. 2.

ivqreg2 can differ from ivqreg and sivqr, though is qualitatively similar;

  3. 3.

ivqreg2 can be significantly slower than sivqr.

As noted earlier, ivqreg does not allow multiple endogenous regressors, which
precludes even a quadratic function of the endogenous regressor, and the
factorial operator ## is supported only by sivqr.

###  4.1 Main comparison

The following revisits Example 4 in [R] ivregress. Briefly, it hopes to
estimate a structural wage model, treating job tenure as endogenous; see [R]
ivregress for details. I focus on the coefficient on job tenure.

For comparison, 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.

  *

        . webuse nlswork , clear
    (National Longitudinal Survey of Young Women, 14-24 years old in 1968)

    . ivregress 2sls ln_wage c.age##c.age birth_yr grade (tenure = union wks_work msp),
    >   vce(cluster idcode)

    Instrumental variables 2SLS regression            Number of obs   =     18,625
                                                      Wald chi2(5)    =    1600.81
                                                      Prob > chi2     =     0.0000
                                                      R-squared       =          .
                                                      Root MSE        =     .48441

                                 (Std. err. adjusted for 4,110 clusters in idcode)

                                 Robust
         ln_wage   Coefficient  std. err.      z    P>|z|     [95% conf. interval]

          tenure     .1060832   .0044835    23.66   0.000     .0972957    .1148707
             age     .0162345   .0069467     2.34   0.019     .0026193    .0298497

     c.age#c.age    -.0005309   .0001155    -4.59   0.000    -.0007573   -.0003044

        birth_yr    -.0091139   .0022602    -4.03   0.000    -.0135438   -.0046839
           grade      .070454   .0031158    22.61   0.000     .0643471    .0765609
           _cons     .9079537   .1681148     5.40   0.000     .5784548    1.237453

    Instrumented: tenure
     Instruments: age c.age#c.age birth_yr grade union wks_work msp

The following shows results from sivqr at the median level. To cluster
standard errors at the individual level, the bootstrap prefix is used after
using xtset appropriately.111For example, see
https://www.stata.com/support/faqs/statistics/bootstrap-with-panel-data/ With
bootstrap, it would only waste time (without changing the reported estimates
or standard errors) to specify the reps option (with a non-zero number) for
sivqr. At the median, the sivqr coefficient estimates (other than the
intercept) as well as the corresponding cluster-robust standard errors are
similar to those from ivregress 2sls.

  *

        . generate newid = idcode

    . xtset newid year

    Panel variable: newid (unbalanced)
     Time variable: year, 68 to 88, but with gaps
             Delta: 1 unit

    . bootstrap , reps(100) cluster(idcode) idcluster(newid) seed(112358) : sivqr ln_
    > wage c.age##c.age birth_yr grade (tenure = union wks_work msp) , quantile(0.50)
    (running sivqr on estimation sample)

    Bootstrap replications (100)
              1          2          3          4          5
    ..................................................    50
    ..................................................   100

    Smoothed instrumental variables quantile regression (SIVQR)   Quantile = .5
    Smoothing bandwidth used = .0600669                 Number of obs =     18,625
                                  (Replications based on 4,110 clusters in idcode)

                     Observed   Bootstrap                         Normal-based
         ln_wage   coefficient  std. err.      z    P>|z|     [95% conf. interval]

          tenure     .1076941   .0046079    23.37   0.000     .0986627    .1167255
             age     .0060803   .0073372     0.83   0.407    -.0083005     .020461

     c.age#c.age    -.0003585   .0001212    -2.96   0.003     -.000596   -.0001209

        birth_yr     -.011967   .0021773    -5.50   0.000    -.0162344   -.0076997
           grade      .065723   .0030378    21.63   0.000      .059769     .071677
           _cons     1.255391   .1643564     7.64   0.000     .9332579    1.577523

    Instrumented:  tenure
    Instruments:   age c.age#c.age birth_yr grade union wks_work msp

The following code produces the results in tables 1 and 2, besides the cluster
bootstrap standard errors computed above (that are repeated for other
quantiles). Compared to the provided do-file, below I omit various timer and
matrix storage commands (and output). I run sivqr with both the default plug-
in bandwidth as well as bandwidth(0), which tries to find the smallest
computationally feasible bandwidth. The squared age variable is generated
manually because ivqreg and ivqreg2 do not support the factorial operator ##.

  *

        . generate agesq = age^2
    (24 missing values generated)

    . forvalues i = 1/3 {
      2.   local q = `i´/4
      3.   sivqr ln_wage age agesq birth_yr grade (tenure = union wks_work msp) ,
    >   quantile(`q´)
      4.   sivqr ln_wage age agesq birth_yr grade (tenure = union wks_work msp) ,
    >   quantile(`q´) bandwidth(0)
      5.   ivqreg ln_wage age agesq birth_yr grade (tenure = union wks_work msp) ,
    >   q(`q´) robust
      6. }

    . ivqreg2 ln_wage tenure age agesq birth_yr grade , quantile(0.25 0.5 0.75)
    >   instruments(union wks_work msp age agesq birth_yr grade)

Table 1: Estimates of tenure coefficient; clustered standard errors. Quantile | sivqr | sivqr, b(0) | ivqreg | ivqreg2
---|---|---|---|---
$0.25$ | $0.0865756$ | $0.0860257$ | $0.0861186$ | $0.0933732$
| (0.0031621) |  |  |
$0.50$ | $0.1076941$ | $0.1080343$ | $0.1032661$ | $0.1084423$
| (0.0046079) |  |  |
$0.75$ | $0.1565857$ | $0.1553029$ | $0.1584188$ | $0.1250514$
| (0.0103318) |  |  |
Note: 2SLS estimate (std. err.) is $0.1061$ ($0.0045$).

Table 1 shows the different tenure coefficient estimates at different
quantiles, with the 2SLS results for comparison. To help gauge the magnitude
of the differences, table 1 also shows cluster-robust standard errors from
2SLS and from using bootstrap with 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 sivqr standard error is also
similar to the 2SLS cluster-robust standard error.

Away from the median, sivqr (either bandwidth) and 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, ivqreg2 differs from the other IVQR estimators. This is especially
seen at the $0.75$ quantile: the 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 ivqreg2 is correctly specified and that this
difference is only due to larger estimation error of sivqr and ivqreg, but
given the consistency of both sivqr and 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 ivqreg2 location–scale model.

As a side note, these ivqreg2 results were using the regressor order that gave
the most economically plausible estimates. Keeping everything the same but
with tenure last in the regressor list, the estimates of the 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 2: Runtime (seconds) of IVQR commands on nlswork data. Quantile | bootstrap:sivqr | sivqr | sivqr, b(0) | ivqreg | ivqreg2
---|---|---|---|---|---
$0.25$ | $861$ | $16$ | $98$ | $130$ | $\smash{\uparrow}\phantom{64}$
$0.50$ | $500$ | $5$ | $68$ | $68$ | $1304$
$0.75$ | $1026$ | $4$ | $47$ | $153$ | $\smash{\downarrow}\phantom{64}$

Table 2 shows the runtime of the different commands at different quantiles.
The fastest is sivqr with the default plug-in bandwidth and analytic standard
errors, taking around $25$ seconds to compute at all three quantiles. Using
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 ivqreg time is yet slower than
bandwidth(0), also multiple minutes total. Naturally, running a cluster
bootstrap with $100$ replications takes even longer, even using the fast sivqr
with plug-in bandwidth, taking around $40$ minutes total. Even without
bootstrap, ivqreg2 takes over $20$ minutes, an order of magnitude more than
the $25$ seconds taken by sivqr.

###  4.2 More computation times

The provided do-file illustrates some additional important considerations
about computation time. Most of the sivqr computation time is from calling the
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 3 illustrates this by running sivqr over quantiles
$0.10,0.15,\ldots,0.90$, using either the plug-in bandwidth (always between
$0.05$ and $0.08$) or 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 bandwidth(0)
because it depends not only on the runtime per 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.

Table 3: Runtime (seconds) across quantile levels, nlswork data. Quantile | plug-in | bandwidth(0)
---|---|---
$0.10$ | $16$ | $99$
$0.15$ | $3$ | $51$
$0.20$ | $6$ | $49$
$0.25$ | $5$ | $88$
$0.30$ | $2$ | $82$
$0.35$ | $2$ | $54$
$0.40$ | $3$ | $46$
$0.45$ | $4$ | $25$
$0.50$ | $4$ | $61$
$0.55$ | $4$ | $29$
$0.60$ | $3$ | $27$
$0.65$ | $2$ | $29$
$0.70$ | $5$ | $38$
$0.75$ | $4$ | $44$
$0.80$ | $13$ | $36$
$0.85$ | $53$ | $44$
$0.90$ | $8$ | $43$

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 quantile(0.5), out of $100$ runs with different seeds,
runtimes for the built-in 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 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 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.

##  5 Methods and formulas

This section contains additional theoretical details. After introducing
notation in section 5.1, section 5.2 concerns the difficulty of estimation and
the solution of Kaplan and Sun (2017), while section 5.3 has some of the
identification arguments from Chernozhukov and Hansen (2005) that characterize
causal parameters as solutions to moment conditions (estimating equations).
Section 5.4 describes various bandwidths. Section 5.6 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 sivqr in
practice.

###  5.1 Notation

Notationally, let $y$ be the scalar outcome variable,
$\boldsymbol{\mathbf{{x}}}$ a column vector of regressors, and $u$ a scalar
unobserved variable. Let $\boldsymbol{\mathbf{{d}}}$ be a subset of
$\boldsymbol{\mathbf{{x}}}$ containing all the endogenous regressors. Let
$\boldsymbol{\mathbf{{z}}}$ be a column vector of all exogenous variables:
both exogenous regressors in $\boldsymbol{\mathbf{{x}}}$ as well as the
excluded instruments (instrumental variables). The scalar random variable $u$
is unobserved. Following Stata Journal convention, vectors are lowercase bold
(like $\boldsymbol{\mathbf{{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.

###  5.2 Estimation

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{{\beta}}}$ in
the structural model

$y=\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}+u$

If the instruments $\boldsymbol{\mathbf{{z}}}$ are independent of the
structural error $u$, then

$\boldsymbol{\mathbf{{0}}}=\operatorname{E}(u\mid\boldsymbol{\mathbf{{z}}})\implies\boldsymbol{\mathbf{{0}}}=\operatorname{E}(\boldsymbol{\mathbf{{z}}}u)=\operatorname{E}[\boldsymbol{\mathbf{{z}}}(y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}})]$
(2)

where the $\implies$ follows from the law of iterated expectations. If the
dimensions of $\boldsymbol{\mathbf{{z}}}$ and $\boldsymbol{\mathbf{{x}}}$ are
the same (“exact identification”), then this can be solved explicitly for
$\boldsymbol{\mathbf{{\beta}}}$:

$\operatorname{E}(\boldsymbol{\mathbf{{z}}}y)=\operatorname{E}(\boldsymbol{\mathbf{{z}}}\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}})\implies\operatorname{E}(\boldsymbol{\mathbf{{z}}}y)=\operatorname{E}(\boldsymbol{\mathbf{{z}}}\boldsymbol{\mathbf{{x}}}^{\prime})\boldsymbol{\mathbf{{\beta}}}\implies\boldsymbol{\mathbf{{\beta}}}=[\operatorname{E}(\boldsymbol{\mathbf{{z}}}\boldsymbol{\mathbf{{x}}}^{\prime})]^{-1}\operatorname{E}(\boldsymbol{\mathbf{{z}}}y)$

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{{z}}})$ in (2) with a conditional
$\tau$-quantile restriction for some $0<\tau<1$:

$0=\operatorname{Q}_{\tau}(u\mid\boldsymbol{\mathbf{{z}}})$ (3)

That is, conditional on any value of $\boldsymbol{\mathbf{{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{{z}}}$. Thus, (3) is equivalent to

$\tau=\Pr(u\leq
0\mid\boldsymbol{\mathbf{{z}}})=\Pr(y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}\leq
0\mid\boldsymbol{\mathbf{{z}}})$

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 (2),

$\displaystyle\tau=\operatorname{E}[\operatorname{\mathbf{1}}\\{y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}\leq
0\\}\mid\boldsymbol{\mathbf{{z}}}]$ $\displaystyle\implies
0=\operatorname{E}[\operatorname{\mathbf{1}}\\{y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}\leq
0\\}-\tau\mid\boldsymbol{\mathbf{{z}}}]$
$\displaystyle\implies\boldsymbol{\mathbf{{0}}}=\operatorname{E}[\boldsymbol{\mathbf{{z}}}(\operatorname{\mathbf{1}}\\{y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}\leq
0\\}-\tau)]$ (4)

Despite some similarities, computing an IVQR estimator based on (4) is much
more difficult than computing the mean IV estimator based on (2). With mean
IV, it was possible to solve for $\boldsymbol{\mathbf{{\beta}}}$ and replace
the mean (expectation) with the sample average. Here, the
$\boldsymbol{\mathbf{{\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
(4). 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\leq 0\\}$, $\tilde{I}(v)$ is a function of $v$
that smoothly decreases from $1$ to $0$ over $-1\leq v\leq 1$ instead of
decreasing discontinuously from $1$ to $0$ at $v=0$.222This notation is
simpler and equivalent to (but different than) that in Kaplan and Sun (2017).
Adding a bandwidth $h$ allows $\tilde{I}(v/h)$ to decrease more steeply over
$-h\leq v\leq h$. The estimator $\widehat{\boldsymbol{\mathbf{{\beta}}}}$
solves the “smoothed estimating equations” (or “smoothed moment conditions”)

$\boldsymbol{\mathbf{{0}}}=\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{\mathbf{{z}}}_{i}[\tilde{I}((y_{i}-\boldsymbol{\mathbf{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}})/h)-\tau]$
(5)

As a bonus, this smoothing also improves estimation precision, as shown both
theoretically and in simulations by Kaplan and Sun (2017, §§5,7). Weights
$w_{i}$ can be inserted into (5) as needed:

$\boldsymbol{\mathbf{{0}}}=\frac{1}{n}\sum_{i=1}^{n}w_{i}\boldsymbol{\mathbf{{z}}}_{i}[\tilde{I}((y_{i}-\boldsymbol{\mathbf{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}})/h)-\tau]$
(6)

Specifically, 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\leq-1$ and $\tilde{I}(v)=0$ for $v\geq 1$, but transitions linearly with
$\tilde{I}(v)=(1-v)/2$ for $-1<v<1$. Kaplan and Sun (2017, §2.2, p. 111) show
how this choice produces the Winsorized mean estimator of the type in Huber
(1964, ex. (iii), p. 79) in the special case with $\tau=0.5$ and an intercept-
only model ($\boldsymbol{\mathbf{{x}}}=\boldsymbol{\mathbf{{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{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}}$,
just as in the mean IV moment conditions; if $\tau\neq 0.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{{z}}}$ is used in (5). Specifically, it is
replaced by the linear projection of the regressor vector
$\boldsymbol{\mathbf{{x}}}$ onto the instruments $\boldsymbol{\mathbf{{z}}}$,
which has the same dimension as $\boldsymbol{\mathbf{{x}}}$ and thus
$\boldsymbol{\mathbf{{\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{{z}}}$ with the linear
projection of $\boldsymbol{\mathbf{{x}}}$ onto $\boldsymbol{\mathbf{{z}}}$.
Although theoretically more efficient estimators may exist, this estimator
remains consistent, reliable, and fast to compute regardless of the degree of
overidentification.

###  5.3 Identification

It remains to motivate the conditional quantile restriction in (3) from a
causal model, as originally done by Chernozhukov and Hansen (2005).333See also
Chapter 6 of Kaplan (2020) for an introductory discussion. Instead of assuming
every individual (or firm, or school, etc.) has the same coefficients
$\boldsymbol{\mathbf{{\beta}}}$, imagine each individual has their own
coefficient vector $\boldsymbol{\mathbf{{b}}}$. This is a “random coefficient”
model, meaning $\boldsymbol{\mathbf{{b}}}$ can differ by individual the same
way that $y$, $\boldsymbol{\mathbf{{x}}}$, and $\boldsymbol{\mathbf{{z}}}$ do:

$y=\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{b}}}$

An additional error term would be redundant; for example, if
$y=b_{0}+b_{1}x+v$, then the random intercept can simply absorb $v$ to become
$b_{0}+v$. To make the model more tractable empirically, imagine
$\boldsymbol{\mathbf{{b}}}$ can be written as a deterministic (but unknown)
vector-valued function $\boldsymbol{\mathbf{{\beta}}}(\cdot)$ applied to a
scalar unobserved $u$:

$y=\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(u)$ (7)

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{{\beta}}}(\cdot)$. Assume a “monotonicity” property that
given any value of $\boldsymbol{\mathbf{{x}}}$,
$\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(u)$ is an
increasing function of $u$, and like before assume the instruments are
independent of $u$. Then,

$\displaystyle\Pr(y\leq\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(\tau)\mid\boldsymbol{\mathbf{{z}}})$
$\displaystyle=\Pr(\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(u)\leq\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(\tau)\mid\boldsymbol{\mathbf{{z}}})$
by (7) $\displaystyle=\Pr(u\leq\tau\mid\boldsymbol{\mathbf{{z}}})$ by
monotonicity $\displaystyle=\Pr(u\leq\tau)$ $\displaystyle\textrm{by
}u\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\mathbf{{z}}}$
$\displaystyle=\tau$ $\displaystyle\textrm{by }u\sim\mathrm{Unif}(0,1)$

That is, for any $u=\tau$, $\boldsymbol{\mathbf{{\beta}}}(\tau)$ is identified
by the conditional quantile restriction from (3), which can be estimated by
the sivqr command using (5). Given $\tau$, the function
$q_{\tau}(\boldsymbol{\mathbf{{x}}})=\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(\tau)$
is also the $\tau$-quantile structural function introduced by Imbens and Newey
(2009, §3.1).

###  5.4 Plug-in bandwidth

This section only applies to the plug-in bandwidth. If instead the bandwidth
option is manually specified, then none of the following applies.

Kaplan and Sun (2017, prop. 2, p. 117) 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 Kaplan and Sun (2017).

A plug-in (data-dependent) version of the optimal bandwidth from Kaplan and
Sun (2017) is implemented in sivqr as follows. First, as in their proposition
2, the bandwidth is simplified by assuming $v\equiv
y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}(\tau)$ is
independent of the full vector of instruments. Second, the smoothed indicator
function is piecewise linear, as mentioned in section 5.2. In the notation of
Kaplan and Sun (2017),

$G(v)=\max\\{0,\min\\{1,(v+1)/2\\}\\}$ (8)

Third, the smallest amount of underlying smoothness in the data-generating
process is assumed, $r=2$ (see Kaplan and Sun, 2017, ass. 3).

With the chosen $G(\cdot)$ in (8) and $r=2$, simplifying the optimal bandwidth
$h^{*}$ from propositon 2 of Kaplan and Sun (2017) with
$v\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}\boldsymbol{\mathbf{{z}}}$,

$\displaystyle
1-\int_{-1}^{1}[G(v)]^{2}\mathop{}\\!\mathrm{d}{v}=1-\int_{-1}^{1}(v+1)^{2}/4\mathop{}\\!\mathrm{d}{v}=1/3$
$\displaystyle\left[\int_{-1}^{1}G^{\prime}(v)v^{r}\mathop{}\\!\mathrm{d}{v}\right]^{2}=\left[\int_{-1}^{1}v^{2}/2\mathop{}\\!\mathrm{d}{v}\right]^{2}=1/9$
$\displaystyle\begin{split}h^{*}&=\left(\frac{(r!)^{2}\left[1-\int_{-1}^{1}[G(v)]^{2}\mathop{}\\!\mathrm{d}{v}\right]f_{v}(0)}{2r\left[\int_{-1}^{1}G^{\prime}(v)v^{r}\mathop{}\\!\mathrm{d}{v}f_{v}^{(r-1)}(0)\right]^{2}}\frac{d}{n}\right)^{1/(2r-1)}\\\
&=n^{-1/3}\left(d\frac{4(1/3)f_{v}(0)}{4(1/9)\left[f_{v}^{\prime}(0)\right]^{2}}\right)^{1/3}=n^{-1/3}\left(3df_{v}(0)/[f_{v}^{\prime}(0)]^{2}\right)^{1/3}\end{split}$

where $f_{v}(\cdot)$ is the probability density function of $v$ and
$f_{v}^{\prime\prime}(\cdot)$ its second derivative, and $d$ is the dimension
(length) of $\boldsymbol{\mathbf{{\beta}}}(\tau)$.

The 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
(Silverman, 1986, §3.4.2) for estimating the density of $v$, as suggested by
Fernandes et al. (2021) 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 $\widehat{v}_{i}$ are then recomputed from the new
$\widehat{\boldsymbol{\mathbf{{\beta}}}}(\tau)$, leading to an updated plug-in
bandwidth, and the resulting IVQR estimate is reported by sivqr.

#### Nonparametrically estimated bandwidth

One approach is to replace $f_{v}(0)$ and $f_{v}^{\prime}(0)$ in $h^{*}$ by
nonparametric kernel estimates. Since $v$ is unobserved, an initial estimate
$\widehat{\boldsymbol{\mathbf{{\beta}}}}(\tau)$ must be used to construct
residuals

$\widehat{v}_{i}=y_{i}-\boldsymbol{\mathbf{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}}(\tau)$

For $f(0)$ (dropping the $v$ subscript for now), given the $\widehat{v}_{i}$,
the usual kernel density estimator is

$\frac{1}{ns}\sum_{i=1}^{n}K(-\widehat{v}_{i}/s)$ (9)

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 (Silverman,
1986, §3.4.2) 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 DasGupta (2008, Ch. 32, p. 536),
for example:

$s^{*}=n^{-1/5}\left(\frac{f(0)}{\\{f^{\prime\prime}(0)\\}^{2}}\right)^{1/5}\left(\frac{\int_{\mathbb{R}}\\{K(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}}{\\{\int_{\mathbb{R}}K(v)v^{2}\mathop{}\\!\mathrm{d}{v}\\}^{2}}\right)^{1/5}$

Following the Gaussian reference approach of Silverman (1986, §3.4.2), 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,

$\tau=\Phi((0-\mu)/\sigma)\implies-\mu/\sigma=\Phi^{-1}(\tau)$

Thus, the density at zero is

$f(0)=\phi((0-\mu)/\sigma)/\sigma=\phi(\Phi^{-1}(\tau))/\sigma$ (10)

Similarly, assuming normality allows us to express
$\widehat{f}^{\prime\prime}(0)$ in terms of $\widehat{\sigma}$. The second
derivative of the $\mathrm{N}(\mu,\sigma^{2})$ density evaluated at zero
is444D[PDF[NormalDistribution[\\[Mu], \\[Sigma]], x], {x,2}] at
http://www.wolframalpha.com for example; then plug in $x=0$ and
$-\mu/\sigma=\Phi^{-1}(\tau)$.

$\displaystyle f^{\prime\prime}(0)$
$\displaystyle=f(0)\\{(0-\mu)^{2}/\sigma^{4}-1/\sigma^{2}\\}=\overbrace{\phi(\Phi^{-1}(\tau))\sigma^{-1}}^{=f(0)\textrm{
by \eqref{eqn:f0}}}\sigma^{-2}[\\{\Phi^{-1}(\tau)\\}^{2}-1]$
$\displaystyle=\sigma^{-3}\phi(\Phi^{-1}(\tau))[\\{\Phi^{-1}(\tau)\\}^{2}-1]$

Thus,

$\displaystyle\frac{f(0)}{\\{f^{\prime\prime}(0)\\}^{2}}=\frac{\phi(\Phi^{-1}(\tau))\sigma^{-1}}{\sigma^{-6}\\{\phi(\Phi^{-1}(\tau))\\}^{2}[\\{\Phi^{-1}(\tau)\\}^{2}-1]^{2}}=\frac{\sigma^{5}}{\\{\phi(\Phi^{-1}(\tau))\\}[\\{\Phi^{-1}(\tau)\\}^{2}-1]^{2}}$

Also, for the Gaussian kernel,555Unfortunately, the value in table 32.1 in
DasGupta (2008, p. 537) is incorrectly stated as $1/\sqrt{2\pi}$.

$\displaystyle\int_{\mathbb{R}}\\{K(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}$
$\displaystyle=\frac{1}{2\pi}\int_{\mathbb{R}}\exp(-v^{2})\mathop{}\\!\mathrm{d}{v}=1/(2\sqrt{\pi}),$
$\displaystyle\int_{\mathbb{R}}K(v)v^{2}\mathop{}\\!\mathrm{d}{v}$
$\displaystyle=1$ (11)

Finally, plugging everything into the formula for $s^{*}$,

$\displaystyle s^{*}$
$\displaystyle=n^{-1/5}\left(\frac{f(0)}{\\{f^{\prime\prime}(0)\\}^{2}}\right)^{1/5}\left(\frac{\int_{\mathbb{R}}\\{K(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}}{\\{\int_{\mathbb{R}}K(v)v^{2}\mathop{}\\!\mathrm{d}{v}\\}^{2}}\right)^{1/5}$
$\displaystyle=n^{-1/5}\sigma\left[\phi(\Phi^{-1}(\tau))[\\{\Phi^{-1}(\tau)\\}^{2}-1]^{2}\right]^{-1/5}\left(\frac{1/(2\sqrt{\pi})}{1^{2}}\right)^{1/5}$
$\displaystyle=0.776n^{-1/5}\sigma\left[\phi(\Phi^{-1}(\tau))[\\{\Phi^{-1}(\tau)\\}^{2}-1]^{2}\right]^{-1/5}$

Following Silverman (1986, eq. (3.30)), to get a feasible bandwidth
$\widehat{s}$, the unknown $\sigma$ is replaced by either the sample standard
deviation of the $\widehat{v}_{i}$ or their sample interquartile range divided
by $1.349$ (the standard normal interquartile range), whichever is smaller.
This $\widehat{s}$ is then used to compute $\widehat{f}_{v}(0)$ using (9).

For the density derivative $f_{v}^{\prime}(0)$ (here simply $f^{\prime}(0)$),
the approach is the same: use a nonparametric kernel estimator with a
Silverman-type bandwidth. The usual estimator is

$\frac{1}{nb^{2}}\sum_{i=1}^{n}K^{\prime}(-\widehat{v}_{i}/b)$ (12)

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
Wand and Jones (1994),666Unfortunately, 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.

$n^{-1}b^{-3}f(0)\int_{\mathbb{R}}\\{K^{\prime}(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}+(1/4)b^{4}\\{f^{\prime\prime\prime}(0)\\}^{2}\left\\{\int_{\mathbb{R}}K(v)v^{2}\mathop{}\\!\mathrm{d}{v}\right\\}^{2}$

where $\int_{\mathbb{R}}K(v)v^{2}\mathop{}\\!\mathrm{d}{v}=1$ as in (11).
Solving the first-order condition (setting the derivative with respect to $b$
equal to zero), this leads to the optimal bandwidth

$b^{*}=n^{-1/7}\left(\frac{3f(0)\int_{\mathbb{R}}\\{K^{\prime}(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}}{\\{f^{\prime\prime\prime}(0)\\}^{2}}\right)^{1/7}$

Plugging in the standard normal density for $K(\cdot)$,

$\int_{\mathbb{R}}\\{K^{\prime}(v)\\}^{2}\mathop{}\\!\mathrm{d}{v}=1/(4\sqrt{\pi})$

Plugging in $f(0)$ from (10) and777D[PDF[NormalDistribution[\\[Mu],
\\[Sigma]], x], {x,3}] at http://www.wolframalpha.com for example; then plug
in $x=0$ and $-\mu/\sigma=\Phi^{-1}(\tau)$.

$f^{\prime\prime\prime}(0)=\sigma^{-3}f(0)\Phi^{-1}(\tau)\left[3-\\{\Phi^{-1}(\tau)\\}^{2}\right]$

yields

$\displaystyle b^{*}$
$\displaystyle=n^{-1/7}\left(\frac{3f(0)/(4\sqrt{\pi})}{\sigma^{-6}\\{f(0)\\}^{2}\\{\Phi^{-1}(\tau)\\}^{2}\left[3-\\{\Phi^{-1}(\tau)\\}^{2}\right]^{2}}\right)^{1/7}$
$\displaystyle=n^{-1/7}\left(\frac{3/(4\sqrt{\pi})}{\sigma^{-6}\sigma^{-1}\phi(\Phi^{-1}(\tau))\\{\Phi^{-1}(\tau)\\}^{2}\left[3-\\{\Phi^{-1}(\tau)\\}^{2}\right]^{2}}\right)^{1/7}$
$\displaystyle=n^{-1/7}\sigma\left(\frac{0.423}{\phi(\Phi^{-1}(\tau))\\{\Phi^{-1}(\tau)\\}^{2}\left[3-\\{\Phi^{-1}(\tau)\\}^{2}\right]^{2}}\right)^{1/7}$

Again, to get a feasible bandwidth $\widehat{b}$, $\sigma$ is replaced by the
sample standard deviation or interquartile range divided by $1.349$, whichever
is smaller. This $\widehat{b}$ is then used to compute
$\widehat{f}_{v}^{\prime}(0)$ using (12).

Finally, $\widehat{f}_{v}(0)$ and $\widehat{f}_{v}^{\prime}(0)$ are plugged
into the expression for $h^{*}$ to get the feasible plug-in bandwidth
$\widehat{h}$.

#### Gaussian reference bandwidth

Here, the Gaussian reference approach is used directly for $h^{*}$, to
simplify the density and density derivative of $v$.

From (10), 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 is888D[PDF[NormalDistribution[\\[Mu], \\[Sigma]], x], {x,1}] at
http://www.wolframalpha.com for example; then plug in $x=0$ and
$-\mu/\sigma=\Phi^{-1}(\tau)$.

$\displaystyle f^{\prime}(0)$ $\displaystyle=-\sigma^{-1}(-\mu/\sigma)f(0)$
$\displaystyle f(0)/\\{f^{\prime}(0)\\}^{2}$
$\displaystyle=\frac{1}{\sigma^{-2}\\{\Phi^{-1}(\tau)\\}^{2}\sigma^{-1}\phi(\Phi^{-1}(\tau))}=\frac{\sigma^{3}}{\\{\Phi^{-1}(\tau)\\}^{2}\phi(\Phi^{-1}(\tau))}$

Plugging this into the expression for $h^{*}$ yields

$\displaystyle
h^{*}=n^{-1/3}\sigma\left(\frac{3d}{\\{\Phi^{-1}(\tau)\\}^{2}\phi(\Phi^{-1}(\tau))}\right)^{1/3}$

and as usual $\sigma$ is replaced by the smaller of the sample standard
deviation of the $\widehat{v}_{i}$ or the sample interquartile range of the
$\widehat{v}_{i}$ divided by $1.349$.

###  5.5 Standard errors

By default or with reps(0), analytic standard errors are computed as follows;
otherwise, bootstrap is used as in section 5.6. The formulas here are special
cases of (3.11) of Chernozhukov and Hansen (2006); (18), (23), and Theorem 5
of de Castro et al. (2019); or Theorem 7 of Kaplan and Sun (2017). For
notational simplicity, I leave the dependence on $\tau$ implicit. Define an
additive error term

$\epsilon\equiv
y-\boldsymbol{\mathbf{{x}}}^{\prime}\boldsymbol{\mathbf{{\beta}}}.$ (13)

The first-order asymptotic distribution of the sivqr estimator is

$\sqrt{n}(\widehat{\boldsymbol{\mathbf{{\beta}}}}-\boldsymbol{\mathbf{{\beta}}})\stackrel{{\scriptstyle
d}}{{\to}}\mathrm{N}(\boldsymbol{\mathbf{{0}}},(\boldsymbol{\mathbf{{J}}}^{\prime}\boldsymbol{\mathbf{{S}}}^{-1}\boldsymbol{\mathbf{{J}}})^{-1}),\;\boldsymbol{\mathbf{{J}}}\equiv\operatorname{E}[f_{\epsilon|\boldsymbol{\mathbf{{z}}},\boldsymbol{\mathbf{{x}}}}(0\mid\boldsymbol{\mathbf{{z}}},\boldsymbol{\mathbf{{x}}})\boldsymbol{\mathbf{{z}}}\boldsymbol{\mathbf{{x}}}^{\prime}],\;\boldsymbol{\mathbf{{S}}}\equiv\tau(1-\tau)\operatorname{E}(\boldsymbol{\mathbf{{z}}}\boldsymbol{\mathbf{{z}}}^{\prime}),$
(14)

where
$f_{\epsilon|Z,X}(0\mid\boldsymbol{\mathbf{{z}}},\boldsymbol{\mathbf{{x}}})$
is the conditional PDF of $\epsilon$ given
$(\boldsymbol{\mathbf{{z}}},\boldsymbol{\mathbf{{x}}})$ evaluated at
$\epsilon=0$. The expression
$\boldsymbol{\mathbf{{J}}}^{-1}\boldsymbol{\mathbf{{S}}}(\boldsymbol{\mathbf{{J}}}^{\prime})^{-1}$
is mathematically equivalent but requires inverting the asymmetric matrix
$\boldsymbol{\mathbf{{J}}}$ by itself. The estimated asymptotic covariance
matrix is

$(\widehat{\boldsymbol{\mathbf{{J}}}}{}^{\prime}\widehat{\boldsymbol{\mathbf{{S}}}}{}^{-1}\widehat{\boldsymbol{\mathbf{{J}}}})^{-1},\quad\widehat{\boldsymbol{\mathbf{{S}}}}=\tau(1-\tau)\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{\mathbf{{z}}}_{i}\boldsymbol{\mathbf{{z}}}_{i}^{\prime},\quad\widehat{\boldsymbol{\mathbf{{J}}}}=\frac{1}{nh}\sum_{i=1}^{n}\phi(\frac{\boldsymbol{\mathbf{{y}}}_{i}-\boldsymbol{\mathbf{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}}}{h})\boldsymbol{\mathbf{{z}}}_{i}\boldsymbol{\mathbf{{x}}}_{i}^{\prime},$
(15)

using a conventional kernel density estimator for
$\widehat{\boldsymbol{\mathbf{{J}}}}$ with a Gaussian kernel $\phi(\cdot)$ and
bandwidth $h$ being Silverman’s rule of thumb (Silverman, 1986, §3.4.2),
specifically $1.06n^{-1/5}$ times the minimum of either the sample standard
deviation of the
$\boldsymbol{\mathbf{{y}}}_{i}-\boldsymbol{\mathbf{{x}}}_{i}^{\prime}\widehat{\boldsymbol{\mathbf{{\beta}}}}$
or their sample interquartile range divided by $1.349$. Alternatively, the
bandwidth in (A.8) of Liu (2019) would also work well. Note
$\widehat{\boldsymbol{\mathbf{{J}}}}$ is a Powell (1986) type estimator,
similar to (3.14) of Chernozhukov and Hansen (2006) but here with a Gaussian
instead of uniform kernel. The reported standard errors are then the square
roots of the diagonals of the matrix
$(\widehat{\boldsymbol{\mathbf{{J}}}}{}^{\prime}\widehat{\boldsymbol{\mathbf{{S}}}}{}^{-1}\widehat{\boldsymbol{\mathbf{{J}}}})^{-1}/n$.

###  5.6 Bayesian bootstrap

When reps(#) specifies multiple bootstrap replications, standard errors are
computed by Bayesian bootstrap (Rubin, 1981), assuming i.i.d. sampling. (With
non-i.i.d. sampling, bootstrap can be used; see [R] bootstrap.) This is a
particular type of frequentist exchangeable bootstrap (for example, van der
Vaart and Wellner, 1996, ex. 3.6.9, p. 354) that also has a nonparametric
Bayesian interpretation (for example, Chamberlain and Imbens, 2003, who
advocate for Bayesian bootstrap and show IV and QR applications). In each
replication, weights are set as $w_{i}=\xi_{i}/\overline{\xi}$ where the
$\xi_{i}$ are i.i.d. standard exponential random variables and
$\overline{\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 (6).
This is repeated many times; the reported standard error is the standard
deviation of the different $\widehat{\boldsymbol{\mathbf{{\beta}}}}$.

##  6 Conclusion

The sivqr command offers smoothed estimation of the instrumental variables
quantile regression model of Chernozhukov and Hansen (2005). The smoothing
improves computation and accuracy, helping sivqr to overcome several
limitations of ivqreg and ivqreg2, while complementing commands for
alternative quantile models with endogeneity (cqiv and ivqte). The sivqr
command can help Stata users reliably estimate heterogeneous effects across a
variety of settings and datasets.

##  7 Acknowledgements

There would be no sivqr Stata command without the work of Yixiao Sun on Kaplan
and Sun (2017). Thanks to Xin Liu for research assistance. Thanks to Joao
Santos Silva for suggesting the initial() option. Thanks to Di Liu for several
helpful suggestions and encouragement, including 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.

## References

  * de Castro et al. (2019) de Castro, L., A. F. Galvao, D. M. Kaplan, and X. Liu. 2019.  Smoothed GMM for quantile models.  _Journal of Econometrics_ 213(1): 121–144. URL https://doi.org/10.1016/j.jeconom.2019.04.008.
  * Chamberlain and Imbens (2003) Chamberlain, G., and G. W. Imbens. 2003.  Nonparametric Applications of Bayesian Inference.  _Journal of Business & Economic Statistics_ 21(1): 12–18. URL https://www.jstor.org/stable/1392346.
  * Cheng et al. (2019) Cheng, X., Z. Liao, and R. Shi. 2019.  On Uniform Asymptotic Risk of Averaging GMM Estimators.  _Quantitative Economics_ 10(3): 931–979. URL https://doi.org/10.3982/QE711.
  * Chernozhukov et al. (2012) Chernozhukov, V., I. Fernández-Val, S. Han, and A. Kowalski. 2012.  CQIV: Stata module to perform censored quantile instrumental variables regression.  Statistical Software Components S457478, Boston College Department of Economics.  Revised 25 Sep 2019. URL https://ideas.repec.org/c/boc/bocode/s457478.html.
  * Chernozhukov et al. (2019) . 2019.  Censored quantile instrumental-variable estimation with Stata.  _Stata Journal_ 19(4): 768–781. URL https://doi.org/10.1177/1536867X19893615.
  * Chernozhukov et al. (2015) Chernozhukov, V., I. Fernández-Val, and A. E. Kowalski. 2015.  Quantile regression with censoring and endogeneity.  _Journal of Econometrics_ 186(1): 201–221. URL https://doi.org/10.1016/j.jeconom.2014.06.017.
  * Chernozhukov and Hansen (2005) Chernozhukov, V., and C. Hansen. 2005.  An IV Model of Quantile Treatment Effects.  _Econometrica_ 73(1): 245–261. URL https://www.jstor.org/stable/3598944.
  * Chernozhukov and Hansen (2006) . 2006.  Instrumental quantile regression inference for structural and treatment effect models.  _Journal of Econometrics_ 132(2): 491–525. URL https://doi.org/10.1016/j.jeconom.2005.02.009.
  * Chernozhukov and Hansen (2008) . 2008.  Instrumental Variable Quantile Regression: A Robust Inference Approach.  _Journal of Econometrics_ 142(1): 379–398.
  * Chernozhukov et al. (2009) Chernozhukov, V., C. Hansen, and M. Jansson. 2009.  Finite sample inference for quantile regression models.  _Journal of Econometrics_ 152(2): 93–103. URL https://doi.org/10.1016/j.jeconom.2009.01.004.
  * Chernozhukov et al. (2017) Chernozhukov, V., C. Hansen, and K. Wüthrich. 2017.  Instrumental Variable Quantile Regression.  In _Handbook of Quantile Regression_ , ed. R. Koenker, V. Chernozhukov, X. He, and L. Peng, chap. 9, 119–144. CRC/Chapman-Hall. URL https://www.routledgehandbooks.com/doi/10.1201/9781315120256.
  * DasGupta (2008) DasGupta, A. 2008.  _Asymptotic Theory of Statistics and Probability_.  New York: Springer. URL https://doi.org/10.1007/978-0-387-75971-5.
  * Fernandes et al. (2021) Fernandes, M., E. Guerre, and E. Horta. 2021.  Smoothing Quantile Regressions.  _Journal of Business & Economic Statistics_ 39(1): 338–357. URL https://doi.org/10.1080/07350015.2019.1660177.
  * Frölich and Melly (2010) Frölich, M., and B. Melly. 2010.  Estimation of quantile treatment effects with Stata.  _Stata Journal_ 10(3): 423–457. URL https://www.stata-journal.com/article.html?article=st0203.
  * Frölich and Melly (2013) . 2013.  Unconditional quantile treatment effects under endogeneity.  _Journal of Business & Economic Statistics_ 31(3): 346–357.
  * Huber (1964) Huber, P. J. 1964.  Robust Estimation of a Location Parameter.  _Annals of Mathematical Statistics_ 35(1): 73–101. URL https://www.jstor.org/stable/2238020.
  * Imbens and Angrist (1994) Imbens, G. W., and J. D. Angrist. 1994.  Identification and Estimation of Local Average Treatment Effects.  _Econometrica_ 62(2): 467–475. URL https://www.jstor.org/stable/2951620.
  * Imbens and Newey (2009) Imbens, G. W., and W. K. Newey. 2009.  Identification and Estimation of Triangular Simultaneous Equations Models Without Additivity.  _Econometrica_ 77(5): 1481–1512. URL https://www.jstor.org/stable/25621368.
  * Kaplan (2020) Kaplan, D. M. 2020.  Distributional and Nonparametric Econometrics.  Textbook draft. URL https://kaplandm.github.io/teach.html.
  * Kaplan and Liu (2021) Kaplan, D. M., and X. Liu. 2021.  $k$-class IV quantile regression.  Working paper available at https://kaplandm.github.io/.
  * Kaplan and Sun (2017) Kaplan, D. M., and Y. Sun. 2017.  Smoothed Estimating Equations for Instrumental Variables Quantile Regression.  _Econometric Theory_ 33(1): 105–157. URL https://doi.org/10.1017/S0266466615000407.
  * Kwak (2010a) Kwak, D. W. 2010a.  ivqreg.ado (Stata code) version 1.0.0.
  * Kwak (2010b) . 2010b.  ivqreg.ado (Stata code) version 1.0.1.  Personal website.  Accessed Dec. 4, 2021. URL https://sites.google.com/site/dwkwak/dataset-and-code.
  * Lee (2007) Lee, S. 2007.  Endogeneity in quantile regression models: A control function approach.  _Journal of Econometrics_ 141(2): 1131–1158. URL https://doi.org/10.1016/j.jeconom.2007.01.014.
  * Liu (2019) Liu, X. 2019.  Averaging estimation for instrumental variables quantile regression.  Working paper available at https://ideas.repec.org/p/umc/wpaper/1907.html.
  * Machado and Santos Silva (2018) Machado, J. A. F., and J. M. C. Santos Silva. 2018.  IVQREG2: Stata module to provide structural quantile function estimation.  Statistical Software Components S458571, Boston College Department of Economics.  Revised 13 Oct 2021. URL https://ideas.repec.org/c/boc/bocode/s458571.html.
  * Machado and Santos Silva (2019) . 2019.  Quantiles via moments.  _Journal of Econometrics_ 213(1): 145–173. URL https://doi.org/10.1016/j.jeconom.2019.04.009.
  * Melly and Wüthrich (2017) Melly, B., and K. Wüthrich. 2017.  Local Quantile Treatment Effects.  In _Handbook of Quantile Regression_ , ed. R. Koenker, V. Chernozhukov, X. He, and L. Peng, chap. 10, 145–164. CRC/Chapman-Hall. URL https://www.routledgehandbooks.com/doi/10.1201/9781315120256.
  * Powell (1986) Powell, J. L. 1986.  Censored regression quantiles.  _Journal of Econometrics_ 32(1): 143–155.
  * Rubin (1981) Rubin, D. B. 1981.  The Bayesian Bootstrap.  _Annals of Statistics_ 9(1): 130–134. URL https://projecteuclid.org/euclid.aos/1176345338.
  * Silverman (1986) Silverman, B. W. 1986.  _Density Estimation for Statistics and Data Analysis_ , vol. 26 of _Monographs on Statistics and Applied Probability_.  London: Chapman & Hall/CRC.
  * StataCorp (2021) StataCorp. 2021.  Stata Statistical Software: Release 17.  College Station, TX: StataCorp LP.
  * van der Vaart and Wellner (1996) van der Vaart, A. W., and J. A. Wellner. 1996.  _Weak Convergence and Empirical Processes: With Applications to Statistics_.  Springer Series in Statistics, New York: Springer. URL https://doi.org/10.1007/978-1-4757-2545-2.
  * Wand and Jones (1994) Wand, M. P., and M. C. Jones. 1994.  _Kernel Smoothing_ , vol. 90 of _Monographs on Statistics and Applied Probability_.  Chapman & Hall/CRC. URL https://www.worldcat.org/title/kernel-smoothing/oclc/711630787/editions.
  * Wüthrich (2020) Wüthrich, K. 2020.  A Comparison of Two Quantile Models With Endogeneity.  _Journal of Business & Economic Statistics_ 38(2): 443–456. URL https://doi.org/10.1080/07350015.2018.1514307.

About the author

David M. Kaplan is an associate professor in the Department of Economics at
the University of Missouri. His primary research interest is econometric
methodology. In particular, he enjoys creating and advancing methods for
understanding changes and treatment effects on entire distributions (instead
of just averages), also a focus of his free open-source textbook and his first
Stata command, distcomp.