EconBase
← Back to paper

Mostly Harmless Machine Learning: Learning Optimal Instruments in Linear IV Models

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

58,289 characters

Mostly Harmless Machine Learning: Learning Optimal Instruments in Linear IV Models


\begin{titlepage}

\maketitle

\begin{abstract}
  We offer straightforward theoretical results that justify incorporating
  machine learning in the standard linear instrumental variable setting. The key
  idea is to use machine learning, combined with sample-splitting, to predict
  the treatment variable from the instrument and any exogenous covariates, and
  then use this predicted treatment and the covariates as technical instruments
  to recover the coefficients in the second-stage. This allows the researcher to
  extract non-linear co-variation between the treatment and instrument that may
  dramatically improve estimation precision and robustness by boosting
  instrument strength. Importantly, we constrain the machine-learned predictions
  to be linear in the exogenous covariates, thus avoiding spurious
  identification arising from non-linear relationships between the treatment and
  the covariates. We show that this approach delivers consistent and
  asymptotically normal estimates under weak conditions and that it may be
  adapted to be semiparametrically efficient \citep{chamberlain1992comment}. Our
  method preserves standard intuitions and interpretations of linear
  instrumental variable methods, including under weak identification, and
  provides a simple, user-friendly upgrade to the applied economics toolbox. We
  illustrate our method with an example in law and criminal justice, examining
  the causal effect of appellate court reversals on district court sentencing
  decisions.
\end{abstract}

\end{titlepage}

\section{Introduction}
\label{sec:intro}

Instrumental variable (IV) designs are a popular method in empirical economics.
Over 30\% of all NBER working papers and top journal publications considered by
\cite{currie2020technology} include some discussion of instrumental variables.
The vast majority of IV designs used in practice are linear IV estimated via
two-stage least squares (TSLS), a familiar technique in standard
introductions to econometrics and causal inference
\citep[e.g.][]{angrist2008mostly}. Standard TSLS, however, leaves on the table
some variation provided by the instruments that may improve precision of
estimates, as it only exploits variation that is linearly related to the
endogenous regressors. In the event that the instrument has a low linear
correlation with the endogenous variable, but nevertheless predicts the
endogenous variable well through a nonlinear transformation, we should expect
TSLS to perform poorly in terms of both estimation precision and inference
robustness. In particular, in some cases, TSLS would provide spuriously precise
but biased estimates \citep[due to weak instruments, see][]{andrews2018weak}.
Such nonlinear settings become increasingly plausible when exogenous variation
includes high dimensional data or alternative data, such as text, images, or
other complex attributes like weather. We show that off-the-shelf machine
learning techniques provide a general-purpose toolbox for leveraging such
complex variation, improving instrument strength and estimate quality.

Replacing the first stage linear regression with more flexible
specifications does not come without cost in terms of stronger identifying
assumptions. The validity of TSLS hinges only upon the restriction that the
instrument is linearly uncorrelated with unobserved disturbances in the response
variable. Relaxing the linearity requires that endogenous residuals are mean
zero conditional on the exogenous instruments, which is stronger. However, it is
rare that a researcher has a compelling reason to believe the weaker
non-correlation assumption, but rejects the slightly stronger mean-independence
assumption. Indeed, whenever researchers contemplate including higher
order polynomials of the instruments, they are implicitly accepting
stronger assumptions than TSLS allows.
In fact, by not exploiting the
nonlinearities, TSLS may accidentally
make a strong instrument weak, and deliver spuriously precise inference:
\cite{dieterle2016simple} and references therein find that several applied
microeconomics papers have conclusions that are sensitive to the specification
(linear vs. quadratic) of the first-stage.


A more serious identification concern with leveraging machine learning in the
first-stage comes from the parametric functional form in the second stage. When
there are exogenous covariates that are included in the parametric structural
specification, nonlinear transformations of these covariates could in principle
be valid instruments, and provide variation that precisely estimates the
parameter of interest.
For example, in the standard IV setup of $Y = D^\t\tau + X^\t\beta
+ U$ where $X$ is an exogenous covariate, imposing $\E[U \mid X] = 0$ would
formally result in $X^2, X^3$, etc. being valid ``excluded'' instruments.
However, given that the researcher's stated source of identification comes
from excluded instruments, such ``identifying variation'' provided by
covariates is more of an artifact of parametric specification than any serious
information from the data that relates to the researcher's scientific inquiry.

One principled response to the above concern is to make the second stage
structural specification likewise nonparametric, thereby including an infinite
dimensional parameter to estimate, making the empirical design a
\emph{nonparametric instrumental variable} (NPIV) design. Significant
theoretical and computational progress have been made in this regard
\citep[inter alia,][]{newey2003instrumental,ai2003efficient,ai2007estimation,horowitz2007nonparametric,severini2012efficiency,ai2012semiparametric,hartford2017deep,dikkala2020minimax,chen2012estimation,chen2015sieve,chernozhukov2018double,chernozhukov2016locally}.
However, regrettably, NPIV has received relatively little attention in applied
work in economics, potentially due to theoretical complications, difficulty in
interpretation and troubleshooting, and computational scalability. Moreover, in
some cases parametric restrictions on structural functions come from theoretical
considerations or techniques like log-linearization, where estimated parameters
have intuitive theoretical interpretation and policy relevance. In these cases
the researcher may have compelling reasons to stick with parametric
specifications.

In the spirit of being user-friendly to practitioners, this paper considers
estimation and inference in an instrumental variable model where the second
stage structural relationship is linear, while allowing for as much nonlinearity
in the instrumental variable as possible, without creating unintended and
spurious identifying variation from included covariates. Our results
provide intuition and justification for
using machine learning methods in instrumental variable designs. We show that
with sample-splitting, under weak consistency conditions, a simple estimator
that uses the predicted values of endogenous and included regressors as
technical instruments is consistent, asymptotically normal, and
semiparametrically efficient. The constructed instrumental variable also readily
provides weak instrument diagnostics and robust procedures. Moreover, standard
diagnostics like out-of-sample prediction quality are directly related to
quality of estimates. In the presence of exogenous covariates that are
parametrically included in the second-stage structural function, adapting
machine learning techniques requires caution to avoid spurious identification
from functional forms of the included covariates. To that end, we formulate and
analyze the problem as a sequential moment restriction, and develop estimators
that utilize machine learning for extracting nonlinear variation from and only
from instruments.

\paragraph{Related Literature.} The core techniques that allow for the construction of our
estimators follow from \cite{chamberlain1987asymptotic,chamberlain1992comment}.
The ideas in our proofs are also familiar in the double machine learning
\citep{chernozhukov2018double,belloni2012sparse} and semiparametrics literatures
\citep[e.g.][]{liu2020deep}; our arguments, however, follow from elementary
techniques that are accessible to graduate students and are self-contained. Our
proposed estimator is similar to the split-sample IV or jackknife IV estimators
in \cite{angrist1999jackknife}, but we do not restrict ourselves to linear
settings or linear smoothers. Using nonlinear or machine learning in the first
stage of IV settings is considered by \cite{xu2021weakiv} (for probit),
\cite{hansen2014instrumental} (for ridge),
\cite{belloni2012sparse,chernozhukov2015post} (for lasso), and
\cite{bai2010instrumental} (for boosting), among others; and our work can be
viewed as providing a simple, unified analysis for practitioners, much in the
spirit of \cite{chernozhukov2018double}. To the best of our knowledge, we are
the first to formally explore practical complications of making the first stage
nonlinear in a context with exogenous covariates. Finally, we view our work as
counterpoint to the recent work by \cite{angrist2019machine}, which is more
pessimistic about combining machine learning with instrumental variables---a
point we explore in detail in \cref{sub:disc}.


