Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,130 characters · 13 sections · 175 citation commands
Nonparametric rich covariateswithout saturation
JEL codes: C14, C26, C45.
Key words: Average treatment effect; Instrumental variables; Local average treatment effect; Semiparametric estimation.
Linear instrumental variables estimators (LIVEs) have become one of the workhorses of the program evaluation literature, being regularly used to estimate the causal effect of policies with heterogeneous effects (see, e.g., angrist2009, and angrist2010). Indeed, based on the pioneering work of Imbens1994 and angrist1995, LIVEs are generally viewed as identifying the average treatment effect for the sub-population of compliers (the local average treatment effect or LATE), or a weighted average of LATEs; see, e.g., Imbens1994, angrist1995, angrist2000, and mogstad2021.
However, several authors have recently pointed out that, in the standard case in which the model includes covariates, the causal interpretation of estimates obtained using LIVEs requires some strong assumptions. In particular, building on the work of abadie2003 and kolesar2013, blandhol2025 emphasized that, in a just-identified model with a single endogenous explanatory variable and covariates, LIVEs only identify a parameter with a causal interpretation if the expectation of the instrument conditional on the covariates is linear. In the words of blandhol2022, when this condition is satisfied the model is said to have rich covariates.\footnote{The introduction of covariates also has implications for the so-called monotonicity condition; see sloczynski2024 and mogstad2024 for details. This is not the focus of the paper and we assume that the necessary monotonicity condition is satisfied.}
kolesar2013 notes that the rich-covariates condition is uncontroversial when the model is saturated, or when the instrument is unconditionally randomly assigned, and blandhol2025 state that, outside of these two special cases, having rich covariates is a parametric assumption that needs to be defended mogstad2024.
In this paper, we propose methods to ensure that the rich-covariates condition is satisfied without imposing parametric assumptions, and without relying on model saturation or randomly assigned instruments. The estimators we propose make use of a preliminary estimate of the expectation of the instrument conditional on the covariates, which can be obtained nonparametrically. We then consider two ways of using this estimate to ensure that the rich-covariates condition is fulfilled, at least asymptotically, even in cases where the instrument is not unconditionally randomly assigned and the model is not saturated. The two approaches we consider are closely related to the two cases that kolesar2013 identified as ensuring rich covariates.
The first method, which we refer to as the instrument-residual approach following lee2021, uses as the instrument the difference between the original instrument and its estimated conditional expectation. This ensures that the rich-covariates condition is satisfied because the difference between the instrument and its conditional expectation is mean independent of all functions of the covariates, and therefore its conditional expectation is constant (more specifically zero). The second method is a control-function approach in which the estimated conditional expectation of the instrument is added to the set of explanatory variables. This ensures that the rich-covariates condition is satisfied because the conditional expectation of the instrument is trivially a linear function of itself.
The estimators we propose have as special cases other estimators considered in the literature. In particular, if the conditional expectation of the instrument is estimated using a saturated linear model, the two methods we propose lead to results that are numerically identical to those obtained by estimating the saturated model discussed, for example, by angrist1995, kolesar2013, mogstad2024, and blandhol2025. The advantage of our approach is that other nonparametric estimators of the conditional expectation can be used, which is particularly convenient when the set of controls includes continuous covariates or discrete covariates with large support. Naturally, like saturation, nonparametric estimators of the conditional expectation may be affected by the curse of dimensionality, and therefore implementation may be challenging in models with many controls; we discuss this issue in Section (ref).
In turn, our instrument-residual estimator has as a special case the estimators proposed by lee2021 and kim2024; see also lee2024. Our approach differs from that of Lee and co-authors in two important ways. First, we argue that the parametric approaches used by these authors to estimate the conditional expectation of the instrument are far too restrictive, and emphasize the importance of using nonparametric estimators. Second, our approach accommodates an extensive set of functions of the control variables in the estimating equation, potentially leading to substantial efficiency gains.
The instrument-residual estimator is also closely related to the double/debiased machine learning (DML) estimators of the partially linear instrumental variables model considered by chernozhukov2018, which also use an instrument residual.\footnote{LeeandLee25 discuss estimators closely related to those suggested by chernozhukov2018.} Their estimators require the nonparametric estimation of several different objects and involve a computationally-expensive cross-fitting step, therefore being more difficult to implement than our estimators. In Section (ref) we provide further details on the relation between our proposed instrument-residual estimator and those suggested by chernozhukov2018.
Crucially, under standard conditions, both of our estimators identify a positively-weighted average of expected treatment effects for compliers, conditional on the set of covariates. Therefore, the parameter that our estimators\ identify is weakly causal, in the sense of blandhol2025, and has a clear interpretation. Moreover, this estimand is the one identified by the LIVE under rich covariates and coincides with the estimand of the estimators proposed by lee2021 and kim2024, and with the estimand of the DML estimator when there are heterogeneous treatment effects.
In short, the estimators we propose are in between the instrument-residual estimators of Lee and co-authors, which are implemented using strong parametric assumptions and can be very inefficient, and the DML estimators of Chernozhukov and co-authors, which require the nonparametric estimation of more objects and rely on cross-fitting which, as we will show, brings its own problems.
Although with a very different motivation, the two approaches we suggest parallel the ones considered by BH23 in a related but distinct context.\footnote{We are grateful to Julian Costas-Fernandez for pointing to us the parallel between our approaches and those of BH23.} In contrast to what happens in BH23, we cannot use simulation methods to estimate the expectation of the instrument conditional on the covariates, and have to rely on other non-parametric estimators such as kernel regression or machine-learning methods.
The use of a nonparametric estimate of the conditional expectation of the instrument raises the possibility that the estimate of the parameter of interest may not be $\sqrt{n}$-consistent. However, it is well known that estimators can converge at the parametric rate even when they depend on a preliminary nonparametric step that does not converge at the same rate NMcF1994. In Section (ref) we show that, under suitable conditions, a similar result holds in our case when the nonparametric step is performed using kernel methods. It must be emphasized that this result is presented to show that our approaches can lead to estimators that are $\sqrt{n}$-consistent, and does not imply that we recommend the use of kernel regression in the preliminary nonparametric step.
The remainder of the paper is organized as follows. Section (ref) presents the setup and notation. Section (ref) introduces our proposed estimators under the assumption that the expectation of the instrument conditional on the covariates is known, and then considers the properties of the estimators when the expectation of the instrument conditional on the covariates needs to be estimated. In Section (ref) we show that the estimators of the parameter of interest can be $\sqrt{n}$ -consistent when the expectation of the instrument conditional on the covariates is estimated by a suitable kernel regression. Section (ref) presents a small simulation study illustrating the performance of the proposed methods and comparing them with alternative approaches. Section (ref) revisits one of the empirical studies considered by blandhol2025 to illustrate the application of the methods we propose and highlight their attractiveness. Finally, Section (ref) concludes and the appendices contain further technical material and the proofs.
As, for example, in kolesar2013 and blandhol2025, we consider the case where we want to learn about the heterogeneous causal effect of a binary treatment $t\in \left\{ 0,1\right\} $ on the outcome $y$. As usual, $y$ can be written as $y(1)t + y(0)\left( 1-t\right)$, where $y(1)$ and $y(0)$ denote the outcome with and without treatment, respectively. The treatment is potentially endogenous, but we assume that we also observe an instrument $z\in $ $\left\{ 0,1\right\} $ that is exogenous, at least when conditioning on a set of $d$ observable covariates collected in the vector $\boldsymbol{c}$. The treatment status depends on $z$, and we use $t(0)$ and $t(1)$ to denote the potential values of the treatment when $z$ equals $0$ or $1$, respectively.
In this framework, we are interested in the conditions under which a researcher can obtain information on the causal effect of $t$ on $y$ by using instrumental variables to estimate a linear equation of the form
where $\boldsymbol{r}$ is a $k\times 1$ vector containing functions of $\boldsymbol{c}$ ($k\geq 0$), $\boldsymbol{x}\coloneqq (t,\boldsymbol{r}^{\prime })^{\prime }$, $\boldsymbol{\beta }\coloneqq (\alpha ,\boldsymbol{\gamma }^{\prime })^{\prime }$ is the corresponding $\left( k+1\right) \times 1$ vector of parameters, and $\varepsilon $ is an error term. We note that we focus on the LIVE of ((ref)) because this is the approach commonly used by practitioners who want to estimate causal relationships, and not because we assume that such linear equation is an adequate representation of the data generating process. For example, we do not assume that the effect of $t$ on $y$ is constant, or that there is linearity of $y$ in $r$.\footnote{It is worth noting, however, that linearity of $y$ in $r$ is used in the moment condition for the LIVE, so we would expect the estimator to be more precise if linearity in $r$ is a good approximation. }
For the instrumental variable estimate of $\alpha $ to have standard properties and a causal interpretation, the instrument needs to satisfy a number of conditions. In particular, the instrument needs to be exogenous, in the sense that $\left( y\left( 0\right) ,y\left( 1\right) ,t\left(0\right) ,t\left( 1\right) \right) \mathrel{\perp\!\!\!\!\perp}z\mid\boldsymbol{c}$, relevant, in the sense that $\mathbb{C}\text{ov}( z,t\mid\boldsymbol{c}) \neq 0$, and a suitable form of monotonicity needs to hold; see sloczynski2024 and mogstad2024 for details.\footnote{Specifically, as we will see in Section (ref), we need that the sign of $\mathbb{C}\text{ov}( z,t\mid\boldsymbol{c}) $ is the same for all $\boldsymbol{c}$. The validity of this assumption has to be considered on a case-by-case basis.} Here, we assume that those other necessary conditions hold and, as in blandhol2025, focus on the rich-covariates condition, which requires that
where $\mathbb{L}\left[ z\mid\boldsymbol{r}\right] $ denotes the linear projection of $z$ on $\boldsymbol{r}$.\footnote{blandhol2025 write the rich-covariates condition as $\mathbb{E} \left[ z\mid\boldsymbol{c}\right] =\mathbb{L}\left[ z\mid\boldsymbol{c}\right] $ but it is useful to distinguish the set of conditioning variables from the set of functions of those variables that are included in the estimating equation.} kolesar2013 calls ((ref)) the linearity assumption, but we prefer to use the rich-covariates label introduced by blandhol2022 because linearity can be misconstrued to relate to the linearity of ((ref)), to the linearity of the data generation process, or to the linearity of the conditional expectation $\mathbb{E}\left[y(t)\mid\boldsymbol{c}\right] $ for $t\in $ $\left\{ 0,1\right\} $.
This section introduces the two approaches we propose. As, for example, in RMN92 and abadie2003, we start by considering the infeasible \textquotedblleft oracle\textquotedblright\ estimators obtained assuming that $\zeta _{0}(\boldsymbol{c})\coloneqq\mathbb{E}[z\mid\boldsymbol{c}]$ is known, and then we consider the implications of replacing the unknown conditional expectation by its estimated counterpart.
As mentioned before, Koles\'{a}r (2013) noted that the rich-covariates condition is uncontroversial when the instrument is unconditionally randomly assigned. In this case the instrument is statistically independent of $\boldsymbol{c}$, and consequently its conditional expectation is constant. However, the instrument does not need to be statistically independent of the controls for its conditional expectation to be constant, it is enough for the instrument to be mean independent of $\boldsymbol{c}$. Therefore, given $z$ and $\zeta _{0}( \boldsymbol{c}) $, we can ensure that the model has rich covariates by using the instrument residual $z^{\ast}( \zeta _{0}) \coloneqq z-\zeta _{0}( \boldsymbol{c}) $ as the instrument for $t$.
We are not the first to propose an instrument-residual approach. Among others, lee2021 and kim2024 have also proposed estimators that use $z^{\ast }( \zeta _{0}) $ as the instrument BH21. Their estimators can be interpreted as the LIVE of $y$ on $t$ using the instrument residual\ $z^{\ast }( \zeta _{0}) $, in an estimating equation that has no additional controls. For the case where $z$ is binary, Lee (2021) shows that his estimator identifies
where, as in blandhol2025, the subscript “rich” indicates that this is the parameter identified by LIVEs when the rich-covariates condition is satisfied, $\mathbb{E}_{\text{C}}$ denotes that expectations are taken over the sub-population of compliers, and $\omega ( \boldsymbol{c}) $ are the so-called overlap weights defined as
with $\mathbb{C}\text{ov}( z,t\mid\boldsymbol{c}) =\Pr \left( \text{Complier}\mid\boldsymbol{c}\right) \zeta _{0}( \boldsymbol{c}) \left( 1-\zeta_{0}( \boldsymbol{c}) \right) $, where $\Pr \left( \text{Complier}\mid\boldsymbol{c}\right) $ denotes the conditional probability of being a complier.\footnote{sloczynski2020,sloczynski2024 also presents this result.} That is, the estimator gives additional weight to individuals who are more likely to comply, and therefore for whom the instrument is stronger, and to individuals for whom $\boldsymbol{c}$ has a weaker relation with the instrument, i.e., for whom the instrument is closer to being randomly allocated lee2021.
Equation ((ref)) shows that $\alpha _{\text{rich}}$ is weakly causal, in the sense of blandhol2025, as long as the sign of $\mathbb{C}\text{ov}( z,t\mid\boldsymbol{c}) $ is the same for all $\boldsymbol{c}$; that is, as long as the strong monotonicity condition holds sloczynski2024, mogstad2024. kim2024 show that this result can be generalized to the case where the instrument is multi-valued, and note that continuous instruments can be discretized. BH21 address a closely related problem and present the estimand for the corresponding instrumental-residual estimator when the instrument is continuous. Crucially, all these results do not depend on ((ref)) being correctly specified in any reasonable sense.
It is easy to see that adding covariates to the model does not change the probability limit of the estimator for the coefficient on $t$ (see BH21, and Appendix (ref)). Therefore, our instrument-residual estimator has the same estimand as the estimators of lee2021 and kim2024, but our estimator can be much more efficient because the additional variables included in $\boldsymbol{r}$ reduce the variance of the error. Indeed, controlling for additional covariates is often motivated as a way to gain efficiency (see, e.g., abadie2003 and BH23).
kim2024 suggest the parametric estimation of $\zeta _{0}\left( \boldsymbol{c}\right) $, but note that a nonparametric alternative is to use the DML estimators of the partially linear instrumental variables model presented by chernozhukov2018; blandhol2025 also consider this estimator as an alternative to the LIVE.
The estimators introduced by chernozhukov2018 are closely related to our instrument-residual estimator in that they also use an instrument residual. Indeed, chernozhukov2018 consider the estimation of a partially linear instrumental variables model of the form
where $\xi ( \boldsymbol{c}) $ is a function of $\boldsymbol{c}$ and $\mathbb{E}\left[ v\mid\boldsymbol{c},z\right] =0$. The authors consider estimators for $\rho $, the parameter of interest, based on two alternative moment conditions, both of which use $z^{\ast }( \zeta _{0}) $ as the instrument. Their preferred approach is based on the score (see their equation 4.8)
where, besides $\zeta _{0}( \boldsymbol{c}) $, $\mathbb{E}\left[ y\mid\boldsymbol{c}\right] $ and $\mathbb{E}\left[ t\mid\boldsymbol{c}\right] $ have to be estimated nonparametrically. The alternative estimator is based on the score (see their equation 4.7)
which is more difficult to operationalize because $\xi \left( \boldsymbol{c} \right) $ cannot be directly estimated.\footnote{chernozhukov2018 also consider a different score, but that is not based on $z^{\ast }( \zeta _{0}) $.}
Both of these scores can be expressed as
where $m( \boldsymbol{c}) $ is a function of $\boldsymbol{c}$. However, because the instrument residual $z^{\ast }( \zeta _{0}) $ is mean independent of any function of $\boldsymbol{c}$, the form of $m( \boldsymbol{c}) $ does not affect the estimand, but it affects the efficiency of the estimator.\footnote{Indeed, $m( \boldsymbol{c}) $ can even be omitted from the score function, as in lee2021 and kim2024.}
Therefore, our estimator can be seen as a simplification of the estimators proposed by chernozhukov2018, in which $m( \boldsymbol{c}) $ is specified by the practitioner as $\boldsymbol{r}^{\prime } \boldsymbol{\gamma }$. Naturally, chernozhukov2018's (chernozhukov2018) approach has the advantage of controlling for $\boldsymbol{c}$ nonparametrically, and that will often lead to increased efficiency, at the cost of being computationally more expensive.\footnote{Additionally, because it relies on estimates of $\mathbb{E}\left[ y\mid\boldsymbol{c}\right] $ and $\mathbb{E}\left[ t\mid\boldsymbol{c}\right] $, chernozhukov2018's (chernozhukov2018) estimator requires assumptions about these functions, which we do not generally need.}
A key feature of the DML estimators introduced by chernozhukov2018 is the fact that they are based on Neyman-orthogonal scores, which reduces the sensitivity of the estimator of $\rho $ to the noise in the estimates of the nonparametric functions LeeandLee25. We show in Appendix (ref) that the score on which our instrument-residual estimator is based would satisfy the Neyman-orthogonality condition if and only if ((ref)) were correctly specified, in the sense that $\mathbb{E}\left[ \varepsilon \mid \boldsymbol{c}\right] =0$, which is an assumption we do not make.\footnote{Note, however, that because the score implied by our estimator nests the one of chernozhukov2018's (chernozhukov2018) preferred estimator, the condition $\mathbb{E}\left[ \varepsilon \mid \boldsymbol{c}\right] =0$ would be satisfied if $\boldsymbol{r}$ were to include estimates of $\mathbb{E}\left[ y\mid\boldsymbol{c}\right] $ and $\mathbb{E}\left[ t\mid\boldsymbol{c}\right] $, in which case the score on which our estimator is based would be Neyman orthogonal.} Because, in general, the score on which our estimator is based is not Neyman orthogonal, we cannot use the results in chernozhukov2018 to establish its properties. As a consequence, our estimator does not use cross-fitting which, besides being computationally expensive, will be shown to be unsuitable in some empirically-relevant scenarios.
The second situation identified by kolesar2013 where the rich-covariates condition is uncontroversial is in saturated models. This shows that the fulfillment of the rich-covariates condition does not depend on $\boldsymbol{c}$ but on how $\boldsymbol{c}$ is used to specify the equation used to estimate the causal effect of $t$ on $y$. We consider an alternative specification to equation ((ref)) that ensures the model has rich covariates.
Observing that
it follows that the rich-covariates condition is satisfied as long as ((ref)) is augmented by the addition of $\zeta \left( \boldsymbol{c} \right) $ as an explanatory variable. That is, the estimate of interest is obtained by applying a LIVE to the extended equation
where $\boldsymbol{x}^{\ast }\coloneqq(t,\boldsymbol{r}^{\prime },\zeta _{0}( \boldsymbol{c}))^{\prime }$ is a $\left( k+2\right) \times 1$ vector, $ \boldsymbol{\beta }^{\ast }\coloneqq(\alpha ^{\ast },\boldsymbol{\gamma }^{\ast \prime },\phi )^{\prime }$ is the corresponding $\left( k+2\right) \times 1$ vector of unknown parameters, $\varepsilon ^{\ast }$ in an error term, and the other symbols are defined as before. As in the instrument-residual approach, equation ((ref)) defines the moment condition upon which our estimator is based, but is not assumed to be correctly specified.
In Appendix (ref) we show that the estimators of $\alpha^{\ast} $ and $\gamma^{\ast}$ using $z$ as an instrument for $t$ yield the same estimands as the instrument-residual approach BH23. As before, it is helpful if $\boldsymbol{r}$ includes an extensive set of functions of $ \boldsymbol{c}$ because that has the potential to improve the efficiency of the estimator. Additionally, we show in Appendix (ref) that the score on which estimation of $\alpha $ is based is Neyman orthogonal only under very strong conditions that generally do not hold.\footnote{For example, the score will be Neyman-Orthogonal when the pseudo-true value of $\phi$ is $0$ and $\mathbb{E}\left[ \varepsilon ^{\ast}\mid\boldsymbol{c}\right] =0$. Once again, these conditions hold when $\boldsymbol{r}$ is composed of estimates of $\mathbb{E}\left[ y\mid\boldsymbol{c}\right] $ and $\mathbb{E}\left[ t\mid\boldsymbol{c}\right] $.}
Given that the two proposed approaches identify the same parameter, it is natural to ask which of them is preferable in practice. We will address this question later, but for now we note that, in contrast to what happens in the instrument-residual approach, in the control-function approach the instrument is binary, which may be relevant in some circumstances. For example, the control-function approach can be used to ensure rich covariates when using the $\kappa $-weighted regression suggested by abadie2003; see also blandhol2025.
We now discuss the consequences of using different estimators of $\zeta _{0}(\boldsymbol{c})$, consider the effects of this choice on the properties of the estimators of $\alpha $, and motivate the estimators used in the simulation exercise performed in Section (ref) and in the empirical illustration presented in Section (ref). However, we emphasize that the optimal choice of the estimator of $\zeta _{0}(\boldsymbol{c})$ may vary from application to application and that future developments in nonparametric estimation may lead to different recommendations on the choice of method to use in the preliminary step.
lee2021 and kim2024 suggest estimating $\zeta _{0}(\boldsymbol{c})$ by parametric methods, and that is the approach implemented in lee2024. The obvious advantages of this approach are its simplicity and the fact that, as in the classic works by pagan1984 and murphy2002, the estimators of $\alpha $ will converge at the usual $\sqrt{n}$ rate (RMN92, RMN92, and abadie2003, abadie2003, for example, also consider parametric approaches in related contexts). The drawback, however, is that misspecification of the parametric model for $\zeta _{0}(\boldsymbol{c})$ implies that the rich-covariates condition is not fulfilled and there is no guarantee that the estimator of $\alpha $ can be given a causal interpretation. Therefore, except possibly in applications where economic theory provides some guidance on its form, we find it difficult to recommend the use of a parametric estimator of $\zeta _{0}(\boldsymbol{c})$.
The additional robustness afforded by using a nonparametric method to estimate $\zeta _{0}(\boldsymbol{c})$ has non-negligible costs. A first potential drawback is that the resulting estimator of $\alpha $ may not converge at the parametric rate. However, it is well known that under suitable regularity conditions, estimators can be $\sqrt{n}$-consistent even when they depend on a preliminary nonparametric step that does not converge at the same rate (see, e.g., robinson1988, newey1994, newey1994, newey1997, NMcF1994, and newey2004). For example, abadie2003 obtains an estimate of $\zeta _{0}(\boldsymbol{c})$ using power series estimators and, using the results of newey1994, shows that the second step estimator is $\sqrt{n}$-consistent. In the next section we show that, under suitable conditions, a similar result holds in our case when $\zeta _{0}(\boldsymbol{c})$ is estimated using kernel methods. In the simulations in Section (ref) we illustrate the performance of our estimators when the preliminary step is performed in this way.
A related cost of using a nonparametric preliminary step is the well-known curse of dimensionality. In empirical applications, $d$ (the dimension of $\boldsymbol{c}$) is often large and that may severely affect the performance of series and kernel estimators. A possible alternative is to estimate $\zeta _{0}(\boldsymbol{c})$ using machine learning methods such as neural networks, which barron1994 showed can be less affected by the curse of dimensionality than series- and kernel-based estimators. Indeed, recent results suggest that neural network estimators may be able to circumvent the curse of dimensionality under certain conditions bach2017,bauer2019,schmidt2020,kohler2021,braun2024. Neural networks have the additional practical advantage of being able to seamlessly handle both continuous and discrete controls.
The properties of neural networks with different architectures is currently a very active area of research, so it is not yet clear what would be the optimal way to obtain a preliminary estimate of $\zeta _{0}(\boldsymbol{c})$ based on such methods. In the simulations in Section (ref), we study the performance of our estimators when $\zeta _{0}(\boldsymbol{c}) $ is estimated using a neural network based on the rectified linear unit (ReLU) activation function, with a single hidden layer with $100$ nodes, as implemented with the default options in the user-written pystacked Stata command ahrens2023.\footnote{pystacked uses scikit-learn's machine learning algorithms, see scikit-learn.} This kind of neural network shares some of the features of the ones considered by bach2017, and therefore it is likely to have a good performance braun2024.\footnote{Alternatively, we could have used deep networks, such as those considered by schmidt2020, kohler2021, and farrell2021, which also have attractive properties. However, the results in Hornik1989 suggest that the number of hidden nodes is more important than the number of layers. Comparing the performance of alternative network architectures is left for future work.}
An additional advantage of using neural network estimators is that, in contrast with lasso, random forests, and other popular machine learning methods, they do not rely heavily on regularization farrell2021. Specifically, by default, the neural network estimator implemented in pystacked uses $L_2$ regularization with a very small penalty parameter, which is then divided by the sample size. Therefore, these estimates will have minimal regularization bias. This is specially important for the methods that we propose because they are based on moment conditions that generally are not Neyman-orthogonal.
Finally, we note that, just like in the case of series and kernel estimators, it is possible to obtain $\sqrt{n}$-consistent estimators for the parametric part of semiparametric models when neural networks or other machine learning methods are used in a preliminary nonparametric step chen1999,chernozhukov2018, farrell2021. Studying the conditions under which it is possible to obtain a $\sqrt{n}$-consistent estimator of the causal effect of the treatment when $\zeta _{0}(\boldsymbol{c})$ is estimated using a neural network is, however, left for future work.
In this section we show that, under the set of assumptions detailed below, the LIVEs of the causal effect of the treatment in equations ((ref)) and ((ref)) are $\sqrt{n}$-consistent and asymptotically normal when $ \zeta _{0}(\boldsymbol{c})$ is estimated using a suitable kernel regression estimator. The main device we use to achieve the parametric rate for the LIVEs, despite the fact that the first-step kernel estimator has slower rate of convergence, is undersmoothing NMcF1994.
Based on a random sample of size $n$, the conditional expectation $\mathbb{E} \left( z\mid\boldsymbol{c}\right) =\zeta _{0}(\boldsymbol{c})$ can be estimated by the Nadaraya-Watson estimator nadaraya1964,watson1964
where $K_{h}(\boldsymbol{u})\coloneqq h^{-d}K(h^{-1}\boldsymbol{u})$, with $ K(\boldsymbol{u})$ being a kernel whose properties are specified below, and $ h=h(n)$ is a bandwidth term. We make the following assumptions about the kernel, the bandwidth, and the covariates.
Assumption (ref) states that the kernel must be of order $ m $. The first condition in Assumption (ref) controls the variance of $\hat{\zeta}( \boldsymbol{c}) $, while the second condition forces undersmoothing, compared to the MSE-minimizing bandwidth (which is proportional to $n^{-1/(2m+d)}$; see scott2015). Taken together, the two conditions require that $m>d$ (so the kernel must be of higher order as long as $d>1$) and that the bandwidth goes to $0$ at the rate $n^{-a}$, with $ 1/(2m)<a<1/(2d)$. Assumption (ref) is standard in kernel regression analysis. It excludes discrete covariates for simplicity, but extensions of our semiparametric estimators to settings where some or all covariates are discrete are possible Gao2015. Assumption (ref) ensures that the bias of the Nadaraya-Watson estimator decays sufficiently fast when using a kernel of order $m$; this level of smoothness is standard in kernel-based semiparametric estimation and is necessary to achieve $\sqrt{n}$-consistency of the second-step estimator under undersmoothing.
We now consider equation ((ref)) and define the vector $\boldsymbol{q} (\zeta )\coloneqq(z^{\ast }(\zeta ),\boldsymbol{r}^{\prime })^{\prime }$. Also, let $\boldsymbol{q}_{0}\coloneqq\boldsymbol{q}(\zeta _{0})$ and, for $ i=1,\ldots ,n$, $\boldsymbol{q}{_{i}}(\hat{\zeta})\coloneqq\bigl(z_{i}^{\ast }(\hat{\zeta}\left( \boldsymbol{c}_{i}\right) ),\boldsymbol{r}_{i}^{\prime } \bigr)^{\prime }$. Based on a random sample of size $n$, the (two-step) LIVE of $\boldsymbol{\beta }$ is
The moment condition defining $\boldsymbol{\beta }_{0}$ is
and $\boldsymbol{\hat{\beta}}$ solves the sample moment condition
To proceed, we make the following additional assumptions.
Assumption (ref) is a standard integrability condition needed for the asymptotic analysis. Assumption (ref) is an identification condition guaranteeing a unique solution to moment condition ( (ref)). The following result establishes that $\boldsymbol{\hat{ \beta}}$, the LIVE of $\boldsymbol{\beta }_{0}$, is $\sqrt{n}$-consistent and asymptotically normal.
The asymptotic variance $\boldsymbol{G}_{\boldsymbol{\beta }}^{-1} \boldsymbol{\Omega }(\boldsymbol{G}_{\boldsymbol{\beta }}^{\prime })^{-1}$ can be estimated by $\boldsymbol{\hat{G}}_{\boldsymbol{\beta }}^{-1} \boldsymbol{\hat{\Omega}}(\boldsymbol{\hat{G}}_{\boldsymbol{\beta }}^{\prime })^{-1}$, where $\boldsymbol{\hat{G}}_{\boldsymbol{\beta }}\coloneqq\frac{1}{ n}\sum_{i=1}^{n}\boldsymbol{q}_{i}(\hat{\zeta})\boldsymbol{x}_{i}^{\prime }$ and $\boldsymbol{\hat{\Omega}}\coloneqq-\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{ \hat{\tau}}_{i}\boldsymbol{\hat{\tau}}_{i}^{\prime }$, with $\boldsymbol{ \hat{\tau}}_{i}\coloneqq\boldsymbol{q}{_{i}}(\hat{\zeta})(y_{i}-\boldsymbol{x }_{i}^{\prime }\boldsymbol{\hat{\beta}})+\boldsymbol{\hat{\varphi}}_{i}$, where $\boldsymbol{\hat{\varphi}}_{i}\coloneqq(z_{i}^{\ast }(\hat{\zeta} \left( \boldsymbol{c}_{i}\right) )\hat{\varepsilon}_{i},\boldsymbol{0} _{k}^{\prime })^{\prime }$ and $\hat{\varepsilon}_{i}\coloneqq y_{i}- \boldsymbol{x}_{i}^{\prime }\boldsymbol{\hat{\beta}}$.\footnote{ Note that $\boldsymbol{q}(\zeta _{0})\varepsilon +\boldsymbol{\varphi}=
\varepsilon -
\mathbb{E}\left[ \varepsilon \mid\boldsymbol{c}\right] =
$.}
The next result establishes consistency of the variance estimator.
We now consider equation ((ref)) and define $\boldsymbol{q}^{\ast }(\zeta )\coloneqq(z,\boldsymbol{r}^{\prime },\zeta (\boldsymbol{c} ))^{\prime }$. Let $\boldsymbol{q}_{0}^{\ast }\coloneqq\boldsymbol{q}^{\ast }(\zeta _{0})$ and, for $i=1,\ldots ,n$, $\boldsymbol{q}{_{i}^{\ast }}(\hat{ \zeta})\coloneqq\bigl(z_{i},\boldsymbol{r}_{i}^{\prime },z_{i}^{\ast }(\hat{ \zeta}\left( \boldsymbol{c}_{i}\right) )\bigr)^{\prime }$. The (two-step) LIVE of $\boldsymbol{\beta }^{\ast }$ is
The moment condition defining $\boldsymbol{\beta }_{0}^{\ast }$ is
and $\boldsymbol{\hat{\beta}}^{\ast }$ solves the sample moment condition
For the control function approach, we need to replace Assumption (ref) with the following assumption.
The control-function counterpart of Theorem 1 is as follows.
The asymptotic variance $\boldsymbol{G}_{\boldsymbol{\beta ^{\ast }}}^{\ast -1}\boldsymbol{\Omega ^{\ast }(}\boldsymbol{G}_{\boldsymbol{\beta ^{\ast }} }^{\ast \prime })^{-1}$ can be estimated by $\boldsymbol{\hat{G}}_{ \boldsymbol{\beta ^{\ast }}}^{\ast -1}\boldsymbol{\hat{\Omega}^{\ast }(} \boldsymbol{\hat{G}}_{\boldsymbol{\beta ^{\ast }}}^{\ast \prime })^{-1}$, where $\boldsymbol{\hat{G}_{\boldsymbol{\beta ^{\ast }}}^{\ast }}=\frac{1}{n} \sum_{i=1}^{n}\boldsymbol{q_{i}^{\ast }}(\hat{\zeta})\boldsymbol{\hat{x}} _{i}^{\ast \prime }$ and $\boldsymbol{\hat{\Omega}^{\ast }}=-\frac{1}{n} \sum_{i=1}^{n}\boldsymbol{\hat{\tau}^{\ast }}_{i}\boldsymbol{\hat{\tau}} _{i}^{\ast \prime }$, where $\boldsymbol{\hat{\tau}^{\ast }}_{i}\coloneqq \boldsymbol{q}{_{i}^{\ast }}(\hat{\zeta})(y_{i}-\boldsymbol{\hat{x}} _{i}^{\ast \prime }\boldsymbol{\hat{\beta}^{\ast }})+\boldsymbol{\hat{\varphi }_{i}^{\ast }}$, with
and $\hat{\varepsilon}_{i}^{\ast }\coloneqq y_{i}-\boldsymbol{\hat{x}} _{i}^{\ast \prime }\boldsymbol{\hat{\beta}}^{\ast }$, with $\hat{\boldsymbol{ x}}_{i}^{\ast }\coloneqq(t_{i},\boldsymbol{r}_{i}^{\prime },\hat{\zeta}( \boldsymbol{c}_{i}))^{\prime }$.
It is interesting to notice that $\hat{\phi}$ plays an important role in the variance of the estimator. Since the value of $\hat{\phi}$ depends on the specification of $\boldsymbol{r}$, the choice of regressors to include in the equation to be estimated affects the asymptotic variance of the control-function estimator both through changing the variance of the errors and through its influence on the value of $\hat{\phi}$. This contrasts with the instrument-residual estimator where the choice of $\boldsymbol{r}$ only affects the variance of the estimator through the variance of the errors.
In this section we present the results of a small simulation study illustrating the behavior of the proposed approaches, focusing particularly on the effects of the curse of dimensionality. The key feature of the design we use is that the data generating process does not vary as we change the number of control variables used in the estimation. Specifically, in all cases considered, we use $d$ control variables, $c_{i1},\ldots ,c_{id}$, $i=1,\ldots ,n$, each of them drawn independently from a standard normal distribution. However, for all values of $d$ that we consider, these control variables enter the data generating process as $\tilde{c}_{i}\coloneqq \sum_{q=1}^{d}c_{iq}/\sqrt{d}$, which also follows a standard normal distribution. This allows us to generate data with the same distribution for any value of $d$, while allowing the number of control variables to vary. Having generated $c_{i1},\ldots ,c_{id}$ and $\tilde{c}_{i}$, for $i=1,\ldots ,n$, the remaining variables are generated as follows.
The instrument $z_{i}$ is generated as independent draws from a Bernoulli distribution with
where $\Phi \left( \cdot \right) $ denotes the normal cumulative distribution function. That is, for $\psi =0$, $\Pr \left(z_{i}=1\mid c_{i1},\ldots ,c_{id}\right) $ can be consistently estimated by a probit using $c_{i1},\ldots ,c_{id}$ as regressors lee2021,kim2024,lee2024, but for other values of $\psi $ such model will be misspecified.
The treatment indicator is obtained as $t_{i}=\left( 1-z_{i}\right) t_{i}(0)+z_{i}t_{i}(1)$, where $t_{i}(0)$ and $t_{i}(1)$ denote potential treatment indicators and are constructed as
where $\mathds{1}[a]$ is the usual indicator function of the event $a$, $u_{i}\sim \mathcal{N}\left( 0,1\right) $, and $0<\kappa _{\text{AT}},\kappa _{\text{NT}}<1$ denote the share of always takers and never takers, respectively, with $\kappa _{\text{AT}}+\kappa _{\text{NT}}\leq 1$. Note that $t_{i}(0)$ and $t_{i}(1)$ are both equal to $0$ with probability $\kappa _{\text{NT}}$ (the probability of being a never taker), are both $1$ with probability $\kappa _{\text{AT}}$ (the probability of being an always taker), and $t_{i}(1)>t_{i}(0)$ with probability $1-\kappa_{\text{AT}}-\kappa _{\text{NT}}$ (the probability of being a complier); there are no defiers.\footnote{Notice that the probability of being a complier determines the strength of the instrument because $t_{i}=z_{i}$ for compliers.}
Finally, the outcome $y_{i}$ is generated as
where $\eta _{i}\sim \mathcal{N}\left( 0,1\right) $ and
where $\alpha_{\text{NT}}$, $\alpha_{\text{C}}$, and $\alpha_{\text{AT}}$ denote the treatment effects for the subpopulations of never takers, compliers, and always takers, respectively. In this setup, endogeneity is caused by the fact that both $\alpha _{i}$ and $t_{i}$ depend on $u_{i}$. Additionally, there is also noise introduced through $\eta _{i}$ and through the misspecification of the equations being estimated.
Using this data generating process, we performed simulations for $d\in \left\{ 1,3,5,7,9\right\} $, $n\in \left\{ 500,2000,8000\right\} $, $\kappa _{\text{AT}}=\kappa _{\text{NT}}=0.25$, $\alpha _{\text{NT}}=-1$, $\alpha _{\text{C}}=0$, $\alpha _{\text{AT}}=1$, and for $\psi =1$, with new sets of variables being drawn independently for each Monte Carlo replication. In all cases, the goal is to estimate $\alpha _{\text{C}}$, which is the LATE, and the models are estimated using $c_{i1},\ldots ,c_{id}$ as controls.
We start by estimating the regression of $y_{i}$ on $t_{i}$ and $c_{i1},\ldots, c_{id}$ using $z_{i}$ as an instrument for $t_{i}$.\footnote{Therefore, in these simulations, $d=k$. We note that, as discussed before, the performance of the estimators may be improved by changing the functions of $c_{i1},\ldots, c_{id}$ used as explanatory variables in the equation, but we do not explore that possibility.} Additionally, we estimate $\alpha _{\text{C}}$ using lee2024's (lee2024) psr Stata command (with the default options), which implements the estimator proposed by lee2021,\footnote{This is essentially the LIVE of $y_{i}$ on $t_{i}$ only, using a residual instrument based on an estimate of $\mathbb{E}\left[ z_{i}\mid c_{i1},\ldots, c_{id}\right] $ obtained with a probit. Notice that, because we set $\psi =1$, this estimator is based on a misspecified model for $\mathbb{E}\left[ z_{i}\mid c_{i1},\ldots, c_{id}\right] $.} as well as using ahrens2024's (ahrens2024) ddml Stata command, which implements the DML estimator of the partially linear instrumental variables model presented by chernozhukov2018. In light of the discussion in Section (ref), the nonparametric estimations needed to implement the DML are performed using the same neural network estimator we use in the implementation of the methods we propose; see more details below. We also run some simulations in which the nonparametric estimates needed for the DML estimator were obtained using an ensemble method of the type used by blandhol2025, but the performance of the estimator was substantially worse in that case.
For the implementation of the proposed estimators, we considered their oracle version in which the preliminary step uses $\zeta _{0}(\boldsymbol{c}) $, which is given by ((ref)), and two feasible versions based on nonparametric estimates of $\zeta _{0}(\boldsymbol{c})$. Specifically, we estimated $\zeta _{0}(\boldsymbol{c})$ both by using the standard Nadaraya-Watson kernel estimator and a neural network regression. In line with the results in Section (ref), the kernel estimator is implemented using multiplicative Epanechnikov kernels of order $d+1$.\footnote{See Hansen2005 for details on how to construct such kernels.} After some experimentation, the bandwidth for each regressor was set to $(1.1+0.725d)\hat\sigma n^{\frac{-1}{2d+1}}$, where $\hat\sigma$ denotes the estimated standard deviation of the regressor.\footnote{We let the scaling constant increase with the dimension of $\boldsymbol{c}$, and therefore with the order of the kernel, to avoid overfitting when $d$ increases. It may be possible to improve the performance of the estimator by using a different scaling constant, but we did not explore that further. In practice, the scaling constant needs to be chosen on a case-by-case basis.} As mentioned before, the neural network estimator we use is the one implemented by ahrens2023 in their pystacked Stata command, with the default options, which is based on the ReLU activation function, has a single hidden layer with $100$ nodes, and uses the Adam optimizer proposed by kingma2014 for weight optimization.\footnote{We also implemented the proposed estimators based on the neural network regressions using the cross-fitting approach suggested by chernozhukov2018. However, that did not improve the performance of the estimators farrell2021.}
Table 1 summarizes the results for each case and for each estimator. Specifically, the table reports the median of the estimates obtained in $1000 $ replications of the simulation procedure, as well as the interquartile range (in parentheses).\footnote{We report the median and the interquartile range because it is well known that the finite sample distribution of the exactly-identified LIVE has no finite moments, see kinal1980.} Because $\alpha _{\text{C}}=0$, these results can also be interpreted as the median and interquartile range of the bias of the estimators. In the table, the results of the instrument-residual and control-function estimators are labeled Oracle when they are obtained using the true conditional expectations, and labeled NW and NN when the conditional expectations are obtained, respectively, with the Nadaraya-Watson and neural network estimators.
The results in Table 1 show that, as expected, the LIVE does not identify the parameter of interest. The results obtained with lee2021's (lee2021) estimator are more interesting. This estimator is based on a misspecified model for $\zeta _{0}(\boldsymbol{c})$, but the method implemented by lee2024 offers some protection against misspecification of this function and, for $d=1$, this approach has substantially smaller bias than the LIVE. However, as $d$ increases, the effects of the misspecification become evident, and the results obtained with this estimator become similar to those obtained with the LIVE.
The DML estimator performs very well for the larger sample sizes used in these experiments. Indeed, for the five values of $d$ that we consider, DML has essentially no bias when $n\in \left\{2000,8000\right\}$. However, for $n = 500$, DML has noticeable biases, especially for larger values of $d$. The performance of the DML estimator is better in these simulations than in the ones reported by blandhol2025, who note that the DML estimator tends to have a sizable bias that vanishes slowly as the sample size increases. This difference may be a consequence of the different simulation designs, but may also result from the different ways in which the estimator is implemented; as noted before, we also performed experiments in which the nonparametric estimates were obtained using an ensemble method as in blandhol2025, but the performance of the estimator was worse in that case.
Turning now to the results obtained with the estimators proposed in this paper, we start by noting that the performance of the oracle estimators does not depend on $d$, which is to be expected because these estimators are not affected by the curse of dimensionality. Additionally, the two oracle estimators have nearly identical performances, presenting essentially no bias and similar levels of dispersion, which drops with the square root of the sample size. As for the feasible estimators, we find that results obtained with the instrument-residual estimator are generally better than the ones based on the control-function approach. This is especially noticeable for $d$ larger than $1$ because the performance of control-function estimator is much more sensitive to the value of $d$ than that of the instrument-residual estimator. Moreover, we find that the estimators that use the neural network estimate of $\zeta _{0}(\boldsymbol{c})$ (whose results are labeled NN) tend to outperform the estimators based on the Nadaraya-Watson kernel regression (whose results are labeled NW), with the difference being particularly clear in the case of the control-function approach. Finally, we note that, as for the oracle estimators, the dispersion of the estimates obtained with the feasible estimators drops with the square root of the sample size.
Overall, these results suggest that, for small $d$, there is little to choose between the four versions of the estimators we introduced. However, for larger values of $d$, the estimators based on the neural network estimates of the conditional expectation of $z_{i}$ appear to dominate their competitors, with the instrument-residual estimator performing particularly well. Indeed, the results obtained with the instrument-residual estimator that uses a neural network estimator in the preliminary step are remarkably insensitive to the value of $d$ and the performance of this estimator, both in terms of bias and dispersion, is comparable to that of the oracle estimators. This estimator is closely related to the DML estimator and has a similar performance in the larger samples, but generally outperforms it for the smaller samples.
The advantage of the neural network-based estimators may be an artifact of the way the different estimators are implemented or of the simulation design used. However, as discussed before, there are reasons to believe that neural network estimators are able to circumvent the curse of dimensionality bach2017,bauer2019,schmidt2020,kohler2021,braun2024, which may explain their good performance in this exercise. As for the choice between the two approaches that we introduced, the instrument-residual estimator appears to be a safer bet, both because it appears to be less sensitive to the choice to the preliminary nonparametric estimator, but also because it generally has a smaller bias than the control-function estimator.
In this section we illustrate the application of the proposed estimators by revisiting the work of dube2020queens, an empirical application highlighted by blandhol2025 in which both the treatment and instrument are binary.
Using data from 1480 to 1913, dube2020queens investigated whether European states experienced more peace under female leadership. In their study, the outcome of interest is a binary indicator for whether a polity was at war in a given year and the treatment variable is a dummy for whether the polity was ruled by a queen in that year.
dube2020queens consider both just- and over-identified estimators, and here we focus on the just-identified case where the instrument is a binary indicator of whether the previous monarch had a legitimate firstborn male child.\footnote{ The results obtained by dube2020queens for the just- and over-identified cases are broadly similar.} To address concerns about the instrument's validity, dube2020queens control for a set of covariates that includes polity and decade fixed effects, an indicator for whether the previous monarchs were unrelated co-rulers, a binary variable for whether the gender of the firstborn of previous monarchs is missing, an indicator for whether previous monarchs had at least one legitimate child and the birth year is known, and a similar indicator for when the birth year is unknown. In total, there are $64$ control variables, all of which are dummies. Both dube2020queens and blandhol2025 emphasize the importance of controlling for the covariates in this analysis.
Column (1) of Table 2 reproduces the estimated effects (and cluster-robust standard errors) reported by blandhol2025, which differ only slightly from those in the original paper. The difference between the LIVE results with and without covariates confirms the importance of accounting for the role of the controls, while the difference between the least squares and LIVE results suggests the need to account for the endogeneity of the treatment. A remarkable feature of these results is that the DML estimate is qualitatively similar to that obtained with the LIVE with covariates. This is surprising because the results of ramsey1969tests's (ramsey1969tests) RESET test reported by blandhol2025 suggest that the LIVE with covariates does not fulfill the rich-covariates condition, and therefore we would expect it to deliver an estimate reasonably different from the one obtained with DML, because DML is supposed to identify $\alpha_{\text{rich}}$ but the LIVE is not. We will return to this point soon.
Column (2) of Table 2 presents our own estimates, which match those in column (1) for the least squares and the LIVE. For the DML, our estimate is not identical to the one reported by blandhol2025, but the difference is small enough to be entirely justified by the fact that the estimator is based on random sample splitting.\footnote{Like blandhol2025, we report the median of the estimates over $100$ repetitions, form folds based on the clustering indicator, and perform the preliminary step using the same ensemble estimator, but the implementation of the estimator may differ in other aspects not explicitly mentioned by blandhol2025.}
The four bottom rows of Column (2) in Table 2 report results obtained with different estimators in which the relevant non-parametric regressions are performed using the same neural network estimator that we used in the simulations. The estimators used to obtain these results are, respectively, a DML-type estimator that does not use cross-fitting,\footnote{Alternatively, this estimator can be seen as an instrumental residual estimator in which $\boldsymbol{r}$ is composed of estimates of $\mathbb{E}\left[ y\mid \boldsymbol{c}\right] $ and $\mathbb{E}\left[ t\mid \boldsymbol{c}\right] $, also obtained using the shallow neural network estimator.} the instrument-residual estimator proposed earlier but without including covariates in the equation, as in lee2021, and the proposed instrument-residual and control-function estimators that include all covariates in the equation. None of these estimators uses cross-fitting and therefore each regression was estimated only once. The reported standard errors should be seen only as indicative because they do not account for the fact that the estimator relies on a preliminary first step.\footnote{Specifically, we report standard instrumental variables clustered standard errors, using the same $176$ clusters considered by dube2020queens.}
The results obtained with these estimators are strikingly different from the ones obtained before because they have the opposite sign. Except for least squares, the standard errors associated with all these estimates are relatively large, and therefore one may think that the different signs are just the result of sampling noise. However, the fact that there are no sign reversals in each half of the table suggests that there may a deeper reason for this.
Because all covariates in this application are binary, it is possible to obtain information on $\alpha_{\text{rich}}$ by estimating a saturated model, and blandhol2022 report that the estimate obtained with the saturated model is $-0.509$ $(0.523)$. This suggests that, with rich-covariates, we should expect negative estimates of the effect, such as those obtained with the methods we proposed and reported in the bottom part of column (2) of Table 2. Therefore, what is puzzling is that the standard DML estimator, whose results are reported in the fourth line of Table 2, leads to estimates that are closer to the one obtained with the LIVE without rich covariates, than to the one obtained with the saturated model. The explanation for this turns out to be quite simple.
The cross-fitting used in DML accounts for the panel nature of the data by forming sample splits based on the variable used to cluster the standard errors. This ensures that the folds are independent. However, because many of the dummies used as controls are rarely equal to $1$, it may not be possible to identify their effect in some of the sub-samples.\footnote{For example, the mean of the dummy indicator of whether the previous monarchs were unrelated co-rulers is $0.008$ and this variable is only equal to $1$ in three clusters. Therefore, if these clusters are not included in a particular sub-sample, it is not possible to identify the effect of this variable.} This implies that the conditional expectations estimated by cross-fitting only control for covariates with variation in the relevant sub-sample. Consequently, these are not valid estimates of the required conditional expectations, and therefore such estimators do not identify $\alpha _{\text{rich}}$. This problem seriously restricts the practical usefulness of estimators based on cross-fitting in applications involving sparse binary covariates.
Another telltale sign of this problem is that the DML estimates based on cross fitting are very sensitive to the sample split, reflecting the fact that the set of variables effectively controlled for varies with the way the sample is split.\footnote{In the $100$ sample splits we used to obtain the value reported in Table 2, the estimates vary between $ -0.202$ and $2.421$.} To mitigate this, and as in blandhol2025, the DML result we report is the median of the estimates obtained in $100$ repetitions of the estimation process. This increases the computational cost of the estimator very substantially: on a standard desktop computer, estimation with the proposed instrument-residual and control-function methods took less than $30$ seconds, while the DML estimates took over $18$ hours. The immense increase in computational cost, however, does not eliminate the fundamental flaw of this method, which still delivers an unreliable estimate.
The point estimates obtained with the four methods that do not use cross-fitting are relatively close to each other, and are also close the estimate obtained with the saturated model, suggesting that in this example all these methods provide reasonable estimates of $\alpha _{\text{rich}}$. Although, as noted before, the reported standard errors do not account for the variability introduced by the first-step estimation, these results point to potential significant efficiency gains from including the covariates in the equation to be estimated, especially if the estimator controls for the covariates in a flexible way.
In summary, in this application, the proposed instrument-residual and control-function estimators deliver results close to the one obtained with a saturated model, which suggests that both methods are successful in ensuring that the rich-covariates condition is satisfied. Therefore, when it is impossible or impractical to estimate saturated models, the methods we propose may be an appealing alternative; these estimators also do not suffer from the drawback of the DML estimators that this application highlighted. Moreover, the results obtained with the proposed estimators suggest that, in this application, violation of the rich-covariates condition leads to an estimate with the “wrong sign,” a phenomenon emphasized by blandhol2025.\footnote{See their numerical illustration, presented in Subsection 2.3.} Therefore, in contrast with the findings of dube2020queens, our results suggest that states led by queens engaged in war less than those led by kings, although this effect is imprecisely estimated. Using over-identified estimators may help to more accurately identify the effect of queens on war, but we do not pursue that avenue here.
For many years, linear instrumental variables estimators have been one of the main tools used by applied economists to estimate causal relationships. Specifically, guided by the findings of Imbens1994 and angrist1995, practitioners often use instrumental variables estimators with the aim of identifying the average treatment effect on the population of compliers, the so-called LATE. However, the ability of such estimators to identify causal relationships has gradually been called into question, and there is a widespread view that, in models with covariates, linear instrumental variables estimators only identify a causal effect if the model is saturated, in the sense that it includes a dummy variable for each possible combination of the values of the covariates mogstad2024, blandhol2025.
Because saturated models are often impractical or even infeasible, it is interesting to consider alternative approaches to the estimation of causal effects in situations where it is important to account for the role of covariates. lee2021 and kim2024 introduced linear instrumental variables estimators that do not rely on saturation and can identify causal relationships in models with covariates. These estimators depend on a preliminary step that ensures that the resulting estimand has a causal interpretation. However, these estimators are based on a strong parametric assumption that may often be invalid, and therefore are not particularly attractive.
kim2024 and blandhol2025 point out that one way to estimate causal effects without the need to rely on parametric assumptions is to use the double/debiased machine learning estimator of the partially linear instrumental variables model considered by chernozhukov2018; see also LeeandLee25. However, the empirical illustration presented in Section (ref) revealed that this estimator may not be able to ensure the fulfillment of rich-covariates condition when the controls include regressors with little variation, such as dummy variables that are rarely equal to $1$. This is a serious problem that restricts the applicability of the double/debiased machine learning estimator.
We propose two versions of the linear instrumental variables estimator, the instrument-residual and control-function estimators, that retain most of the simplicity of the methods proposed by lee2021 and kim2024 but, like the estimators of chernozhukov2018, do not rely on parametric assumptions. BH23 propose similar estimators to address a different but related problem. Importantly, the estimators we propose have an estimand with a clear causal interpretation, are easy to implement, and are closely related to estimators practitioners are familiar with.
The two estimators that we propose depend on a preliminary step, in which the expectation of the instrument conditional on the covariates is estimated nonparametrically. We show that the estimators of the parameter of interest can converge at the usual parametric rate, even when the estimator in the preliminary step converges at a slower rate. We prove this result for the case where the preliminary step is performed using a kernel regression, but it may be possible to obtain similar results when other nonparametric methods are used in this step.
Traditional methods to estimate conditional expectations, such as kernel and series regressions, can be severely affected by the curse of dimensionality, and therefore may be unsuitable in many practical situations. Indeed, the simulation evidence we present in Section (ref) shows that, when a kernel regression is used in the preliminary step, the performance of the estimators deteriorates quickly with the number of controls, with the control-function estimator being particularly affected by the curse of dimensionality.
The current research on machine learning methods may provide attractive alternatives to the use of traditional nonparametric methods and, as discussed in Section (ref), there is evidence that neural network estimators can offer some robustness to the curse of dimensionality. This is borne out by our simulation results in which the estimators that use a neural network in the first step are remarkably insensitive to the number of covariates in the regression---especially the instrument-residual estimator. The empirical illustrations in Section (ref) also suggest that the method works well in practice.
Given how active the research on machine learning methods is, it would not be prudent to advocate the use of a particular estimation method to perform the preliminary step. However, based on our reading of the literature and on the simulation results we report, we recommend that, at least for now, practitioners perform the estimation of the preliminary step using a shallow neural network based on the ReLU activation function, with a large number of nodes, and using the Adam optimizer for weight optimization.
There are a number of extensions of the methods we propose that would be interesting to explore. Our simulation results suggest that the proposed estimators perform particularly well when the preliminary step uses a neural network, but we only provide results on the asymptotic distribution of the estimators whose preliminary step is performed using kernel regression. Therefore, establishing asymptotic properties for estimators that use other nonparametric methods is a priority. Moreover, we only explicitly considered the case of just identified models with a binary treatment and a binary instrument. Using the results in kim2024 and BH23, our results can easily be extended to other types of instrument and to overidentified models, but it may also be possible to formally extend them to the case where there are multiple treatments. Finally, it would be also very interesting to have more precise guidance on how to perform the nonparametric estimation in the preliminary step, and on the role that machine learning methods can have on reducing the impact of the curse of dimensionality.