Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
80,884 characters · 14 sections · 89 citation commands
Average Marginal Effects in One-Step Partially Linear Instrumental Regressions
\def\spacingset#1{ {#1}} \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\it Keywords:} Instrumental variables, Semiparametric estimation, Reproducing Kernel Hilbert Spaces, Average Marginal Effects, Bootstrap.
\spacingset{1.3}
We consider the partially linear instrumental variables (IV) model
where \(Y \in \mathbb{R}\) is the outcome, \(Z \in \mathbb{R}\) is the endogenous continuous treatment possibly correlated with the error term \(\varepsilon \in \mathbb{R}\), \(X \in \mathbb{R}^{p}\) is a vector of exogenous covariates, and \(W \in \mathbb{R}\) is a continuous instrument. $\beta_0\in\mathbb{R}^p$ is a vector of coefficients, and the treatment or regression function \(h_0\) is nonparametrically specified. We develop an estimation and inference method for the average marginal effect (AME) of the treatment, which is defined as
The AME is a key object in many empirical applications, such as estimating the returns to education or measuring demand elasticities in industrial organization. If the treatment function $h_0$ were linear, the AME of the treatment would equal the slope coefficient of this function, and could be easily estimated by a two-stage least squares regression. However, a linearity or parametric assumption can rarely be justified on economic grounds and entails the risk of misspecification. When linearity is violated, the AME is not consistently estimated by two-stage least squares, and the resulting causal inference or counterfactual analyses can be misleading. Assuming a semiparametric IV model allows to overcome the risk of misspecification and to learn the treatment effects flexibly from the data, while mitigating the curse of dimensionality, see florens2003inverse, florens2012instrumental, and ai2003efficient.
Although non- and semiparametric approaches allow reducing the misspecification risk, inference for a nonparametric IV function is substantially more difficult than inference on a scalar parameter. Moreover, an estimated nonparametric regression remains less easy to interpret than an estimated linear model, where the coefficients can directly be viewed as average marginal effects of the variables and the treatment. For applied work, this is important because researchers often must report results in a form that policymakers and practitioners can interpret. Focusing on the AME in (ref) is appealing for two reasons: (i) it provides an easy-to-interpret index of the treatment effect while addressing misspecification issues; (ii) it enables inference on a scalar parameter rather than on a nonparametric function.
In this paper, we contribute to the literature by providing an inference procedure for the AME in (ref) that has two distinctive features. First, it is based on a single regularization parameter. Second, it builds on a framework widely used in the machine learning and Support Vector Machines literature, namely the Reproducing Kernel Hilbert Space setting.
The first distinctive feature of our procedure is that it is based on a single regularization parameter. Most of the existing literature estimates functionals of $h_0$, such as the AME, by running multi-step regressions. In the first step, they estimate the conditional expectation operator $\operatorname{E}\{\cdot|W, X\}$; in the second step, they use this estimate to obtain an estimator of the function $h_0$, which is then employed to construct an AME estimator. See, e.g., ai2003efficient and ai2007estimation. Such a multi-step regression approach requires selecting multiple regularization parameters, one for each regression. This complicates the practical implementation of these methods, as the fine-tuning of each of these parameters might be difficult in practice. Moreover, each estimation step introduces an estimation error that impacts the finite sample behavior of the final estimator. Unlike existing methods, we provide an estimation and inference procedure for the AME based on a one-step regression. Its main advantage is that it requires a single regularization parameter, which simplifies practical implementation.
The second distinctive feature of our inference method is that it builds on a popular framework from the machine learning literature, the Reproducing Kernel Hilbert Space (RKHS) setting; see, e.g., wainwright2019high, berlinet2011reproducing, and steinwart2008support. Tools and algorithms from this framework have been extensively used in Support Vector Machine and classification methods. The advantages of constructing an inference procedure for the AME based on RKHSs are twofold. First, RKHS tools allow solving difficult computational problems in a tractable manner. In our context, they allow us to provide an estimator of the AME and a test statistic with simple and easy-to-compute expressions. Second, they have excellent finite-sample performances, as we show in our simulations, see Section (ref).
We show that our AME estimator based on RKHSs is asymptotically normal, but the variance of the limiting distribution has a complex analytical form. Thus, we develop an inference procedure relying on the Bayesian bootstrap and establish its validity. Our method is simple to implement and we provide an R package that practitioners can readily use. In addition to its simplicity, we show in simulations that it yields good size control under the null hypothesis and exhibits good power under the alternative in small- and moderate-sample settings. We also illustrate the potential of our method with three empirical applications based on angrist1999using, frankel1999, and sokullu2016semi, with respective samples of 2,024, 150, and 117 observations. These applications reveal the good performance of our estimation and inference procedure on real data, including small samples. Thanks to its computational simplicity, its reliance on a single regularization parameter, and its good performance in finite samples, our method can be readily adopted by practitioners.
\paragraph{Related literature.} Series-based methods for estimating AMEs and, more generally, functionals of IV regressions have been developed in ai2003efficient, ai2007estimation, and santos2011instrumental. In a recent paper, chen2023efficient investigates the properties of several AME estimators based on Artificial Neural Networks, and study their finite sample performance. All these papers are based on a multi-step regression approach and require selecting multiple tuning parameters, one for each regression estimation, which may complicate the practical implementation. breunig2016adaptive considers adaptive estimation of known functionals in nonparametric IV models, but does not study inference on them. chen2018optimal and chen2025adaptive provide inference procedures for known functionals of nonparametric IV regressions. Our context differs from theirs, as in our case the AME is an unknown functional of the nonparametric IV regression, since it involves the expectation of \(h_0'\) with respect to the unknown distribution of \(Z\). beyhum2023one and zhang2023instrumental propose one-step estimators of nonparametric IV models, but do not study inference on the AME. RKHS estimation of nonparametric IV models is also studied in singh2019kernel and zhang2023instrumental. However, they neither consider a partially linear specification nor develop inference for the AME.
\paragraph{Organization of the paper.} In Section (ref), we introduce our AME estimator and derive its expression. The asymptotic behavior of our estimator is established in Section (ref). Section (ref) presents the bootstrap test and obtains its validity. The implementation of our inference procedure as well as the selection of our unique regularization parameter are discussed in Section (ref). In Section (ref), we provide simulation evidence for the small-sample behavior of our bootstrap test, demonstrating its excellent finite-sample performance. Section (ref) contains three empirical applications of our method. Finally, Section (ref) concludes. The Appendix collects the mathematical proofs of our results, auxiliary lemmas, and additional simulations. The R package for implementing our inference procedure can be downloaded at \url{https://github.com/lucasgirardh/rkhsiv}.
We assume that \(h_0 \in \mathcal{H}\), with \(\mathcal{H}\) being a space of functions defined on the support of $Z$. Below we define precisely the space $\mathcal{H}$. For now, let us present the general features of the estimation method. From Equation (ref), we have the moment condition $\operatorname{E}\{Y-h_0(Z) - X^T \beta_0 |W,X\}=0$. We assume that identification holds on $\mathcal{H} \times \mathbb{R}^p$, that is,
Necessary and sufficient conditions for the identification of the partially linear IV model can be found in florens2012instrumental. For example, a sufficient condition for the identification of $(\beta_0, h_0)$ is that (i) $\operatorname{E}\{h(Z)|X,W\}=0$ implies $h=0$, (ii) $\operatorname{E}\{X X^T\}$ is full-rank, and (iii) if $\operatorname{E}\{h(Z)|X,W\}=X^T\beta$ for some $(h,\beta)\in\mathcal{H}\times \mathbb{R}^p$, then $X^T\beta=0$. Part (i) is the classical completeness condition in nonparametric IV models, see newey2003instrumental, d2011completeness, freyberger2017completeness. Part (ii) is the standard absence of perfect collinearity among control variables \(X\). Part (iii) expresses a similar idea: the only way for \(\operatorname{E}\{h(Z)|W, X\}\) to be a linear function of \(X\) is to be the null function. See florens2012instrumental for details. For our purpose, we simply assume the identification condition in Equation (ref).\\ The logic underlying the construction of our estimator follows the ideas in beyhum2023one. Building an estimator based on the conditional moment restriction in (ref) would require the preliminary estimation of the conditional expectation given $(X,W)$, as done by most of the literature. Differently, we consider an equivalent reformulation of (ref) that allows us to avoid such a preliminary estimation. From bierens2016econometric,
where $\bm{i}$ denotes the imaginary root. Let us define
with $\mu$ a probability measure supported on $\mathbb{R}^{p+1}$ having a symmetric characteristic function. Notice that $M(\beta,h)\geq 0$, and, by (ref) and the identification condition (ref), $M(\beta,h)=0\Leftrightarrow (\beta,h)=(\beta_0,h_0)$. Thus,
We will now use the above equation to build an estimator of $(h_0,\beta_0)$. Given an independent and identically distributed (i.i.d.) sample \(\{Y_i, Z_i, X_i, W_i\}_{i=1, \ldots, n}\), we denote by \(\operatorname{E}_n\) the empirical-mean operator, that is, \(\operatorname{E}_n{f(Y_i, Z_i, X_i, W_i)}:= n^{-1} \sum_{i=1}^n f(Y_i, Z_i, X_i, W_i)\) for any function \(f\). To build the empirical counterpart of $M(\beta,h)$, we replace the population expectation \(\operatorname{E}\) with the empirical mean operator \(\operatorname{E}_n\):
where $\mathcal{F}_\mu(\cdot)=\int_{ }\exp(\textbf{i}t^T \cdot)d\mu(t)$ denotes the symmetric characteristic function of $\mu$. \\ Now, minimizing $M_n(\beta,h)$ with respect to $\beta$ and $h$ would lead to overfitting, as $h$ is a nonparametric function. Thus, we regularize the problem by minimizing a penalized version of $M_n$, that is,
where \(\lambda > 0\) is a penalty parameter, and \(\|\cdot\|_{\mathcal{H}}\) is the norm on the space $\mathcal{H}$. In the next section, we discuss the definition of $\mathcal{H}$ and show that the above minimization problem admits a unique solution. The practical choice of the penalty parameter \(\lambda\) is discussed in Section (ref). Given $\widehat h$, we estimate the AME of the endogenous treatment as follows:
We choose the space $\mathcal{H}$ to be a Reproducing Kernel Hilbert Space (RKHS). This allows us (i) to have enough flexibility to recover nonparametrically the treatment function $h_0$, and (ii) to obtain an estimator with a simple closed form expression. For a comprehensive treatment of RKHSs, we refer the reader to, e.g., wainwright2019high, berlinet2011reproducing, and steinwart2008support. To clearly present our estimator, we briefly recall the notion of an RKHS. We denote by $\mathcal{Z}$ the support of $Z$, and let $K:\mathcal{Z}\times \mathcal{Z}\to \mathbb{R}$ be a symmetric function, that is, $K(z_1,z_2)=K(z_2,z_1)$ for any $z_1,z_2\in\mathcal{Z}$. $K$ is said to be positive semidefinite if, for any $m\in\mathbb{N}$ and any collection $z_1,\ldots,z_m\in\mathcal{Z}$, the $m\times m$ matrix with entries $\{K(z_i,z_j):i,j=1\ldots,m\}$ (the associated Gram matrix) is positive semidefinite. A symmetric and positive semidefinite function $K$ is called a kernel. The formal definition of an RKHS of functions is the following.
From steinwart2008support, every symmetric positive semidefinite kernel has a unique RKHS, and every RKHS has a unique reproducing kernel. For the sake of clarity, we provide two examples of RKHSs. For further examples, see wainwright2019high, berlinet2011reproducing, and steinwart2008support. \\
Example 1 (Sobolev RKHS). Let $\mathcal{Z}=[0,1]$, and for any function $f:[0,1]\to \mathbb{R}$, denote with $f^{(\kappa)}$ its $\kappa$th derivative on $(0,1)$. Define
Then, $\mathcal{H}$ is an RKHS of functions with reproducing kernel
where $(a)_+:=\max(a,0)$ for any $a\in\mathbb{R}$. See wainwright2019high. \nopagebreak $\blacksquare$\\
The RKHS of the following example is very popular in the statistics and machine learning literature.\\
Example 2 (Gaussian RKHS). Let $H$ be the space of functions $g:\mathbb{C}\to \mathbb{C}$ such that $g$ is analytic and $\int |g(u)|^2 \exp(|u-\overline{u}|^2/2)du<\infty$, where $\overline{u}$ denotes the complex conjugate of $u$ and $du$ stands for the complex Lebesgue measure on $\mathbb{C}$. Consider the space
where $\text{Real}(g(z))$ is the real part of $g(z)$. Then, $\mathcal{H}$ is an RKHS with reproducing kernel $K(z,u)=\exp(-|z-u|^2/2)$. See steinwart2008support. \nopagebreak $\blacksquare$\\
Having set $\mathcal{H}$ to be an RKHS of functions, we now show that the empirical program (ref) admits a unique solution, and give a closed form expression to our estimator $\widehat \theta$. To this end, let us introduce the following notation. Recall that the characteristic function of the measure $\mu$ is given by \(\mathcal{F}_\mu(\cdot):=\int_{ }\exp(\bm{i} t^T \cdot)d\mu(t)\). We define the following \(n \times n\) matrix via its generic \((i,j)-\)th entry, for \(i, j = 1, \ldots, n\),
We also denote
where $\bm{K}$ is the \(n \times n\) Gram matrix associated with the kernel \(K\). We let \(\bm{I}_n\) be the identity matrix of order \(n\), and we use the shortcut notations for the \(p \times p\) matrix
Let $\mathcal{N}(\bm{K}):=\{\bm{\gamma}\in\mathbb{R}^n:\bm{K}\bm{\gamma}=0\}$ denote the null space of $\bm{K}$ and let $\mathcal{N}(\bm{K})^\perp$ be its orthogonal complement. Formally, $\mathcal{N}(\bm{K})^\perp:=\{\bm{\gamma}_\perp\in\mathbb{R}^n:\bm{\gamma}_\perp^T \bm{\gamma}=0 \text{ for all } \bm{\gamma}\in \mathcal{N}(\bm{K})\}$.
The proof of Proposition (ref) can be found in Appendix (ref). singh2019kernel obtain a similar result for a fully nonparametric IV model without the linear component. Let us make two remarks about the above proposition. First, when the matrix $\bm{K}$ is invertible, the solution to (ref) is unique. Second, from the proof of Proposition (ref), any solution $\bm{\alpha}=(\alpha_1,\ldots,\alpha_n)^T$ to (ref) gives rise to the same function $\sum_{i=1}^n\alpha_i K(\cdot,Z_i)$. Thus, the specific solution to (ref) used to compute $\widehat h$ in (ref) has no effect on $\widehat h$.
As our parameter of interest involves the derivative of \(h_0\), we assume that \(K\) is a differentiable kernel and denote by $K'(\cdot,Z_i)$ the first derivative of $K(\cdot,Z_i)$. From Equation (ref), we estimate the first derivative of $h_0$ by
We then estimate the targeted quantity \(\theta_0\) using Equation (ref), with $\widehat h'$ computed as above. By construction, our estimator of $\theta_0$ is readily obtained by solving the linear system of equations in (ref), and is straightforward to compute. Moreover, it does not rely on multi-step nonparametric estimators, and involves only a single regularization parameter.
To study the asymptotic properties of our estimator, we first reformulate the integral equation (ref) that identifies $(\beta_0,h_0)$. Let $L^2_\mu$ be the space of functions square-integrable with respect to $\mu$, with inner product $\left<g_1,g_2\right>_\mu=\int_{ }g_1(t)\overline{g_2(t)}d \mu(t)$ for $g_1,g_2\in L^2_\mu$, where $\overline{g_2(t)}$ denotes the complex conjugate of $g_2(t)$. Denote by $\|\cdot\|_\mu$ the norm on $L^2_\mu$ induced by $\left<\cdot,\cdot\right>_\mu$, that is, $\|g\|^2_\mu=\left<g,g\right>_\mu$. We define
Below we provide the integrability conditions ensuring that $s\in L^2_\mu$, and that both $\operatorname{A}$ and $\operatorname{B}$ are valued in $L^2_\mu$. Equation (ref) can be reformulated as follows:
Thus, $(\beta_0,h_0)$ can be identified as the unique solution to the following population program
An advantage of this reformulation for our proof is that $\theta_0$ can be expressed in terms of the operators $(\operatorname{A}, \operatorname{B})$. In our proofs, we also consider the empirical (penalized) version of the above minimization program, which allows us to express our estimator $\widehat \theta$ in terms of the empirical counterparts of $\operatorname{A}$ and $\operatorname{B}$. This facilitates the analysis of the asymptotic behavior of $\widehat \theta$.
Part (i) of the above assumption imposes square integrability conditions ensuring that $s\in L^2_\mu$ and that $\operatorname{B}$ is valued in $L^2_\mu$. Parts (ii) and (iii) impose smoothness and boundedness assumptions that are common in the semiparametric literature.
Part (i) of the above assumption formally defines $\mathcal{H}$ and implies that the operator $\operatorname{A}$ is valued in $L^2_\mu$. Part (ii) states that $h_0$ must belong to the RKHS $\mathcal{H}$. The inclusion of $f_Z'$ in $\mathcal{H}$ is needed to handle a bias term that appears in the expansion of $\widehat \theta$.
The following assumption ensures the identification of $(h_0,\beta_0)$. We let \(\mathcal{R}(\operatorname{A}):=\{ \operatorname{A} h\) such that \(h\in\mathcal{H}\}\) and $\mathcal{R}(\operatorname{B}):=\{ \operatorname{B} \beta\text{ such that }\beta\in\mathbb{R}^p\}$ be the ranges of $\operatorname{A}$ and $\operatorname{B}$, respectively.
florens2012instrumental show that Assumption (ref) is a necessary and sufficient condition for identification of the partially linear IV model. Examples of Assumption (ref) can be found in florens2012instrumental. In the lines below Proposition (ref), we discuss how our results can be extended to the case where the operator $\operatorname{A}$ is not injective.
For the following assumption, given any class of functions $\mathcal{G}$, we define its bracketing entropy $N_{[\,]}(\epsilon,\mathcal{G},L^2(P))$ as the minimal number of brackets of $L^2(P)$ size $\epsilon$ required to cover $\mathcal{G}$. For a formal definition, see van2000asymptotic.
Parts (i) and (ii) of the above assumption are needed to obtain the asymptotic stochastic equicontinuity of an empirical process appearing in the expansion of $\widehat \theta$. Part (i) requires the uniform consistency of $\widehat h'$, without imposing a convergence rate. Part (ii) requires that $\widehat h'$ must be included asymptotically in a class of functions with a bounded entropy. This condition is satisfied if $\widehat h$ is twice differentiable and its first two derivatives are bounded in probability. Part (iii) is a boundary condition commonly used in semiparametric frameworks, see, e.g., powell1989semiparametric, escanciano2010testing. It is satisfied when, as we approach the boundaries of $\mathcal{Z}$, $f_Z$ decays to zero and $\widehat h -h_0$ does not explode. In Appendix (ref), we provide mild primitive conditions guaranteeing Assumption (ref). Such primitive conditions impose some smoothness on the kernel $K$ of the RKHS.
To introduce the next condition, we need some additional pieces of notation. Let $\operatorname{P}:L^2_\mu\to \mathcal{R}(\operatorname{B})$ be the projection operator onto $\mathcal{R}(\operatorname{B})$. As the linear operator $\operatorname{B}:\mathbb{R}^p\to L^2_\mu$ is defined on a finite-dimensional space, its range $\mathcal{R}(\operatorname{B})$ is a linear finite-dimensional space (see kreyszig1991introductory). From kreyszig1991introductory, linear finite-dimensional spaces are closed. Thus $\mathcal{R}(\operatorname{B})$ is linear and closed. Since from kress1999linear projection operators onto linear and closed spaces are well defined, $\operatorname{P}$ is a well defined operator. Let $\operatorname{T}:=(\operatorname{I}-\operatorname{P})\operatorname{A}:\mathcal{H}\to L^2_\mu$, where the operator $\operatorname{I}$ is the identity operator, and let us denote with $\operatorname{T}^*$ the Hilbert adjoint of $\operatorname{T}$. By definition, the Hilbert adjoint of $\operatorname{T}$ is the operator $\operatorname{T}^*:L^2_\mu\to \mathcal{H}$ such that $\left<\operatorname{T}\varphi,\psi\right>_{\mu}=\left<\varphi,\operatorname{T}^*\psi\right>_{\mathcal{H}}$ for any $(\varphi,\psi)\in\mathcal{H}\times L^2_\mu$. See kress1999linear. The next assumption is a source condition common in the literature on inverse problems.
The above source condition is common in the nonparametric IV literature and can be interpreted as a smoothness assumption, see, e.g., carrasco2007linear, beyhum2023one, babii2022high, darolles2011nonparametric, florens2012instrumental. It essentially implies a requirement that is necessary for the asymptotic normality of the estimator, see severini2012efficiency. In our proofs establishing the asymptotic behavior of $\widehat \theta$, we employ the above condition to deal with a bias term appearing in the expansion of $\widehat \theta$.
The proof of Proposition (ref) can be found in Appendix (ref). The condition $n \lambda \rightarrow \infty$ guarantees that several variance terms appearing in the expansion of $\sqrt{n}( \widehat \theta - \theta_0)$ are asymptotically negligible, while the condition $n \lambda^2=o(1)$ ensures that a (regularization) bias term vanishes. The result can be extended to accommodate a data-dependent choice of $\lambda$. However, for simplicity, we consider only a deterministic regularization parameter.\\ The asymptotic normality result of Proposition (ref) is obtained under Assumption (ref), which guarantees the identification of $h_0$. Such an assumption requires the injectivity of $\operatorname{A}$, that is, the completeness of the distribution of $Z$ conditional on $(X,W)$. We point out that the convergence in distribution of our statistic can also be obtained without requiring the injectivity of $\operatorname{A}$. Under Assumption (ref), the conditions in severini2012efficiency are satisfied, and the AME $\theta_0$ is identified also if identification of $h_0$ does not hold. We can then extend our results and obtain the asymptotic normality of $\widehat \theta$ by replacing $h_0$ with the projection of $h_0$ onto $\mathcal{N}(\operatorname{T})^\perp$, the orthogonal complement of the null space of $\operatorname{T}$. For simplicity of presentation, we choose not to provide such a more general version of our results, and to work under the completeness condition in Assumption (ref). \\ We finally note that our estimator is based on a moment condition that is not locally robust with respect to the preliminary estimation of $h_0$; see chernozhukov2022locally and bennett2023source. A locally robust approach would require estimating additional correction terms and hence selecting multiple regularization parameters. Differently, our approach relies on a single regularization parameter, which simplifies its implementation.
Although from Proposition (ref) \(\widehat{\theta}\) is asymptotically normal, the variance of the limiting distribution has a rather involved expression. This makes it inconvenient to use the quantiles of the asymptotic distribution for hypothesis testing on \(\theta_0\). Therefore, in the next section, we construct a simple bootstrap test.
We construct a Bayesian bootstrap test for the null hypothesis
with \(\theta_{\mathcal{H}_0} \in \mathbb{R}\) a user-chosen tested value of the AME. For the sake of brevity, we focus on a two-sided test, but our procedure can be easily adapted to implement a one-sided test as well. In view of Proposition (ref), we construct a bootstrap test for $\theta_0$ by bootstrapping the Wald statistic $\sqrt{n}(\widehat \theta - \theta_0)$. To this end, we first define the Bayesian bootstrap versions of $\widehat h$ and $\widehat \beta$, and then we construct the Bayesian bootstrap version of $\widehat \theta$.
Let $\{\xi_i\}_{i=1, \ldots, n}$ be i.i.d. random variables with $\operatorname{E}\{\xi\} = \operatorname{Var}\{\xi\} = 1$, and independent from the sample data. For example, we can sample each $\xi_i$ from an exponential distribution with density $f(\xi)=\exp(-\xi)$. Let $\overline{\xi}:= n^{-1} \sum_{i=1}^n \xi_i$. Note that we use $\overline{\cdot}$ for both complex conjugation and empirical means to avoid introducing additional notation. The distinction will be clear from the context. We define the bootstrap version of $(\widehat{\beta}, \widehat{h})$ as
The above objective function is similar to the one used for $(\widehat{\beta}, \widehat{h})$, with the exception that the $i$th observation has weight $\xi_i/\overline{\xi}$, so each observation is represented in the bootstrap sample with that weight. From the minimization program in (ref), the bootstrapped estimators $\widehat h_b$ and $\widehat \beta_b$ can be computed in the same way as their sample counterparts in Equations (ref) and (ref), provided that the weighting matrix $\bm{F}=\{\mathcal{F}_\mu((X_i^T,W_i)-(X_j^T,W_j)) : i,j=1,\ldots,n\}$ is replaced with
Similarly, we define \(\bm{C}_b:= \bm{X}^T \bm{F}_b \bm{X}\). Thus, with $\bm{\widehat \alpha}_b=(\widehat \alpha_{b,1},\ldots,\widehat \alpha_{b,n})^T$ minimum norm solution to
we compute
and define the bootstrap version of $\widehat \theta$ as
Hence, the bootstrap version of the statistic $\sqrt{n}(\widehat \theta - \theta_0)$ is
We implement the test by using symmetric bootstrap p-values, but our procedure can be easily modified to implement equal-tail bootstrap p-values. To run a test of nominal size $\alpha$, we set the critical value to the $(1-\alpha)$ quantile of the bootstrap distribution of the bootstrapped statistic $|\sqrt{n}(\widehat \theta_b - \widehat \theta)|$. Formally,
where $\Pr_{\xi}$ represents the probability that considers as random only the bootstrap weights $\{\xi_i:i=1,\ldots,n\}$ and as fixed the sample data. We then reject the null hypothesis if and only if $|\sqrt{n}(\widehat \theta - \theta_{\mathcal{H}_0})|>\widehat c_{1-\alpha}$.
To show the validity of the bootstrap test, we introduce some regularity conditions. To avoid introducing additional notation, we make a slight abuse of notation and denote by $\Pr$ the probability treating both the bootstrap weights $\{\xi_i:i=1,\ldots,n\}$ and the sample data $\{(Y_i,Z_i,X_i,W_i):i=1,\ldots,n\}$ as random.
The following assumption is the bootstrap counterpart of Assumption (ref). Mild primitive conditions guaranteeing this assumption are provided in Appendix (ref).
The above result establishes the validity of the bootstrap test, and its proof can be found in Appendix (ref). Part (i) obtains the bootstrap influence function representation of $\widehat \theta_b$, showing that such a representation is a reweighed version of the expansion of $\widehat \theta$ from Proposition (ref). The new weights are the bootstrap weights $\{\xi_i/\overline{\xi}-1\;:\;i=1,\ldots,n\}$. Such a reweighting ensures that the leading term of the bootstrap expansion estimates consistently the distribution of the statistic under the null, while remaining bounded in probability under the alternative. Parts (ii) and (iii) of Proposition (ref) establish that the test has the correct size asymptotically and is consistent, with power converging to one when the sample size goes to infinity.
This section describes the implementation of our method. We first discuss the choice of our single regularization parameter. We select $\lambda$ by minimizing a cross-validated version of the objective function in (ref). Specifically, let us focus on a two-fold cross-validation. We partition the sample into two folds, denoted $S_1$ and $S_2$. We first use $S_1$ for estimation and $S_2$ for validation, and then interchange their roles. Let $\widehat \beta_{S_1,\lambda}$ and $\widehat h_{S_1,\lambda}$ denote the estimators obtained from the first fold for a given value of $\lambda$. Define $\widehat \beta_{S_2,\lambda}$ and $\widehat h_{S_2,\lambda}$ analogously for the second fold. The cross-validated version of (ref) is given by
where $\sum_{i\in S_{j}}$ denotes the sum over observations in fold $S_j$, for $j=1,2$. Finally, we minimize the above criterion with respect to $\lambda$, and we take the minimizer as the selected regularization parameter.\\ Compared to the usual out-of-sample mean squared error criterion, the above cross-validation objective employs the appropriate weights, those corresponding to our objective function \(M_n(\beta, h)\) in Equation (ref). In our simulations and empirical applications, we implement the above two-fold cross-validation method with good results. \\ To reduce the computational burden of the bootstrap procedure, we fix the regularization parameter at the value selected from the original sample and use this same value across all bootstrap replications. As shown in our simulation study, this is enough to provide a good performance of our procedure. Algorithm (ref) summarizes the practical steps required to implement our test.
In this section, we study the finite-sample performance of our test in two settings. We first consider a fully nonparametric model corresponding to Equation (ref) without \(X\). We then consider the general case of a partially linear model, where we also include the regressor $X$.
\paragraph{Fully nonparametric model.} We use a data-generating process consistent with model (ref) without \(X\), where we draw the variables \(Z\), \(W\), and \(\varepsilon\) as follows:
with \((W, V, U)\) mutually independent standard Gaussians. In this setup, the marginal distributions of $\varepsilon$ and $Z$ are also standard Gaussians. For a non-zero value of $a$, the regressor $Z$ is endogenous and correlated with the error $\varepsilon$. The parameter $\rho_{\varepsilon V}$ measures the degree of such correlation, and hence the level of endogeneity of $Z$. For a non-null value of $b$, $W$ is a valid instrument for $Z$.
We implement our estimator with $K(z,u)=\exp(-|z-u|^2/2)$, which corresponds to a Gaussian kernel. This is a well-established choice for kernel ridge regression in the machine learning literature (see steinwart2008support) and enables the approximation of a wide range of functions. Specifically, it has a universal approximation property, in the sense that any continuous function can be arbitrarily well approximated by a function in the RKHS with Gaussian reproducing kernel (see steinwart2008support). We set $\mu$ to a Laplace distribution with mean zero and unit variance, thus its characteristic function $\mathcal{F}_\mu$ is a Cauchy density. We choose the penalty parameter $\lambda$ by the two-fold cross-validation described in Section (ref), and perform a grid search to minimize the cross-validated objective function. For bootstrapping, we draw the bootstrap weights $\xi$ from an Exponential distribution with expectation one.
In this nonparametric setting, we compare our procedure to an inference method for the average derivative based on two-stage series regressions ai2007estimation. The method is implemented as follows. We set the series basis to B-Splines and select the number of series terms in the first and second stages using the method proposed in \citet*{chen2025adaptive}. Given the two-stage series estimator of the IV regression, we compute the AME similarly to Equation (ref). We then obtain the p-values by pairwise bootstrap. Specifically, we resample with replacement from the original sample $\{Y_i,Z_i,W_i\}_{i=1}^n$ and, for each bootstrap sample, compute the AME, keeping the same number of series terms selected in the original sample.
We consider two functional forms for the treatment effects: $h_{0,1}(z)=z^2/\sqrt{2}$, corresponding to a quadratic function, and $h_{0,2}(z)=\sqrt{3 \sqrt{3}}\exp(-z^2/2)$, corresponding to a non-polynomial function. Both expressions are normalized so that $h_{0,1}(Z)$ and $h_{0,2}(Z)$ have unit variance. With these functions and the marginal distribution of $Z$, the corresponding true values of the AMEs (\(\theta_0\)) are, respectively, 0 in the quadratic case, and approximately 0.81 in the non-polynomial case. Finally, we set $\rho_{ZW} = 0.8$ and consider two levels of endogeneity: $\rho_{\varepsilon V} = 0.5$, corresponding to a moderate level of endogeneity, and $\rho_{\varepsilon V} = 0.8$, corresponding to a high level of endogeneity.
We ran 5,000 simulations with the sample sizes of $n=100, 400$. To speed up computations, we use the warp-speed method of davidson2007improving and giacomini2013warp. Thus, for each simulated sample, we draw one bootstrap sample and use the collection of bootstrapped statistics to compute the bootstrap p-value for each original statistic. Table (ref) reports the rejection probabilities under the null hypothesis at the usual nominal levels of 5% and 10% for the different functions \(h_0\) and levels of endogeneity \(\rho_{\varepsilon V}\).\footnote{ We also consider the nominal level of 1%, but the rejection rates for the two-stage series regression method are very close to 0. } We label our method “one-step” and the two-stage series regression method “two-step”. Under the null hypothesis, both methods control Type I error, with rejection probabilities below the nominal levels. However, our method delivers rejection rates closer to the targeted nominal levels in most configurations -- it is systematically the case with \(n = 100\) and with the non-polynomial function \(h_{0,2}\). The improvement is notably large for the smaller sample size \(n = 100\).
We also study the behavior of our test under the alternative. To do so, we compute the rejection rates of the test of the hypothesis \(\theta_{\mathcal{H}_0} = \theta_0 + \gamma\), for a deviation \(\gamma \in (0, 1)\) from the null, using 5,000 Monte Carlo replications. In Figures (ref) and (ref), we compare our one-step regression test (displayed in blue) with the two-step series regression test (displayed in red). Figure (ref) presents the results for the quadratic function $h_{0,1}$ and a moderate level of endogeneity, while Figure (ref) presents the results for the non-polynomial function $h_{0,2}$ and a high level of endogeneity. The other settings display similar results and are reported in Appendix (ref), Figures (ref) and (ref). They show that our method has considerably higher power for the small sample size \(n = 100\). For \(n = 400\), both methods have comparable power: the two-step test slightly outperforms in the quadratic case \(h_{0, 1}\), whereas ours outperforms in the non-polynomial case \(h_{0, 2}\).
Overall, when considering both size and power, our test performs best in almost all cases, often substantially outperforming the two-step series regression test, while in the few scenarios where it is slightly outperformed, the difference is minor.
\paragraph{Partially linear model.}
We now study the finite-sample performance of our test in the general partially linear setting corresponding to Equation (ref). Similarly to the fully nonparametric model, ($Z$, $\varepsilon$) are drawn from Equation (ref). $U$ and $V$ are mutually independent standard Gaussians, and are jointly independent of the covariate $X$ and the instrument $W$. The pair $(X,W)$ is drawn from a zero-mean bivariate Gaussian, with $\operatorname{Var}(X)=\operatorname{Var}(W)=1$, and a correlation coefficient between $X$ and $W$ equal to 0.5. In the simulations, we set the slope coefficient of \(X\) to 1. We consider the same two scenarios as in the fully nonparametric model, depending on \(\rho_{\varepsilon V}\) to monitor the magnitude of the endogeneity, as well as the two cases of a quadratic and a non-polynomial functional form for \(h_0\). Our test is implemented as described in the fully nonparametric case.
Table (ref) reports the rejection rates of our test under the null hypothesis \(\theta_{\mathcal{H}_0} = \theta_0\) at the nominal levels of 5% and 10%, based on 5,000 Monte Carlo replications. The results show that our test controls Type I error with rejection rates below the nominal levels and approaching them when the sample size increases. Like in the fully nonparametric model, we also study the behavior of our test under the alternative hypothesis by testing \(\theta_{\mathcal{H}_0} = \theta_0 + \gamma\), with \(\gamma\) ranging between 0 and 1. Figure (ref) considers the quadratic case (\(h_{0, 1}\)) with a moderate level of endogeneity (\(\rho_{\varepsilon V} = 0.5\)). The results are qualitatively the same for the other settings reported in Appendix (ref) (Figures (ref) to (ref)). The power curves illustrate the consistency of our test. Overall, these simulations indicate that the solid performance of our test in the fully nonparametric model is preserved when the model is extended to a partially linear specification with additional exogenous covariates.
We apply our method to three empirical datasets to illustrate its potential.
Our first application revisits the seminal study by angrist1999using on the effects of class size on student achievement. We briefly recall the setting and refer to the original article for further details. The data come from a national testing program conducted in Israeli elementary schools in June 1991, with test scores in mathematics and reading comprehension scaled from 0 to 100. angrist1999using exploits Maimonides' Rule -- an administrative regulation that caps class sizes at 40 students -- to construct the predicted class size, \(W = \text{Enrollment} / (\lfloor(\text{Enrollment} - 1)/40\rfloor + 1)\). The discontinuities at enrollment multiples of 40 provide exogenous variation in class size, making \(W\) a valid instrument for the actual class size, denoted \(Z\) in our notation. We focus on fifth-grade students, using a sample of \(n = 2{,}024\) classes. Following the original study, we also consider an exogenous covariate \(X\), the percentage of disadvantaged students (PD index). The outcome variable \(Y\) is the average grade across students for a given class, either in mathematics or reading.
We estimate the partially linear IV model (ref). Our parameter of interest is the AME \(\theta_0 = \operatorname{E}[h'_0(Z)]\), which here measures the average impact of a one-student increase in class size on test scores. In comparison, we recall that a classical application of two-stage least squares assumes a linear relationship between average test score \(Y\) and class size \(Z\), so that the AME equals the corresponding coefficient. We implement our RKHS-based method as described in Sections (ref) and (ref). Inference is conducted using the Bayesian bootstrap test with \(B = 499\) replications. To report the confidence intervals, we invert the bootstrap test of Algorithm (ref). Specifically, for any \(\theta\), we do not reject the null hypothesis {\(\theta_{0} = \theta\)} if and only if \(|\widehat{\theta} - \theta| \leq \widehat{q}_{1-\alpha}\), with \(\widehat{q}_{1-\alpha}\) the \(1 - \alpha\) quantile of the bootstrap distribution of \(|\widehat{\theta}_b - \widehat{\theta}|\). This leads to the interval \([\widehat{\theta} \pm \widehat{q}_{1-\alpha}]\) that uses the quantiles of the bootstrapped statistic \(|\widehat{\theta}_b - \widehat{\theta}|\) as critical values.
Table (ref) presents the results. Applying our method, we find that the effects of class size on test scores are not statistically significant in either mathematics or reading. This finding contrasts with the statistically significant negative effects estimated by two-stage least squares in the original study: \(-0.372\) for reading and \(-0.355\) for mathematics (see Table IV of angrist1999using). Our results are consistent with horowitz2011applied, Section 5.2, who notes that “the [fully] nonparametric [IV] model does not support the conclusion drawn from the linear model that increases in class size are associated with decreased test scores” (page 374). Our results confirm the findings of horowitz2011applied, which are based on an informal bootstrap inference procedure for the nonparametric IV model. He notes that “the [fully] nonparametric [IV] model does not support the conclusion drawn from the linear model that increases in class size are associated with decreased test scores” (page 374). Importantly, our analysis shows that, absent strong functional form assumptions, the data provide limited information about the precise nature of class-size effects. This highlights the value of our inference procedure, which properly accounts for uncertainty in a flexible semiparametric model with theoretical guarantees, while remaining computationally simple through a single regularization parameter.
The previous application uses samples of about 2,000 observations. We now turn to two applications with smaller sample sizes to illustrate the small-sample performance of our method. As we show below, our procedure proves able to detect significant effects while allowing a flexible semiparametric relationship even with small samples.
First, we revisit frankel1999, who investigate whether international trade causes economic growth using a sample of \(n = 150\) countries. The endogenous treatment, the trade share \(Z\) of a given country, is measured as imports plus exports divided by GDP. frankel1999 address the endogeneity issue by a constructed trade share based on geographic factors (bilateral distances, country sizes, and whether countries share borders or are landlocked). We refer to the original article for further details. Here, we use their data and instrumental variable strategy while allowing for a flexible semiparametric relationship between trade and income. More precisely, we estimate the following partially linear model
where the outcome $Y_i$ is GDP per capita in country $i$, the endogenous treatment $Z_i$ is the trade share, the exogenous covariates are population (\(N_i\)) and area (\(A_i\)), and \(W_i\) is the constructed geographic component of trade obtained from the original study. Equation (ref) corresponds to the main specification of frankel1999, with the difference that the trade share enters nonlinearly through \(h_0\), which is nonparametrically specified. Our application thus relaxes the linearity assumption, allowing for a flexible relationship between trade and income.
Table (ref) presents the estimation and inference results for the AME using our method. We find a positive and statistically significant effect of the trade share on GDP per capita, which illustrates that our method can detect effects even with a few hundred observations. The resulting estimate implies that a one-percentage-point increase in trade share raises GDP per capita by 1.15 percent. In comparison, frankel1999 find that a one-percentage-point increase raises income per capita by 1.97 percent (see their Table 3, column 2, page 387). This suggests that allowing for nonlinearities in the relationship changes the estimate, reducing the effects of trade on income, although the effects remain of the same order of magnitude and statistically significant. Based on frankel1999 data and their identification strategy, our semiparametric method thus confirms a positive causal effect of trade on income while relaxing the linearity assumption. This application illustrates the robustness of our flexible inference procedure, which avoids imposing functional form restrictions while delivering reliable inference even in small samples.
Our third application is based on sokullu2016semi, who proposes an empirical semiparametric model for two-sided markets applied to local daily newspapers in the US. Two-sided markets (where a platform serves two distinct groups of users who benefit from each other's participation) are a central object of study in empirical industrial organization. Examples include newspaper platforms, connecting readers and advertisers rysman2004competition, chandra2009mergers, and payment card networks, connecting merchants and cardholders rysman2007empirical; see, e.g., rochet2006two or armstrong2006competition for theoretical contributions and sanchez2021multisided for a recent review. A key feature of these markets is the presence of indirect network effects: the value that agents on one side derive from the platform depends on the level of participation on the other side. In the newspaper context, advertisers value newspapers with larger readership, and readers may be attracted or deterred by the volume of advertising. Correctly estimating these network effects is important for competition policy analysis.
Several empirical studies of two-sided markets have estimated network effects under linearity assumptions. sokullu2016semi relaxes this assumption by specifying the network effect functions nonparametrically. Specifically, she estimates a nonparametric IV model with Tikhonov regularization on a sample of $n = 117$ monopoly newspaper markets. The main finding is that network effects are neither linear nor monotonic: readers' benefits initially increase with advertising but decrease beyond a threshold, and similarly for advertisers with respect to circulation.
Here, to illustrate our method, we consider the readers' side and estimate the demand function. sokullu2016semi uses the results from the nonparametric IV model to select a parametric specification and re-estimates the model. On the readers' side, the nonparametric estimates suggest approximating the network effect of advertising on readers by a third-order polynomial, leading to the following equation
where all variables are defined at the level of each local newspaper market: $N^r$ is the share of readers buying the newspaper; $N^a$ is the share of advertisements in the newspaper; $P^r$ is the daily cover price; and $U$ is an error term capturing unobserved newspaper characteristics affecting readers' demand. The left-hand side is the log-odds transformation of the readers' market share, which arises from the parametric specification of the distribution of readers' and advertisers' benefits (see Section 4.3.1 in sokullu2016semi). The polynomial $\alpha_1 N^a + \alpha_2 (N^a)^2 + \alpha_3 (N^a)^3$ captures the indirect network effect of advertising on reader demand: it measures how the presence of advertisers on the platform affects readers' propensity to buy the newspaper. We apply our RKHS-based method to nonparametrically estimate this network effect by considering
Compared to (ref), the advertising share $N^a$ enters through the unknown function $h_0$, which we estimate nonparametrically while maintaining the linear specification for the price. Our goal is to estimate and perform inference on the average marginal effect $\theta_0 = \operatorname{E}\{h_0'(N^a)\}$, which summarizes how an increase in advertising share affects readers' demand on average. We provide a formal inference procedure for the AME, complementing sokullu2016semi's nonparametric estimation results.
Both $N^a$ and $P^r$ are endogenous in this model: the advertising share is jointly determined with readership on the platform, and the cover price is set by the newspaper in response to unobserved demand characteristics. Our method can be straightforwardly adapted to handle such cases with multiple endogenous regressors. We require two instruments $(W_1, W_2)$ and adapt the matrix $\boldsymbol{F}$ as $[\mathcal{F}_\mu((W_{1i}, W_{2i})^\top - (W_{1j}, W_{2j})^\top) : i, j = 1, \ldots, n]$. Following sokullu2016semi, we use the area of the city (in square miles) as $W_1$ and the average wages in the printing industry at the county level as $W_2$. The first instrument is correlated with advertising rates but plausibly independent of unobserved reader demand characteristics; the second is an income-related variable correlated with circulation but independent of unobserved newspaper characteristics in the advertising equation.
Table (ref) presents the results obtained from our method, implemented as described in the previous sections. As a comparison, sokullu2016semi obtains the following estimates for (ref): $\widehat{\alpha}_1 = -48.34$, $\widehat{\alpha}_2 = 96.75$, and $\widehat{\alpha}_3 = -62.06$. Combined with the descriptive statistics for $N^a$ reported in Table 1 of sokullu2016semi, the average marginal effect of $N^a$ implied by the cubic polynomial specification is thus $ \widehat{\alpha}_1 + 2 \, \widehat{\alpha}_2 \, \operatorname{E}_n\{N^a\} + 3 \, \widehat{\alpha}_3 \, \operatorname{E}_n\{(N^a)^2\} = -8.19. $ In contrast, our nonparametrically estimated AME is $-5.53$ and is statistically significant. The estimate is smaller, although both remain of the same order of magnitude; in particular, $-8.19$ falls within our confidence intervals. Recall that the cubic specification in (ref) is itself guided by sokullu2016semi's prior nonparametric estimation. The similarity between our AME estimate and that implied by the cubic specification therefore yields supporting evidence for the validity of our method.
In conclusion, the three applications presented in this section, spanning sample sizes from 117 to over 2,000 observations, illustrate the value of our method: it provides a formal inference procedure for the AME that delivers economically meaningful results and maintains computational simplicity with a single regularization parameter.
We propose a bootstrap inference procedure for the average marginal effect of an endogenous treatment in the partially linear IV model. Our procedure builds on the RKHS framework and requires a single regularization parameter, which simplifies its practical implementation relative to existing approaches. Our approach can be extended to consider linear functionals of nonparametric IV regressions. In our simulations, we show that our procedure has good size control and power in small and moderate samples, while three empirical applications illustrate its effectiveness on real datasets.
The Appendix is divided into four sections. In Appendix (ref), we provide additional simulation results. Appendix (ref) contains the proofs of our main results, namely Propositions (ref), (ref), and (ref). In Appendix (ref), we collect the auxiliary lemmas and their proofs. Finally, Appendix (ref) contains primitive conditions that satisfy the high-level assumptions.