\section{Main theoretical results}

We consider the standard cross-sectional setup where the data $(R_i)_
{i=1}^N = (Y_i, D_i, X_i, W_i)_ {i=1}^N \iid P$ are sampled from some
infinite population. $Y_i$ is some outcome variable, $D_i$ is a set of
endogenous treatment variables, $X_i$ is a set of exogenous controls, and
$W_i$ is a set of instrumental variables. The researcher is willing to argue
that $W_i$ is exogenously or quasi-experimentally assigned. Moreover, the
researcher believes that $W_i$ provides a source of variation that
``identifies'' an effect $\tau$ of $D_i$ on $Y_i$. We denote the endogenous
variables and covariates as $T_i\equiv  [1, D_i^\t, X_i^\t]^\t$ and the
excluded instrument and covariates as the \emph{technical instruments} $Z_i
\equiv [1, W_i^\t, X_i^\t]^\t$.


A typical specification in empirical economics is the linear instrumental
variables specification:
\begin{align}
Y_i &= \alpha + D_i^\t \tau + X_i^\t \beta + U_i \quad \E[W_i U_i] = 0.
\label{eq:linIV}
\end{align}
We believe that often the researcher is willing to assume more than that $U_i$
is uncorrelated with $(X_i, W_i)$. Common introductions of instrumental
variables \citep{angrist2008mostly,angrist2001instrumental} stress that
instruments induce variation in $D_i$ and are otherwise unrelated to $U_i$, and
that a common source of instruments is natural experiments. We argue that these
narratives imply a stronger form of exogeneity than TSLS requires. After all, a
symmetric mean-zero random variable $S$ is uncorrelated with $S^2$, but one
would be hard pressed to say that $S^2$ is unrelated to $S$. Moreover, the condition $\E[W_i U_i] = 0$, strictly
speaking, does not automatically make polynomial expansions of $W_i$ valid
instruments, yet using higher order polynomials is common in empirical
research, indicating that the conditional restriction $\E[U_i \mid W_i] =
0$ more accurately captures the assumptions imposed in many empirical
projects.  With this in mind,
we will assume mean independence throughout the paper: $\E[U_i
\mid W_i] = 0$.
This stronger exogeneity assumption allows researcher to extract more
identifying variation from instruments, but doing so calls for more flexible
machinery for dealing with the first stage.\footnote{Moreover, under
conditions such that our proposed estimator is
efficient, we can test mean independence
assuming $\E[W_iU_i] = 0$, since TSLS and our proposed
estimator are two
estimators that generate a Hausman test.}


\subsection{No covariates}
\label{sub:simple}
Let us first consider the case in which we do not have exogenous
covariates $X_i$. Our mean-independence restrictions give rise to a conditional
moment restriction, $
\E[Y_i - T_i^\t \theta \mid W_i] = 0
$, where $\theta = (\alpha, \tau^\t, \beta^\t)^\t$. The conditional moment
restriction encodes an infinite set of unconditional
moment restrictions: \[
\text{For all square integrable $\tilde \upsilon$}: \E[\tilde \upsilon(W_i)(Y_i - T_i^\t\theta)] = 0.
\]
\cite{chamberlain1987asymptotic} finds that all relevant statistical information
in a conditional moment restriction is contained in a single unconditional
moment restriction involving an \emph{optimal instrument} $\Upsilon^\star$, and the
unconditional moment restriction with the optimal instrument delivers
semiparametrically efficient estimation and inference. In our case, $\Upsilon^\star(W_i) =
\frac{1}{\sigma^2(W_i)}[1, \mu(W_i)^\t]^\t$, where $
\mu(W_i) \equiv \E[D_i \mid W_i]$ and $\sigma^2(W_i) = \E[U_i^2 \mid
W_i]$. We estimate $\Upsilon^\star$
 with $\hat\Upsilon$ and form a plug-in estimator for $\theta$:
\begin{equation}
    \hat\theta_N = \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(W_i) T_i^\t}^{-1} \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(W_i)
    Y_i}.
\label{eq:plugin}
\end{equation}
This is numerically equivalent to estimating \eqref{eq:linIV} with two-stage
weighted least-squares with $\hat \mu(W_i)$ as an instrument and weighting
with $1/\sigma^2(W_i)$. In particular, if $U_i$ is homoskedastic, the
optimal
instrument is simply $[1, \mu(W_i)^\t]^\t$, and two-stage least-squares with an
estimate of $\mu(W_i)$ returns an estimate
$\hat\theta_N$.\footnote{This approach should not be confused with what many applied researchers think of when they think of two-stage least squares, namely directly regressing $Y_i$
on the estimated instrument $\hat\Upsilon(W_i)$ by OLS - i.e.     $\hat\theta_N = \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(W_i) \hat\Upsilon(W_i)^\t}^{-1} \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(W_i)
    Y_i}.$  This is what \cite{angrist2008mostly}
term the ``forbidden regression'', and it will not generally return
consistent estimates of $\theta$.} Under
heteroskedasticity, this instrument is no longer optimal (in the sense of
semiparametric efficiency) but remains valid. Therefore, we shall refer to the
instrument with the weighting $1/\sigma^2 (W_i)$ as the \emph{optimal instrument
under efficient weighting} and the instrument without $1/\sigma^2(W_i)$ as the
\emph{optimal instrument under identity weighting}.

Under identity-weighting, estimating $\Upsilon^\star$ amounts to learning $\mu
(W_i) \equiv \E[D_i \mid
W_i]$, which is well-suited
to machine learning techniques; this is only slightly complicated by the
estimation of $\sigma^2(W_i) \equiv \E[U_i^2 \mid W_i]$ under efficient
weighting. One might worry that the preliminary estimation of $\Upsilon^\star$ complicates
asymptotic analysis of $\hat \theta_N$. Under a simple sampling-splitting
scheme, however, we state a high-level condition for consistency, normality, and
efficiency of $\hat \theta_N$. Though it simplifies the proof and potentially
weakens regularity conditions, sample-splitting does reduce the amount of data
used to estimate the optimal instrument $\Upsilon^\star$, but such problems can be
effectively mitigated by $k$-fold sample-splitting: 20-fold sample-splitting,
for instance, limits the loss of data to $5\%$ at the cost of 20 computations
that can be effectively parallelized. Such concerns notwithstanding, we focus
our exposition to two-fold sample-splitting.

Specifically, assume $N=2n$ for simplicity and let $S_1, S_2
\subset [N]$ be the two subsamples with size $n$. Under
identity-weighting, for $j \in \{1,2\}$, form
$\hat\Upsilon^{(j)}$ by estimating $\mu(W_i)$ with data from the
other sample, $S_{-j}$. An estimator for $\mu$ may be a neural
network or a random forest trained via empirical risk minimization, or a
penalized linear regression such as elastic net.\footnote{With $k$-fold
sample-splitting, $S_{-j}$ is the union of all sample-split folds other than the
$j$-th one.} The estimated instrument $\hat\Upsilon_i$ is
then formed by evaluating $\hat\Upsilon^{(j)} (W_i)$ for all $ i
\in S_j$. We may then use \eqref{eq:plugin} to form an (identity-weighted)
estimator of $\theta$ by plugging in $\hat\Upsilon$. Under efficient weighting, on
each $S_{-j}$, we would use the identity-weighted estimator of $\theta$ as
an initial estimator to obtain an estimate of $U_i$, and similarly predict
$U_i^2$ with $W_i$ to form an estimate of $\sigma^2(W_i)$. The resulting
estimated optimal instrument under efficient weighting may then be plugged into
\eqref{eq:plugin} and form an efficient-weighted estimator. We term such
estimators the machine learning split-sample (MLSS) estimators. The pseudocode for
major procedures considered in this paper is collected in \cref{algo}.


\Cref{thm:simple} shows that the MLSS estimator is consistent and asymptotically
normal when the first-stage estimator $\hat\Upsilon^{(j)}$ converges to  a strong
instrument. Moreover, it is semiparametrically efficient when $\hat\Upsilon^{(j)}$ is
consistent for the optimal instrument $\Upsilon^\star(W_i) \equiv [1, \mu
(W_i)^\t]/\sigma^2(W_i)$ in $L^2(W)$ norm. The $L^2$ consistency
condition\footnote{In many-instrument settings under
\cite{bekker1994alternative}-type asymptotic sequences, there may be no
consistent estimator of the optimal instrument in the absence of sparsity
assumptions \citep{raskutti2011minimax}.} is not strong---in particular, it is
weaker than the $L^2$ consistency at $o(N^ {-1/4})$-rate conditions commonly
required in the  double machine learning and semiparametrics literatures
\citep{chernozhukov2018double},\footnote{We are not claiming that the MLSS
procedure has any advantage over the double machine learning literature, but
simply that the statistical problem here is sufficiently well-behaved such that
we enjoy weaker conditions than is typically required.} where such conditions
are considered mild.\footnote{The nuisance parameter $\E[D_i \mid W_i]$ in this
setting enjoys higher-order orthogonality property described in
\cite{mackey2018orthogonal}. In particular, it is infinite-order orthogonal,
thereby requiring no rate condition to work. Intuitively, estimation error in
$\Upsilon(\cdot)$ has no effect on the moment condition $\E[\Upsilon (Y_i -
\alpha - T_i^\t \theta)] = 0$ holding, and this feature of the problem makes the
estimation robust to estimation of $\Upsilon$.}

Formally, regularity conditions are stated in \cref{cond:simple}. The first
condition simply states that the nuisance estimation attains \emph{some} limit
as sample size tends to infinity, which is a similar requirement as \emph{sampling
stability} in \cite{lei2018distribution}. The second condition states that the
limit is a strong instrument. The third condition assumes bounded moments so as
to ensure a central limit theorem. The last condition, which is only required
for semiparametric efficiency, states that the nuisance estimation is consistent
for the optimal instrument in $L^2$ norm. For consistency of standard error
estimates, we assume more bounded moments in \cref{as:variance}.
\begin{as}
\label{cond:simple}
Recall that $Z_i = [1, W_i^\t]$, and so $\Upsilon(W_i)$ and $\Upsilon
(Z_i)$
denote the same object.
\begin{enumerate}
    \item ($\hat\Upsilon^{(j)}$ attains a limit $\Upsilon$ in $L^2$ distance) There
    exists some measurable function $\Upsilon(Z_i)$ such
    that \[
    \E{\norm{\hat\Upsilon^{(j)}(Z_i) - \Upsilon(Z_i)}^2 } \to 0 \quad \text{ for $j
    = 1,2$,}
    \]
    where the expectation integrates over both the randomness in $\hat\Upsilon^{(j)}$
    and in $Z_i$, but $\hat\Upsilon^{(j)}$ and $Z_i$ are assumed to be
    independent.

    \item (Strong identification) The matrix $G \equiv \E[\Upsilon(Z_i)
    T_i^\t]$ exists and is full
    rank.

    \item (Lyapunov condition) (i) For some $\epsilon > 0$, the following
    moments
    are finite:
    $\E|U_i|^{2+\epsilon} < \infty$,
    $\E\norm{U_i \Upsilon(Z_i)}^{2 + \epsilon} <
    \infty$, $\E[
  \norm{T_i}^2] < \infty$,
  and (ii) The variance-covariance matrix $\Omega \equiv \E[U_i^2 \Upsilon
  (Z_i)\Upsilon(Z_i)^\t]$
  exists, and (iii) the conditional variance is uniformly bounded: For some
  $M$, $\E[U_i^2 \mid Z_i] < M < \infty$ a.s.
  \item (Consistency to the optimal instrument) We may
  take the optimal instrument $\Upsilon^\star(Z_i)$ as the limit $\Upsilon(Z_i)$ in
  condition 1.
\end{enumerate}
\end{as}

\begin{as}[Variance estimation]
\label{as:variance}
Let $\Upsilon$ be the object defined in \cref{cond:simple}. Assume that the following fourth moments are bounded:
    \[\max\br{\E[\norm{T_i}^4], \E[U_i^4], \E\norm{\Upsilon(Z_i)}^4,
    \limsup_
    {N\to\infty} \E[\norm{\hat\Upsilon^{(j)}(Z_i)} ^4]} <
    \infty.\]
\end{as}


\begin{theorem}
\label{thm:simple}
  Let $\hat\theta^{\mathrm{MLSS}}_{N}$ be the MLSS estimator
  described above.  Under conditions 1--3 in \cref{cond:simple}, \[
 \sqrt{N} \pr{\hat\theta^{\mathrm{MLSS}}_{N} - \theta} \wto
  \Norm(0, V) \quad\quad V \equiv (G\Omega^{-1}G^\t)^{-1} = G^{-1} \Omega G^
  {-\t},
  \]
  where $G, \Omega$ are defined in \cref{cond:simple}. Moreover, if condition
  4 in \cref{cond:simple} holds, then the asymptotic variance $V$ attains
  the semiparametric efficiency bound. Moreover, if we additionally assume
  \cref{as:variance}, then the sample counterparts of $G, \Omega$ are
consistent for the two matrices.
\end{theorem}


\begin{proof}[Proof of \cref{thm:simple}]
We may compute that the scaled estimation error is\[
    \sqrt{N} \pr{\hat\theta^{\mathrm{MLSS}}_{N} - \theta} = \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(Z_i)
    T_i^\t}^{-1}
    \frac{1}{\sqrt{N}}\sum_{i=1}^N \hat\Upsilon(Z_i) U_i.
    \]
We verify in \cref{lemma:expansion} that, for $\Upsilon$ defined in
condition
1 of
\cref{cond:simple}, the following expansions hold:
\begin{equation}
    \label{eq:oracle}
    \hat G \equiv \frac{1}{N} \sum_{i=1}^N \hat\Upsilon(Z_i) T_i^\t = \frac{1}{N} \sum_{i=1}^N \Upsilon(Z_i)
    T_i^\t + o_p(1) \qquad \frac{1}{
\sqrt{N}}\sum_{i=1}^N \hat\Upsilon(Z_i) U_i = \frac{1}{\sqrt{N}}\sum_ {i=1}^N
\Upsilon(Z_i) U_i + o_p(1).
\end{equation}
Expansion \eqref{eq:oracle} implies that $\hat\theta^{\mathrm{MLSS}}_{N}$ is first-order equivalent
to the oracle estimator that plugs in $\Upsilon$: \[\hat \theta^\star_N
\equiv \pr{\frac{1}{N} \sum_{i=1}^N \Upsilon(Z_i)
    T_i^\t}^{-1}
    \frac{1}{N} \sum_{i=1}^N \Upsilon(Z_i) Y_i,\]
     whose consistency and asymptotic normality follows from usual
     arguments under condition 3 of \cref{cond:simple}.
Given \eqref{eq:oracle}, then we have a law of large numbers $\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(Z_i)
T_i^\t \pto G$ by condition
2 of \cref{cond:simple};
 and we obtain a central limit theorem $
\frac{1}{
\sqrt{N}}\sum_
{i=1}^N \Upsilon(Z_i) U_i \wto \Norm(0, \Omega)$ by condition 3. Lastly, by
Slutsky's theorem and the fact that $G$ is nonsingular, we obtain the
desired convergence $\smash{\sqrt{N}
\pr{\hat\theta^{\mathrm{MLSS}}_{N}
- \theta}  \wto \Norm(0,V)}$.

If, additionally, we assume the consistency condition 4, then
$\hat\theta^\star_N$ is exactly the efficient optimal instrument estimator
\citep{chamberlain1987asymptotic}, and hence $V$ attains the semiparametric
efficiency bound.
Finally, \eqref{eq:oracle} implies
that $\hat G \pto G$ via a weak law of large numbers, and
\cref{lemma:variance} implies $\hat \Omega \equiv \frac{1}{N} \sum_{i=1}^N
(Y_i - T_i^\t\hat\theta^{\mathrm{MLSS}}_{N})^2 \hat\Upsilon_i\hat\Upsilon_i^\t \pto \Omega$, and thus the variance
can be consistently estimated.
\end{proof}

\subsection{Exogenous covariates}

The presence of covariates $X_i$ complicates the analysis considerably. Under
the researcher's model, both $W_i$ and $X_i$ are considered exogenous, and
thus we may assume $\E[U_i \mid Z_i] = 0$ and use it as a conditional moment
restriction, under which the efficient instrument is $\var(U_i \mid Z_i)^ {-1}
\E[T_i \mid Z_i]$ and our analysis from the previous section continues to
apply \emph{mutatis mutandis}. However, if the researcher maintains a linear
specification $Y_i = T_i^\t \theta + U_i$, estimating $\theta$ based on the
conditional moment restriction $\E[U_i \mid Z_i] = 0$ may inadvertently
``identify'' $\theta$ through nonlinear behavior in $X_i$ rather than the
variation in $W_i$. Such a specification may allow the researcher to precisely
estimate $\theta$ even when the instrument $W_i$ is completely irrelevant, when,
say, higher-order polynomial terms in the scalar $X_i$, $X_i^2, X_i^3$,
are strongly correlated with $D_i$, perhaps due to misspecification of the
linear moment condition. There may well be compelling reasons why these
nonlinear terms in $X_i$ allow for identification of $\tau$ under an economic
or causal model; however, they are likely not the researcher's stated source of
identification, and allowing their influence to leak into the estimation
procedure undermines credibility of the statistical exercise.

One idea to resolve such a conundrum is to make the structural function
nonparametric as well, and convert the model to a nonparametric instrumental
variable regression
\citep{newey2003instrumental,ai2003efficient,ai2007estimation,ai2012semiparametric,chen2012estimation}
(See \cref{asec:npiv} for  discussion).\footnote{More recently,
\cite{chernozhukov2018double} derive Neyman-orthogonal moment conditions
assuming a partially linear second stage in Section 4.2 of their paper.} Another
idea, which we undertake in this paper, is to weaken the  moment condition and
rule nonlinearities in $X_i$ as inadmissible for inference.

To that end, we analyze the statistical restrictions implied by the
model and consider relaxations. The conditional moment restriction $\E[U_i
\mid Z_i] = 0$ is equivalent to the following orthgonality constraint
\begin{equation}
\label{eq:cmm}
  \text{For all (square integrable) $\Upsilon$},\, \E[\Upsilon(W_i, X_i)
(Y_i-T_i^\t
\theta)] = 0.
\end{equation}
Condition \eqref{eq:cmm} is too strong, since it allows nonlinear transforms of
$X_i$ to be
valid instruments. A natural idea is to restrict the class of allowable
instruments $\Upsilon (W_i, X_i)$ to those that are partially linear in
$X_i$, $\Upsilon(W_i, X_i) = h(W_i) + X_i^\t\ell$, thereby deliberately
discarding information from nonlinear transformations of $X_i$. Doing so
yields the following family of orthogonality constraints:
\begin{equation}
    \text{For all (square integrable) $\Upsilon$},\, \E[\Upsilon(W_i)
(Y_i-T_i^\t \theta)] = \E[X_i (Y_i - T_i^\t \theta)] = 0.
\label{eq:restriction}
\end{equation}
We may view \eqref{eq:restriction} as imposing an orthogonality condition
on the structural errors $U_i$ that is intermediate between that of TSLS
and that of \eqref{eq:cmm}. In particular, if we
define $\E_{\mathrm{(PL)}}[\cdot \mid
X_i, W_i]$ as a projection operator that projects onto partially linear
functions of $(X_i, W_i)$: \[
\E_{\mathrm{(PL)}}[U_i \mid X_i, W_i]\equiv \argmin_{\substack{h(X_i, W_i) \\ h
(X_i,
W_i) =
X_i^\t \ell + g(W_i)}} \E\pr{U_i - h(X_i, W_i)}^2,
\]
then requiring \eqref{eq:restriction} is equivalent to requiring
orthogonality under this partially linear projection operator:
\begin{equation}
\E_{\mathrm{(PL)}}[U_i \mid X_i, W_i] = 0.
\label{eq:pl_restriction}
\end{equation} In contrast, the $\cov(U_i, Z_i) = 0$ orthogonality
requirement of TSLS can be written as $\E_{\mathrm{(L)}} [\cdot \mid Z_i] = 0$, where
$\E_{\mathrm{(L)}}[\cdot \mid Z_i]$ is analogously defined as a projection operator onto
\emph{linear} functions of $Z_i$. We see that \eqref{eq:pl_restriction} is a
natural interpolation between the respective orthogonality structures on the
errors $U_i$ induced by the TSLS and the conditional moment restrictions.

The moment restriction corresponding to
\eqref{eq:restriction} is the following \emph{sequential moment restriction}
\begin{equation}
  \E[X_i (Y_i - T_i^\t \theta)] =  \E[Y_i - T_i^\t \theta \mid W_i] = 0.
\label{eq:seqmom}
\end{equation}
We see that \eqref{eq:seqmom} is a natural interpolation between
the usual
unconditional moment
condition, $\E[Z_i U_i] =
0$, and the conditional moment restriction that may be spurious $\E[U_i
\mid Z_i] =0$, by only allowing nonlinear information in $W_i$ to be
used for estimation and inference.

Having set up the estimation problem as (equivalently) characterized by
\eqref{eq:restriction},
\eqref{eq:pl_restriction}, or \eqref{eq:seqmom}, efficient estimation is
discussed by
\cite{chamberlain1992comment}. In particular, the optimal instrument under
identity weighting  takes the convenient
form
\begin{equation}
    \Upsilon^\star(Z_i) = \E_{\mathrm{(PL)}}\bk{T_i \mid X_i, W_i} = \colvecb{3}{1}{\E_{\mathrm{(PL)}}
\bk{D_i \mid X_i, W_i}}{X_i},
\label{eq:opt_covariate_homo}
\end{equation}
which is simply $(1, X_i)$, along with the best \emph{partially
linear prediction} of the endogenous
treatment $D_i$ from $W_i, X_i$. Observe that the only difference between
\eqref{eq:opt_covariate_homo} and
\cite{chamberlain1987asymptotic}'s optimal instrument under
homoskedasticity is modifying $\E$ into $\E_{\mathrm{(PL)}}$.
Implementing
\eqref{eq:opt_covariate_homo} is straightforward, as by
\cite{robinson1988root},\footnote{Moreover, it is easy to impose partial
linear structure on certain estimators, including series regression and
feedforward neural
networks, and in those cases we may minimize squared error directly
without Robinson's transformation.} partially linear regression is
reducible
to
fitting two nonparametric regressions $\E[D_i \mid W_i]$ and $\E[X_i
\mid W_i]$, and forming the following prediction function \[
\E_{\mathrm{(PL)}}\bk{D_i \mid X_i, W_i} = \E[D_i \mid W_i] + \E_{\mathrm{(L)}}[D_i \mid
\pr{X_i - \E[X_i \mid W_i]}].
\]

\begin{table}[tb]
    \caption{List of nonparametric nuisance parameters that require
    estimation. Note that nuisance parameters that require the unobserved
    error $U_i$ require additional preliminary consistent estimators of
    $\theta$.}
    \label{tab:nuisance}
    \vspace{1em}
    \centering

    \begin{tabular}{ccc}
    \toprule
        Covariates $X_i$ & Identity weighting? & Nonparametric
        nuisance
        parameters \\
        \midrule
         No  & Yes & $\E[D_i \mid W_i]$ \\
         No  & No & $\E[D_i \mid W_i]$, $\E[U_i^2 \mid W_i]$ \\
         Yes & Yes & $\E[D_i \mid W_i]$, $\E[X_i \mid W_i]$ \\
         Yes & No & $\E[D_i \mid W_i]$, $\E[X_i U_i^2 \mid W_i]$, $\E
         [U_i^2 \mid W_i]$
         \\
        \bottomrule
    \end{tabular}
\end{table}

Efficient estimation in the heteroskedastic case is more complex.
The optimal instrument is the vector \[
\Upsilon^\star(Z_i) = \frac{\E[T_i \mid W_i]}{\sigma^2(W_i)} + \E\bk{T_i \tilde
X_i^\t} \E
    \bk{U_i^2 \tilde X_i \tilde X_i^\t}^{-1}\tilde X_i  , \quad \tilde X_i \equiv X_i - \frac{\E[X_i U_i^2 \mid
    W_i]}
    {\sigma^2(W_i)}
\]
and the associated set of unconditional moment restrictions are \begin{equation}
    \E\bk{ U_i\cdot \Upsilon^\star(Z_i)
    } = 0.
    \label{eq:opt_cov_hetero}
\end{equation}
The intuition for \eqref{eq:opt_cov_hetero} is the following: the two
moment conditions $\E[U_i X_i] = \E[U_i \mid W_i] = 0$ provide
non-orthogonal information for $\theta$ that prevents us from applying the
optimal instrument on each moment condition. However, we may orthogonalize
one against the other.\footnote{Orthogonal here does not refer to Neyman orthogonality \citep{chernozhukov2018double}, but simply means that the two moments are uncorrelated.} In particular, the moment condition $\E[\tilde X_i
U_i] = 0$ is orthogonal to $\E[U_i \mid W_i]$ in the sense that $\E[\tilde X_i
U_i \cdot U_i \mid W_i] = 0$. Indeed, the term $\frac{\E[X_i U_i^2 \mid
    W_i]} {\sigma^2(W_i)}U_i \equiv \frac{\ip{X_iU_i, U_i}}{\ip{U_i,
U_i}} U_i $ is constructed to be the projection of $X_i U_i$
onto $U_i$ under the inner product $\ip{A,B} = \E[AB \mid W_i]$.

As before, complications in nonparametric estimation can be avoided by
sample splitting, where nuisance parameters are estimated on $K-1$ folds of
the data and the moment condition is evaluated on the remaining fold.
As a summary across our settings, we collect the nuisance parameters that
require a first-step estimation in
\cref{tab:nuisance}. The estimator \[
\hat\theta^{\mathrm{MLSS}}_{N} = \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(Z_i) T_i^\t }^{-1} \frac{1}{N} \sum_{i=1}^N \hat\Upsilon(Z_i) Y_i
\]
remains the same as \eqref{eq:plugin} and is subjected to the same
analysis in \cref{thm:simple}---under conditions 1--3 in
\cref{cond:simple}, $\hat\theta^{\mathrm{MLSS}}_{N}$
is consistent and asymptotically normal, and additionally, it is
semiparametrically efficient if the $L^2$-limit of $\hat\Upsilon$ coincides with
the optimal instrument $\Upsilon^\star$ in their respective settings.

Lastly, we make two remarks about the case with exogenous
covariates, assuming identity weighting for tractability. First, it is possible
for $\E[D_i \mid W_i] = 0$ and for \eqref{eq:opt_covariate_homo} to generate
precise estimates of the coefficient $\tau$. The reason is that it is possible
for the partially linear specification $D_i = h(W_i) +
X_i^\t \ell + V_i$ to generate nonzero $h (W_i)$ but zero conditional
expectation, in much the same way that some regression coefficients may be zero
without adjusting for $X_i$, but nonzero when adjusted for $X_i$. Whether or
not this makes $W_i$ a plausibly exogenous and strong instrument is likely to be
context specific. A robustness check may be generated by replacing $\E_{\mathrm{(PL)}}[D_i
\mid X_i, W_i]$ with $\E[D_i \mid W_i]$, which delivers consistent and
asymptotically normal estimates (assuming strong instrument) at the cost of
efficiency. Second, a Frisch-Waugh-Lovell- or double machine learning-like
procedure of first partialling out $X_i$ from $Y_i,
D_i$ and then treating \begin{equation}
    Y_i - \E_{\mathrm{(L)}}[Y_i \mid X_i] = \tau\pr{D_i - \E_{\mathrm{(L)}}[D_i \mid X_i]} + U_i \quad
\E[U_i \mid W_i] = 0
\label{eq:frisch-waugh}
\end{equation}
as a conditional moment restriction also delivers consistent and
asymptotically normal estimates. However, using the ``optimal'' instrument
for \eqref{eq:frisch-waugh}---the predicted residual $\E[D_i - \E_{\mathrm{(L)}}[D_i \mid
X_i] \mid W_i]$---\emph{does not} achieve semiparametric efficiency, since it
uses the information in the sequential moment restriction \eqref{eq:seqmom}
separately, without considering them jointly and orthogonalizing one against the
other, resulting in efficiency loss.


\begin{algorithm}[htb]
    \caption{Machine learning split-sample estimation and inference}
    \label{algo}
    \begin{algorithmic}
        \Require A subroutine PredictInstrument($S_{-j}, S_j$)
        that returns the estimated instrument $\{\hat\Upsilon(Z_i) : i\in
        S_j\}$, where $\hat\Upsilon$ is
        a function of $S_{-j}$.

         \vspace{1em}

        \Procedure{GenerateInstrument}{$K$, Data}
        \State Randomly split data into $S_1,\ldots,S_K$
        \For {$j$ in 1,\ldots, $K$}
            \State $\hat\Upsilon^{(j)} \gets \mathrm{PredictInstrument}(S_{-j},
            S_j)$
        \EndFor
        \State Combine $\hat\Upsilon^{(j)}$ into $\hat\Upsilon$
        \State Return $S_1,\ldots,S_K, \hat\Upsilon$
        \EndProcedure

         \vspace{1em}

        \Procedure{MLSSEstimate}{$K$, Data}
        \State $S_1,\ldots,S_K, \hat\Upsilon \gets
        \mathrm{GenerateInstrument}(K, \mathrm{Data})$

        \State For the full parameter vector, return \[\hat\theta = \pr{\frac{1}{N} \sum_{i=1}^N
        \hat\Upsilon(W_i) T_i^\t}^{-1} \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon(W_i)
    Y_i}\] and variance estimate \[
    \hat V = \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon_i T_i^\t}\cdot \pr{\frac{1}{N} \sum_{i=1}^N (Y_i - T_i^\t
    \hat\theta)^2 \hat\Upsilon_i \hat\Upsilon_i^\t}
    \cdot
    \pr{\frac{1}{N} \sum_{i=1}^N \hat\Upsilon_i T_i^\t}^{-\t}
    \]



        \State For the subvector $\hat\tau$, residualize $\hat\Upsilon_i, Y_i,
        D_i$ against $X_i$ to obtain $\tilde \upsilon_i, \tilde Y_i,
        \tilde D_i$, and compute $\hat\tau = \pr{\frac{1}{N} \sum_{i=1}^N \tilde \upsilon_i
        \tilde D_i^\t}^{-1} \pr{\frac{1}{N} \sum_{i=1}^N \tilde \upsilon_i\tilde Y_i}.$

        \Comment{\textbf{Assuming identity weighting}}

        \EndProcedure

        \vspace{1em}

        \Procedure{WeakIVInference}{$K$, Data, $\alpha$}

            \Comment{\textbf{Assuming identity weighting}}

            \Comment{Assuming a routine AndersonRubin($\alpha,
            \mathrm{Data}$) that
            returns the $1-\alpha$ Anderson--Rubin CI}

            \State $S_1,\ldots,S_K, \hat\Upsilon \gets
        \mathrm{GenerateInstrument}(K, \mathrm{Data})$

            \For {$j$ in 1,\ldots, $K$}
            \State On $S_j$, residualize $\hat\Upsilon^{(j)}(Z_i), Y_i,
        D_i$ against $X_i$ to obtain $\tilde \upsilon_i, \tilde Y_i,
        \tilde D_i$
            \State $\mathrm{CI}_j \gets \mathrm{AndersonRubin}\pr{\alpha /
            K,
            \br{\tilde \upsilon_i, \tilde Y_i,
        \tilde D_i}_{i\in S_j}}$
            \EndFor
            \State Return $\mathrm{CI} = \bigcap_j \mathrm{CI}_j$
        \EndProcedure
    \end{algorithmic}
\end{algorithm}


\section{Discussion}
\label{sub:disc}

\textbf{``Forbidden regression.''} Nonlinearities in the first stage are often
discouraged due to a ``forbidden regression,'' where the researcher regresses
$Y$ on \smash{$\hat D$} estimated via nonlinear methods, motivated by a
heuristic explanation for TSLS. As \cite{angrist2001instrumental} point out,
this regression is inconsistent, and consistent estimation follows from using
$\hat D$ as an instrument for $D$, as we do, rather than replacing $D$ with
$\hat D$---in the case where the first-stage is linear, the two estimates are
numerically equivalent, but not in general.

\paragraph{Interpretation under heterogeneous treatment effects.}
Assume $D_i$ is binary and suppose $Y_i = D_i Y_i(1) + (1-D_i)Y_i(0)$. Suppose
the treatment follows a Roy model, $D_i = \one
\pr{\mu(W_i) \ge V_i}$, where $V_i \sim \Unif[0,1]$. In this
setting, the conditional moment restriction \eqref{eq:linIV} is
misspecified, since it assumes constant treatment effects, and different
choices of the instrument would generate estimators that converge to
different population quantities. Nevertheless, the results of
\cite{heckman2005structural} (Section 4) show that different choices of the
instrument generate estimators that estimate different weightings of
\emph{marginal treatment effects} (MTEs); moreover, optimal instruments,
whether under identity weighting or optimal weighting, correspond to convex
averages of MTEs, whereas no such guarantees are available for linear IV
estimators with $W_i$ as the instrument, without assuming that $\E[D_i \mid
W_i]$ is linear. The weights on the MTEs are explicitly stated in
\cref{asec:hte}.

\paragraph{Weak IV detection and robust inference.}

A major practical motivation for our work, following \cite{bai2010instrumental},
is to use machine learning to rescue otherwise weak instruments due to a lack of
linear correlation; nonetheless, the instrument may be irredeemably weak, and
providing weak-instrument robust inference is important in practice. Relatedly,
\cite{xu2021weakiv} and \cite{antoine2019identification} also consider weak IV
inference with nonlinear first-stages; the benefits of split-sampling in the
presence of many or weak instruments are recently exploited by
\cite{mikusheva2020inference} and date to \cite{dufour2003identification},
\cite{angrist1999jackknife}, \cite{staiger1994instrumental}, and references
therein; \cite{kaji2019theory} proposes a general theory of weak identification
in semiparametric settings.

On weak IV detection, our procedure produces estimated optimal instruments,
which result in just-identified moment conditions. As a result, in models with
\emph{a single endogenous treatment variable}, the \cite{stock2005testing}
$F$-statistic rule-of-thumb has its exact interpretation\footnote{Namely, the
worst-case bias of TSLS exceeds 10\% of the worst-case bias of OLS
\citep{andrews2018weak}.} regardless of homo- or heteroskedasticity
\citep{andrews2018weak}, and the first stage $F$-statistic may be used as a tool
for detecting weak instruments.

Pre-testing for weak instruments distorts downstream inferences. Alternatively,
weak IV \emph{robust} inferences, which are inferences of $\tau$ that are valid
regardless of instrument strength, are often preferred. The procedure we
propose, under identity-weighting, is readily compatible with simple robust
procedures. In particular, on each subsample $S_j$, we may perform the
Anderson--Rubin test \citep{anderson1949estimation} and combining the results
across subsamples via Bonferroni correction. For a confidence interval at the
95\% nominal level with two-fold sample-splitting, this amounts to intersecting
two 97.5\%-nominal AR confidence intervals on each subsample, and these
confidence intervals may be computed using off-the-shelf software
implementations of AR confidence intervals.

More formally, consider the null hypothesis $H_0: \tau = \tau_0$. Consider
a Frisch--Waugh--Lovell procedure that partials out the covariates $X_i$.
Let the residuals be $U_i
(\tau_0)
\equiv Y_i
-D_i^\t \tau_0$, and let $\tilde U_i(\tau_0) \equiv U_i(\tau_0) -
 \tilde\delta^\t \bar X_i$ be the residual $U_i
(\tau_0)$ after partialling out $\bar X_i \equiv [1, X_i^\t]^\t$. Suppose the
estimated instrument takes the form $\hat\Upsilon(Z_i) = [1,
\hat\upsilon(Z_i)^\t,
X_i^\t]^\t$, where $\dim\hat\upsilon(Z_i) = \dim D_i$; this requirement is satisfied under
identity weighting. Similarly, partial out $\bar X_i$ from the estimated
instrument to obtain $\tilde \upsilon^{(j)}
\equiv \hat\upsilon^{(j)} - \tilde \lambda^\t \bar X_i$.
Finally, consider the covariance $V_{n,j}
(\tau_0)$ between
the residual and the instrument, after partialling out the covariates
$\bar X_i$, and let $\Omega_{n,j}$ be an estimate of $V_{n,j}$'s variance matrix:
i.e.,
\[
V_{n,j}(\tau_0) \equiv \frac{1}{\sqrt{n}} \sum_{i\inS_j}
\hat\upsilon^{(j)}
(Z_i)\tilde U_i(\tau_0) \text{ and } \Omega_{n,j}(\tau_0) \equiv \frac{1}
{n}\sum_{i\in S_j} \tilde U_i (\tau_0)^{2}
\tilde \upsilon^{(j)} (Z_i)\tilde \upsilon^{(j)} (Z_i)^\t.
\]
The Anderson--Rubin
statistic on the $j$-th subsample is then defined as the normalized
magnitude of $V_{n,j}$: $\operatorname{AR}_j(\tau_0) \equiv V_{n,j}^\t
\Omega_{n,j}^{-1} V_{n,j}.$ Under $H_0$, by virtue of the exclusion restriction, we
should expect $V_{n,j}$ be mean-zero Gaussian, and thus $\operatorname{AR}_j$ should be
$\chi^2$.
Indeed, \Cref{thm:weakiv} shows that on each subsample,
under mild bounded moment conditions that ensure convergence
(\cref{as:weakiv}), $\operatorname{AR}_j(\tau_0)$ attains a limiting $\chi^2$
distribution. Under weak IV asymptotics, it is not necessarily the case that
the AR statistics are asymptotically uncorrelated across subsamples, and so
we resort to the Bonferroni procedure in outputting a single confidence
interval.
\begin{as}[Bounded moments for the AR statistic]
    \label{as:weakiv}
    Without loss of generality and normalizing if necessary, assume the
    estimated instruments are normalized: $\sum_
    {i\inS_j}\hat\upsilon^{(j)}_k
    (Z_i)^2 = 1$ for all $k = 1, \ldots,\dim D_i$.
    Let $\lambda_n \equiv \E[\bar X_i
\bar X_i]^{-1} \E
[\hat\upsilon^{(j)}(Z_i) \bar X_i^\t \mid \hat\upsilon^{(j)}]$ be the projection coefficient of
$\hat\upsilon^{(j)}(Z_i)$ onto $\bar X_i$.
Assume that with probability 1, the sequence $\hat\upsilon^{(j)} =
    \hat\upsilon^{(j)}_n$ satisfies the Lyapunov conditions
    \begin{enumerate}[label=(\roman*)]
        \item $\E[U_i^{4} \norm{\hat\upsilon^{(j)}(Z_i) -
    \lambda_n \bar X_i}^
    {4} \mid \hat\upsilon^{(j)}] < C_1 < \infty$ for some $C_1 > 0$
    \item $\E[U_i^2 (\hat\upsilon^{(j)}(Z_i) - \lambda_n \bar X_i) (\hat\upsilon^{(j)}
        (Z_i)- \lambda_n \bar X_i)^\t \mid \hat\upsilon^{(j)}] \pto
    \Omega$.
    \end{enumerate}
    Moreover, assume (iii) $\max\br{\E[\norm{\hat\upsilon^{(j)}}^4], \E[U_i^4], \E
    [\norm{X_i}^4]} < C_2 < \infty$ and that (iv) $\E
    [\bar X_i\bar X_i^\t]$
    is invertible.
\end{as}
\begin{theorem}
\label{thm:weakiv}
    Under \cref{as:weakiv}, $\operatorname{AR}_j(\tau_0) \wto \chi^2_{\dim
    D_i}$.
\end{theorem}
\begin{proof}
    We relegate the proof, which amounts to checking convergences $V_{n,j}
    \wto \Norm(0,\Omega)$ and $\Omega_{n,j} \pto \Omega$ under \cref{as:weakiv}, to
    the appendix.
\end{proof}

















\paragraph{Connection between first-stage fitting and estimate quality.}
By using machine learning in the first stage, one may be able to improve the quality of the first-stage fit, as measured by out-of-sample $R^2$.
We now offer an argument as to why improving that fit may improve the mean
squared error of the estimator.

Consider a setting with no covariates\footnote{In the presence of covariates
under identity weighting, we may, without loss of generality, partial out the
covariates after estimating the optimal instrument.} and i.i.d. $U_i$ and an
estimated instrument $\hat\upsilon(Z_i)$,\footnote{i.e. $\hat\Upsilon(Z_i) = [1,
\hat\upsilon(Z_i)]^\t.$} meant to approximate $\E[D_i \mid Z_i]$. In linear IV, $\hat\upsilon(Z_i)$ is the linear projection of
$D_i$ onto $Z_i$. Define the \emph{extra-sample error}\footnote{See, for
instance, \cite{friedman2001elements} Section 7.4.} of an estimator $\hat\tau$
based on $\hat\upsilon (Z_i)$ to be the random quantity \begin{align*}
\mathrm{Err}(\hat\upsilon) &\equiv { n \cdot \pr{\frac{\cov_n (\hat\upsilon(Z_i) , Y_i)}{ \cov_n
   (\hat\upsilon(Z_i),
   D_i)} - \tau}^2 }
\end{align*} where $(Y_i,D_i,Z_i)_ {i=1}^n$ is a new and independent sample
unrelated to the estimate $\hat\upsilon$, and we hold $\hat\upsilon$ fixed. The subscript
$n$ denotes sample quantities such as sample variances and covariances. The
quantity $\mathrm{Err} (\hat\upsilon)$ is an optimistic measure of the quality of using
$\hat\upsilon$ as the instrument, as construction of $\hat\upsilon$ without
sample-splitting typically introduces a bias term since
$\cov(\hat\upsilon, U_i)$ cannot be assumed to be small. The following
calculation shows
that  $\mathrm{Err}(\hat\upsilon)$ scales with the inverse out-of-sample $R^2$ of $\hat\upsilon$
as a predictor of $D_i$:
\begin{align*}
\mathrm{Err}(\hat\upsilon) &= { n \frac{\cov_n (\hat\upsilon(Z_i) , U_i)^2}{ \cov_n (\hat\upsilon(Z_i),
   D_i)^2}
} \\ &= {
    \frac{1}{R^2_{n}\pr{\hat\upsilon(Z_i), D_i} } \cdot \var_n (D_i)^{-1} \pr{
    \frac{1}
    {\sqrt{n}}\sum_{i=1}^n \frac{\hat\upsilon(Z_i)}{\sqrt{\var_n (\hat\upsilon
    (Z_i))}} (U_i - \bar U)}^2
}\\
&\wto \frac{1}{R^2\pr{\hat\upsilon(Z_i), D_i}}\frac{\var(U_i)}{\var(D_i)}
\cdot \chi^2_1 \quad \text{as $n \to \infty$},
\end{align*}
where $R^2_n$ is the out-of-sample $R^2$ of predicting $D_i$ with
$\hat\upsilon(Z_i)$, and $R^2(\cdot,\cdot)$ is its limit in
probability.\footnote{The $F$-statistic is a monotone transformation of
the $R^2$, which also serves as an indication of estimation quality.}

The out-of-sample $R^2$, which can be readily computed from a split-sample
procedure, therefore offers a useful indicator for quality of estimation. In
particular, if one is comfortable with the strengthened identification
assumptions, there is little reason not to use the model that achieves the best
out-of-sample prediction performance on the split-sample. In some settings, this
best-performing model will be linear regression, but in many settings it may
not be, and attempting more complex tools may deliver considerable
benefits.

Moreover, much of the discussion on using machine learning for
instrumental variables analysis focuses on \emph{selecting} (relevant or
valid) instruments
\citep{belloni2012sparse,angrist2019machine} assuming some level of sparsity,
motivated by statistical difficulties encountered when the number of instruments
is high.  In light of the heuristic above, a more precise framing is perhaps
\emph{combining} instruments to form a better prediction of the endogenous
regressor, as noted by \cite{hansen2014instrumental}.





\paragraph{(When) is machine learning useful?} We conclude this section by
discussing our work relative to \cite{angrist2019machine}, who note that using
lasso and random forest methods in the first stage does not seem to provide
large performance benefits in practice, on a simulation design based on the data
of \cite{angrist1991does}. We note that, per our discussion above in the
connection between first-stage fitting and estimate quality, a good heuristic
summary for the estimation precision is the $R^2$ between the fitted instrument
and the true optimal instrument---$\E[D_i\mid W_i]$ in the homoskedastic case.
It is quite possible that in some settings, the conditional expectation $\E [D_i
\mid W_i]$ is estimated well with linear regression, and lasso or random forest
do not provide large benefits in terms of out-of-sample prediction quality.
Since \cite{angrist1991does}'s instruments are quarter-of-birth interactions and
are hence binary, it is in fact likely that predicting $D$ with  linear
regression performs well relative to nonlinear or complex
methods\footnote{Indeed, in some of our  experiments calibrated to the
\cite{angrist1991does} data, using a simple gradient boosting method
(\texttt{lightgbm}) does not outperform linear regression in terms of
out-of-sample $R^2$ (0.039\% vs. 0.05\%).} in the setting. Whether or not
machine learning methods
work well relative to linear methods is something that the researcher may verify
in practice, via evaluating performance on a hold-out set, which is standard
machine learning practice but is not yet widely adopted in empirical economics.
Indeed, we observe that in both real (\cref{sec:emp}) and Monte Carlo
(\cref{asec:monte_carlo}) settings where the out-of-sample prediction quality of
more complex machine learning methods out-perform linear regression, MLSS
estimators perform better than TSLS.


\section{Empirical Application}
\label{sec:emp}

We consider an empirical application in the criminal justice setting of
\cite{deepivinlaw}, where we consider the causal effect of appellate court
decisions at the U.S. circuit court level on lengths of criminal sentences at
the U.S. district courts under the jurisdiction of the circuit court.
\cite{deepivinlaw} exploit the fact that appellate judges are randomly assigned,
and use the characteristics of appellate judges---including age, party
affiliation, education, and career backgrounds---as instrumental variables. In
criminal justice cases, plaintiffs rarely appeal, as it would involve trying the
defendant twice for the same offense---generally not permitted in the United
States; therefore, an appellate court reversal would
typically be in favor of defendants, and we may posit a causal channel in which
such reversals affect sentencing; for instance, the district court may be more
lenient as a result of a reversal, as would be naturally
predicted if the reversal sets up a precedent in favor of the defendant.

To connect the empirical setting with our notation, the outcome variable $Y$ is
the change in sentencing length before and after an appellate decision, measured
in months, where positive values of $Y$ indicates that sentences become longer
after the appellate court decision. The endogenous treatment variable $D$ is
whether or not an appellate decision reverses a district court ruling. The
instruments $W$ are the characteristics of the randomly assigned circuit judge
presiding over the appellate case in question, and covariates $X$ contain
textual features from the circuit case, represented by Doc2Vec embeddings
\citep{le2014distributed}.


We compute two estimators of the optimal instrument under identity weighting
based on flexible methods that are often characterized as machine learning
(random forest and LightGBM, for light gradient boosting machine), as well as a
variety of linear or polynomial regression estimators, with or without
sample-splitting. We present our results in \cref{fig:judge} (see the notes to the figure for precise definitions of each estimator). For all of the split-sample estimators we have three sets of point estimates and confidence intervals, corresponding to the different splits.  We also have specifications that exclude (panel (a)) and include covariates (panel (b)); except where discussed below the results are consistent across both.

The machine learning estimators perform similarly across splits, reporting
mildly negative point estimates
of between $1$ and $2$ months reduction in sentencing, with confidence intervals that are reasonably tight but include a zero effect. Moreover, the Wald and Anderson--Rubin
confidence intervals are similar, suggesting that the instrument constructed by the ML methods is sufficiently strong
to result in inferences that are not distorted.

A natural benchmark to compare these results to is TSLS, with the instruments
entering either linearly or quadratically (but not interacted). Linear TSLS
estimates a slightly positive effect without covariates and a slightly negative
one with them. Though it has a tight Wald confidence interval, the AR interval
is quite large ($[-3.3, 23.8]$). This indicates a weak instruments
problem, and indeed the first-stage $F$-statistic is only 1.5.
Quadratic TSLS returns
a point estimate that is close to the MLSS estimates, with a very tight Wald
confidence interval; however, the corresponding Anderson--Rubin interval is
empty. Anderson--Rubin tests
test model misspecification jointly with point nulls of the structural
coefficient, and may report an empty interval if the model is
misspecified.\footnote{In our case, the endogenous treatment is binary,
and so the only source of model misspecification is heterogeneous treatment
effects. In that case, TSLS continues to estimate population objects
that are (possibly nonconvex) averages of marginal treatment effects, and
arguably researchers would nonetheless like non-empty confidence sets. One
benefit of split-sample approaches is that the power of the
Anderson--Rubin test is wholly directed to testing the structural parameter
rather than testing overidentification, since the estimated instrument always
results in a just-identified system. As a result, Anderson--Rubin intervals
under split-sample approaches will never be empty. Another benefit of our
approach, which
we discuss in \cref{sub:disc} and \cref{asec:hte}, is that MLSS
consistently estimates a convex average of marginal treatment effects
assuming the first-stage is consistent for the conditional mean of the
endogenous treatment on the exogenous instruments.}
Moreover, we should still be concerned with weak instruments: The
first-stage $F$-statistic is only $2.2$
(excluding covariates) and $2.3$ (including covariates).


The sample-splitting estimators based on traditional polynomial expansions rather than machine learning all perform poorly, with out-of-sample $R^2$ close to zero and consequently huge confidence intervals (the point estimates also vary wildly across splits). Overall, the MLSS estimators successfully extract more variation from the instruments than the alternatives, and consequently deliver more statistical precision.

\begin{figure}[tb]
  {\centering
    (a) Excluding covariates
    \includegraphics[width=\textwidth]{exhibits/judge_nocov.pdf}



    \vspace{1em}
    (b) Including covariates

    \includegraphics[width=\textwidth]{exhibits/judge_cov.pdf}}

  Notes: Point estimate and confidence intervals across three sample
  splits, represented by the three horizontal panels. TSLS and TSLS
  (Quadratic) are direct estimates without sample-splitting. Out-of-sample
  $R^2$ of instruments on endogenous treatment in annotation in panel (a).

  \quad \textsf{TSLS} is a standard TSLS estimator without sample
  splitting,
  using the instruments directly from the dataset.
  \textsf{TSLS (Quadratic)} includes second-order terms (but not
  interactions) for the instruments, again without sample-splitting---in
  particular, it results in an empty Anderson--Rubin interval.
  \textsf{LightGBM} and \textsf{RandomForest} are MLSS estimators, where
  LightGBM is an algorithm for gradient boosted trees. Finally,
  \textsf{Linear}, \textsf{Quadratic}, and \textsf{Saturated Quadratic} are
  split-sample estimators with linear regression, quadratic regression
  (without interactions), and quadratic regression with interactions as the estimators for the
  instrument, respectively.


  \caption{IV estimation of the effect of appellate court
  reversal on district court sentencing decisions}
  \label{fig:judge}
\end{figure}


\section{Conclusion}
In this paper, we provide a simple and user-friendly analysis of incorporating
flexible prediction into instrumental variable analysis in a manner that
is familiar to applied researchers. In particular, we document via elementary
techniques that a split-sample IV estimator with machine learning methods as the
first stage inherits classical asymptotic and optimality properties of usual
instrumental regression, requiring only weak conditions governing the
consistency of the first stage prediction. In the presence of covariates, we
also formalize moment conditions for instrumental regression that continues to
leverage nonlinearities in the excluded instrument without creating spurious
identification from the nonlinearities in the included covariates. Leveraging
such nonlinearity in the first stage allows the user to extract more identifying
variation from the instrumental variables and can have the potential of rescuing
seemingly weak instruments into strong ones, as we demonstrate with simulated
data and real data from a criminal justice context. Conventional components
of an instrumental variable analysis, such as identification-robust confidence
sets, extend seamlessly in the presence of a machine learning first stage. We
believe that machine learning in IV settings is a mostly harmless addition to
the empiricist's toolbox.




\newpage


\bibliographystyle{ecca}
\bibliography{main.bib}

\newpage