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.
44,443 characters · 9 sections · 23 citation commands
Minimax Estimation of Conditional Moment Models
\etocdepthtag.toc{mtchapter} \etocsettagdepth{mtchapter}{subsection} \etocsettagdepth{mtappendix}{none}
\begingroup \footnote{A very preliminary version of this work appeared as Adversarial Generalized Method of Moments (see \url{https://arxiv.org/abs/1803.07164})} \addtocounter{footnote}{-1} \endgroup
Understanding how policy choices affect social systems requires an understanding of the underlying causal relationships between them. To measure these causal relationships, social scientists look to either field experiments, or quasi-experimental variation in observational data. Most observational studies rely on assumptions that can be formalized in moment conditions. This is the basis of the estimation approach known as generalized method of moments (GMM) Hansen1982.
While GMM is an incredibly flexible estimation approach, it suffers from some drawbacks. The underlying independence (randomization) assumptions often imply an infinite number of moment conditions. Imposing all of them is infeasible with finite data, but it is hard to know which ones to select. For some special cases, asymptotic theory provides some guidance, but it is not clear that this guidance translates well when the data is finite and/or the models are non-parametric. Given the increasing availability of data and new machine learning approaches, researchers and data scientists may want to apply adaptive non-parametric learners such as reproducing kernel Hilbert spaces, high-dimensional regularized linear models, neural networks and random forests to these GMM estimation problems, but this requires a way of finding solutions to the moment conditions within complex hypothesis classes imposed by the learner and selecting moment conditions that are adapted to the hypothesis class of the learner.
Most recent theoretical developments in machine learning and high-dimensional statistics are founded on statistical learning theory: formulate a loss function (typically strongly convex with respect to the output of the hypothesis), whose minimizer over the hypothesis space is the desired solution; typically referred to as an $M$-estimator. Being able to frame the problem as an $M$-estimation problem with a strongly convex function, leads to many desirable properties : i) tight generalization bounds and mean squared error rates based on localized notions of statistical complexity can be invoked to provide tight and fast finite sample rates with minimal assumptions bartlett2005local,wainwright2019high, ii) regularization can be invoked to make the estimation adaptive to the complexity of the true hypothesis space, without knowledge of that complexity Lecue2018,lecue2017regularization,negahban2012, iii) the computational problem can be typically efficiently solved via first order methods that can scale massively agarwal2014reliable,rahimi2008random,le2013building,sra2012optimization,Bottou2007. This formulation is seemingly at odds with the method of moments language, as many times the moment conditions do not correspond to the gradient of some loss function and this problem is exacerbated in the case of non-parametric endogenous regression problems (i.e. when the instruments in the observational study does not coincide with the treatments). The problem is: Can we develop an analogue of modern statistical learning theory of $M$-estimators, for non-parametric problems defined via moment restrictions?
Our starting point is a set of conditional moment restrictions:
where $y$ is an outcome of interest, $x$ is a vector of treatments and $z$ is a vector of instruments.
To obtain a criterion function, we first move to an unconditional moment formulation, where the moment restrictions are products of the moment conditions and test functions in the instruments. We then take as our criterion function the maximum moment deviation over the set of test functions, where the set of test functions is potentially infinite.
We show that as long as the set of test functions ${\mathcal F}$ contains all functions of the form $f(z) = \mathbb{E}[h(x) - h'(x) \mid z]$ for $h,h'\in {\mathcal H}$, then such an estimator achieves a projected MSE rate that scales with the critical radius of the function classes ${\mathcal F}$, ${\mathcal H}$ and their tensor product class (i.e. functions of the form $f(z)\cdot h(x)$, with $f\in {\mathcal F}$ and $h\in {\mathcal H}$). The critical radius captures information theoretically optimal rates for many function classes of interest and thereby this main theorem can be used to derive tight estimation rates for many hypothesis spaces. Moreover, if the regularization terms relate to the squared norms of $h, f$ in their corresponding spaces, then the estimation error scales with the norm of the true hypothesis, without knowledge of this norm.
We offer several applications of our main theorems for several hypothesis spaces of practical interest, such as reproducing kernel Hilbert spaces (RKHS), sparse linear functions, functions defined via shape restrictions, neural networks and random forests. For many of these estimators, we offer optimization algorithms with performance guarantees. As we illustrate in extensive simulation studies, different estimators are best in different regimes.
\paragraph{Related work} The non-parametric IV problem has a long history in econometrics newey2003instrumental,blundell2007semi,chen2012estimation,chen2018optimal,hall2005nonparametric,horowitz2007asymptotic,horowitz2011applied,darolles2011nonparametric,chen2009efficient. Arguably the closest to our work is that of chen2012estimation, who consider estimation of non-parametric function classes and estimation via the method of sieves and a penalized minimum distance estimator of the form: $\min_{h \in {\mathcal H}} \mathbb{E}[\mathbb{E}[y-h(x)\mid z]^2] + \lambda R(h)$, where $R(h)$ is a regularizer. As we show in (ref), our estimator can be interpreted asymptotically as a minimum distance estimator, albeit our estimation method applies to arbitrary function classes and non just linear sieves. There is also a growing body of work in the machine learning literature on the non-parametric instrumental variable regression problem deepiv,bennett2019deep,singh2019kernel,muandet2019dual,muandet2020kernel. Our work has several features that draw connections to each of these works, e.g. bennett2019deep,muandet2019dual,muandet2020kernel also use a minimax criterion and bennett2019deep,muandet2019dual also impose some form of variance penalty on the test function. We discuss subtle differences in (ref). Moreover, singh2019kernel,muandet2019dual also study RKHS hypothesis spaces and deepiv,bennett2019deep also study neural net hypothesis spaces. None of these prior works provide finite sample estimation error rates for arbitrary hypothesis spaces and typically only show consistency for the particular hypothesis space analyzed (with the exception of singh2019kernel, who provide finite sample rates for RKHS spaces, under further conditions on the smoothness of the true hypothesis). In (ref) we offer a more detailed exposition on the related work and how it relates to our main results.
We consider the problem of estimating a flexible econometric model that satisfies a set of conditional moment restrictions presented in (ref) (see also Appendix (ref)), where $z\in {\mathcal Z}\subseteq \mathbb{R}^{d}$, $X \in {\mathcal X} \subseteq \mathbb{R}^p$, $y\in \mathbb{R}$, $h\in {\mathcal H} \subseteq ({\mathcal X} \to \mathbb{R})$ for ${\mathcal H}$ a hypothesis space. For simplicity of notation we will also denote with $\psi(y; h(x))=y - h(x)$. The truth is some model $h_0$ that satisfies all the moment restrictions.
We assume we have access to a set of $n$ i.i.d. sample points $\{v_i:=(y_i,x_i, z_i)\}_{i=1}^n$ drawn from some unknown distribution $\cD$ that satisfies the moment condition in Equation (ref). We will analyze estimators that optimize an empirical analogue of the minimax objective presented in the introduction, potentially adding norm-based penalties $\Phi: {\mathcal F}\to \mathbb{R}_+$, $R: {\mathcal H}\to \mathbb{R}_+$:
where $\Psi_n(h, f) := \frac{1}{n} \sum_{i=1}^n \psi(y_i; h(x_i))\, f(z_i)$.
We assume that ${\mathcal H}$ and ${\mathcal F}$ are classes of bounded functions on their corresponding domains and, without loss of generality, their image is a subset of $[-1, 1]$. Similarly, we will also assume that $y\in [-1,1]$. The results of this section hold for a general bounded range $[-b, b]$ via standard re-scaling arguments with an extra multiplicative factor of $b$. Moreover, we will assume that ${\mathcal F}$ is a symmetric class, i.e. if $f\in {\mathcal F}$ then $-f \in {\mathcal F}$. Moreover, we will assume that ${\mathcal H}$ and ${\mathcal F}$ are equipped with norms $\|\cdot\|_{{\mathcal H}}, \|\cdot\|_{{\mathcal F}}$ and we will define the norm-constrained classes and for any function class ${\mathcal G}$ we let ${\mathcal G}_B= \{g\in {\mathcal G}: \|g\|\leq B\}$, be the $B$ bounded norm subset of the class.
Our estimation target is good generalization performance with respect to the projected residual mean squared error (RMSE), defined as the RMSE projected onto the space of instruments:
where $T:{\mathcal H} \to {\mathcal F}$ is the linear operator defined as $T h := \mathbb{E}[h(X)\mid Z=\cdot]$. This performance metric is appropriate given the ill-posedness problem well known in this setting; imposing further conditions on the strength of the correlation between the treatments and instruments (instrument strength) allows one to, translate bounds on the projected RMSE to bounds on the RMSE (see e.g. chen2012estimation and other references in the applications below).
We start by defining some preliminary notions from empirical process theory that are required to state our main results. Let ${\mathcal G}$ a class of uniformly bounded functions $g: {\mathcal V} \to [-1, 1]$ from some domain ${\mathcal V}$ to $[-1, 1]$. The localized Rademacher complexity of the function class is defined as: ${\mathcal R}_n(\delta; {\mathcal G}) = \mathbb{E}_{\{\epsilon_i\}_{i=1}^n, \{v_i\}_{i=1}^n}\left[\sup_{\substack{g\in {\mathcal G}\\ \|g\|_2\leq \delta}} \left| \frac{1}{n} \sum_{i=1}^n \epsilon_i g(v_i)\right| \right]$, where $\{v_i\}_{i=1}^n$ are i.i.d. samples from some distribution $D$ on ${\mathcal V}$ and $\{\epsilon_i\}_{i=1}^n$ are i.i.d. Rademacher random variables taking values equiprobably in $\{-1, 1\}$. We will also denote with ${\mathcal R}_n({\mathcal G})$, the un-restricted Rademacher complexity, i.e. $\delta=\infty$.
We denote with $\|\cdot\|_2$ the $\ell_2$-norm with respect to the distribution $D$, i.e. $\|g\|_2=\sqrt{\mathbb{E}_{v\sim D}[g(v)^2]}$, and analogously we define the empirical $\ell_2$-norm as $\|g\|_{2,n}=\sqrt{\frac{1}{n}\sum_i g(v_i)^2}$. In our context, where $v=(y, x, z)$, when functions take as input subsets of the vector $v$, then we will overload notation and let $\|\cdot\|_2$ and $\|\cdot\|_{2,n}$ denote the population and sample $\ell_2$ norms with respect to the marginal distribution of the corresponding input, e.g., if $h$ is a function of $x$ alone and $f$ a function of $z$ alone, we write $\|h\|_2=\sqrt{\mathbb{E}_{x}[h(x)^2]}$, $\|f\|_2=\sqrt{\mathbb{E}_{z}[f(z)^2]}$, and $\|h f\|_2=\sqrt{\mathbb{E}_{x, z}[h(x)^2\, f(z)^2]}$.
A function class ${\mathcal G}$ is said to be symmetric if $g\in {\mathcal G} \implies -g \in {\mathcal G}$. Moreover, it is said to be star-convex if: $g\in {\mathcal G} \implies r\, g \in {\mathcal G}, \forall r\in [0,1]$. The critical radius $\delta_n$ of the function class ${\mathcal G}$ is any solution to the inequality ${\mathcal R}_n(\delta; {\mathcal G}) \leq \delta^2$.
We show that, if the function space ${\mathcal F}_{U}$ contains projected differences of hypothesis spaces $h\in {\mathcal H}_B$, with some benchmark hypothesis $h_*\in {\mathcal H}_B$, i.e. $T(h-h_*)\in {\mathcal F}_U$, then a regularized minimax estimator can achieve estimation rates that are of the order of the projected root-mean-squared-error of the benchmark hypothesis $h_*$ and the critical radii of (i) the function class ${\mathcal F}_{3U}$ and (ii) a function class ${\mathcal G}$ that consists of functions of the form: $q(x)\cdot Tq(z)$, for $q=h-h_*$. The projected root mean squared error of the benchmark class can be understood as the approximation error or bias of the hypothesis space ${\mathcal H}_B$, and the critical radius can be understood as the sampling error or variance of the estimate. If $h_0\in {\mathcal H}_B$, then the approximation error is zero. We present a slightly more general statement, where we also allow for ${\mathcal F}_U$ to not exactly include $T(h-h_*)$, but rather functions that are close to it with respect to the $\ell_2$ norm. For this reason, we will need to define the following slightly more complex hypothesis space, in order to state our main theorem:
where $f_{h}^U=\operatorname*{arg\,min}_{f\in {\mathcal F}_U} \|f-T(h-h_*)\|_2$. If $T(h-h_*)\in {\mathcal F}_U$, then this simplifies to the class of functions of the form: $(h-h_*)(x)\, T(h-h_*)(z)$.
Observe that if the classes ${\mathcal H},{\mathcal F}$ already are norm constrained, then the theorem directly applies to the estimator that solely penalizes the $\ell_{2,n}$ norm of $f$, i.e.:\footnote{By setting $\lambda=\delta^2/U$, $\mu=2\lambda \left(4L^2 + 27U/B\right)$ using an $\ell_{\infty}$ norm in both function spaces and taking $U, B\to \infty$. Observe that we can also take $L=1$, since $\|Th\|_{\infty} \leq \|h\|_{\infty}$ for any $T$.}
However, as we show below, imposing norm regularization as opposed to hard norm constraints leads to adaptivity properties of the estimator.
\paragraph{Adaptivity of regularized estimator} Suppose that we know that for $B, U=1$, we have that functions in ${\mathcal H}_B, {\mathcal F}_U$ have ranges in $[-1, 1]$ as their inputs range in ${\mathcal X}$ and ${\mathcal Z}$ correspondingly. Then our Theorem requires that we set: $\lambda \geq \delta^2$ and $\mu\geq 2\lambda (4L^2 + 27)$, where $\delta^2$ depends on the critical radius of the function class ${\mathcal F}_1$ and ${\mathcal G}_1$. Observe that none of these values depend on the norm of the benchmark hypothesis $\|h_*\|_{{\mathcal H}}$, which can be arbitrary and not constrained by our theorem (see also Appendix (ref)).
For some function classes ${\mathcal H}$ that admit sparse representations, we can get an improved performance if instead of testing for classes of functions ${\mathcal F}$ that contain $T(h-h_*)$, we test functions whose linear span contains $T(h-h_*)$, i.e. that $T(h-h_*)=\sum_i w_i f_i$, assuming the weights required in this linear span have small $\ell_1$ norm. The reason being that the generalization error of linear spans with bounded $\ell_1$ norm can be prohibitively large to get fast error rates, i.e. the Rademacher complexity of the span of ${\mathcal F}$ can be much larger than ${\mathcal F}$, thereby introducing large sampling variance to our sup-loss objective. To state the improved result, we define for any function space ${\mathcal F}$: $\ensuremath{\text{span}}_{\kappa}({\mathcal F}) := \left\{\sum_{i=1}^p w_i f_i: f_i\in {\mathcal F}, \|w\|_1\leq \kappa, p\leq\infty\right\}$, i.e. the set of functions that consist of linear combinations of a finite set of elements ${\mathcal F}$, with the $\ell_1$ norm of the weights bounded by $R$. To get fast rates in this second result, we will require that the $\ell_2$-normalized $T(h-h_*)$ belongs to the span. We present the theorem in the well-specified setting, but a similar result holds in the case where $h_0\notin {\mathcal H}_B$, with the extra modification of adding a second moment penalty on $f$.
In (ref) we provide further discussion related to our main theorems: i) we provide further discussion on the adaptivity of our estimators, ii) we provide connections between the critical radius and the entropy integral and how to bound the critical radius via covering arguments, iii) we provide generic approaches to solving the optimization problem, iv) we show how to combine our main theorem on the projected MSE with bounds on the ill-posedness of the inverse problem in order to achieve MSE rates, v) we offer a discussion on the optimality of our estimation rate.
In this section we describe how (ref) applies to the case where $h_0$ lies in a Reproducing Kernel Hilbert space (RKHS) with kernel $K_{\mathcal H}:{\mathcal X}\times {\mathcal X} \to \mathbb{R}$, denoted with $\bbH_K$ and $T h_0$ lies in another RKHS with kernel $K_{{\mathcal F}}: {\mathcal Z} \times {\mathcal Z} \to \mathbb{R}$ (see (ref) for more details). We outline here the main ideas behind the three components required to apply our general theory and defer the full discussion to (ref).
First we characterize the set of test functions that are sufficient to satisfy the requirement that $T(h-h_0)\in {\mathcal F}_U$. We show (see (ref)) that if the conditional density function $p(x\mid z)$ satisfies that the function $p(x\mid \cdot)$ falls in an RKHS $\bbH_{K_{{\mathcal F}}}$, then $Th\in \bbH_{K_{{\mathcal F}}}$. Moreover, we show that under the stronger conditions (see (ref)) that $p(x\mid z)=\rho(x-z)$ and $K_{\mathcal H}(x,y) = k(x-y)$, for $k$ positive definite and continuous, then $Th\in \bbH_K$, i.e. $Th$ falls in the same RKHS as $h$. These two theorems give conrete guidance in terms of primitive assumptions, on what RKHS should be used as a test function space, so that the condition that $T(h-h_0)\in {\mathcal F}$ is satisfied.
Second, by recent results in statistical learning theory, the critical radius of any RKHS-norm constrained subset of an RKHS class with kernel $K$ and norm bound $B$, can be characterized as a function of the eigen-decay of the empirical kernel matrix ${\bf K}$ defined as ${\bf K}_{ij}=K(x_i, x_j)/n$. More concretely, it is the solution to: $B\sqrt{\frac{2}{n}}\sqrt{\sum_{j=1}^{n} \min\{\lambda_j^S, \delta^2\}} \leq \delta^2$, where $\lambda_j^S$ are the empirical eigenvalues. In the worst-case is of the order of $n^{-1/4}$. In the context of (ref), the function classes ${\mathcal F}$ and ${\mathcal G}_B$ are kernel classes, with kernels $K_{{\mathcal F}}$ and $K_{\times}((x, z), (x',z'))=K_{\mathcal H}(x,x')\cdot K_{{\mathcal F}}(z,z')$. Thus we can bound the critical radius required in the theorem as a function of the eigendecay of the corresponding empirical kernel matrices, which are data-dependent quantities.
Combining these two facts, we can then apply (ref), to get a bound on the estimation error of the minimax or regularized minimax estimator. Moreover, we show that for this set of test functions and hypothesis spaces, the empirical min-max optimization problem can be solved in closed form. In particular, the estimator in Equation (ref) takes the form:
where $K_{{\mathcal H},n} = (K_{{\mathcal H}}(x_i,x_j))_{i,j=1}^n$ and $K_{{\mathcal F},n} = (K_{{\mathcal F}}(z_i,z_j))_{i,j=1}^n$, are empirical kernel matrices, and $M = K_{{\mathcal F},n}^{1/2}({\textstyle\frac{U}{n\delta^2}}K_{{\mathcal F},n} + I)^{-1} K_{{\mathcal F},n}^{1/2}$ (where $A^\dagger$ is the Moore-Penrose pseudoinverse of $A$). Moreover, in (ref), we discuss how ideas from low rank kernel matrix approximation (such as the Nystrom method) can avoid the $O(n^3)$ running time for matrix inverse computation in the latter closed form. Finally, we show (see (ref)) that if we make further assumptions on the rate at which the operator $T$ distorts the orthonormality of the eigenfunctions of the kernel $K_{\mathcal H}$, then we can show that our estimator also implies mean-squared-error rates.
In this section we deal with high-dimensional linear function classes, i.e. the case when ${\mathcal X}, {\mathcal Z}\subseteq \mathbb{R}^p$ for $p\gg n$ and $h_0(x) = \langle \theta_0, x \rangle$ (see (ref) for more details). We will address the case when the function $\theta_0$ is assumed to be sparse, i.e. $\|\theta_0\|_{0}:=\{j\in [p]: |\theta_j|>0\}\leq s$. We will be denoting with $S$ the subset of coordinates of $\theta_0$ that are non-zero and with $S^c$ its complement. For simplicity of exposition we will also assume that $\mathbb{E}[x_i\mid z]=\langle \beta, z \rangle$, though most of the results of this section also extend to the case where $\mathbb{E}[x_i\mid z]\in {\mathcal F}_i$ for some ${\mathcal F}_i$ with small Rademacher complexity. Variants of this setting have been analyzed in the prior works of gautier2011high,fan2014endogeneity. We focus on the case where the covariance matrix $V:=\mathbb{E}[\mathbb{E}[x\mid z]\mathbb{E}[x\mid z]^\top ]$, has a restricted minimum eigenvalue of $\gamma$ and apply (ref). We note that without the minimum eigenvalue condition, our (ref) provides slow rates of the order of $n^{-1/4}$, for computationally efficient estimators that replace the hard sparsity constraint with an $\ell_1$-norm constraint.
Notably, observe that in the case of $\|\beta_0^i\|_2\leq U$, we note that if one wants to learn the true $\beta$ with respect to the $\ell_2$ norm or the functions $\mathbb{E}[x_i\mid z]$ with respect to the RMSE, then the best rate one can achieve (by standard results for statistical learning with the square loss), even when one assumes that $\sup_{z\in {\mathcal Z}} \|z\|_2\leq R$ and that $\mathbb{E}[zz^{\top}]$ has minimum eigenvalue of at least $\gamma$, is: $\min\left\{\sqrt{\frac{p}{n}}, \left(\frac{U\, R}{n}\right)^{1/4}\right\}$. For large $p\gg n$ the first rate is vacuous. Thus we see that even though we cannot accurately learn the conditional expectation functions at a $1/\sqrt{n}$ rate, we can still estimate $h_0$ at a $1/\sqrt{n}$ rate, assuming that $h_0$ is sparse. Therefore, the minimax approach offers some form of robustness to nuisance parameters, reminiscent of Neyman orthogonal methods (see e.g. Chernozhukov2018double).
In (ref) we also provide first-order iterative and computationally efficient algorithms with provable guarantees for solving the optimization problem. Moreover, we show that recent advances in online learning theory can be utilized to get fast iteration complexity, i.e. achieve error $\epsilon$ after $O(1/\epsilon)$ iterations (instead of the typical rate of $O(1/\epsilon^2)$ for non-smooth functions). Finally, in (ref), we also show if we assume that the minimum eigenvalue of $V$ is at least $\gamma$ and the maximum eigenvalue of $\Sigma=\mathbb{E}[xx^\dagger]$ is at most $\sigma$, then the same rate as the one presented in (ref) holds for the MSE, multiplied by the constant $\sqrt{\sigma/\gamma}$.
In this section we describe how one can apply the theoretical findings from the previous sections to understand how to train neural networks that solve the conditional moment problem. We will consider the case when our true function $h_0$ can be represented (or well-approximated) by a deep neural network function of $x$, for some given domain specific network architecture, and we will represent it as $h_0(x)=h_{\theta_0}(x)$, where $\theta_0$ are the weights of the neural net (see (ref) for more details). Moreover, we will assume that the linear operator $T$, satisfies that for any set of weights $\theta$, we have that $T h_{\theta}$ belongs to a set of functions that can be represented (or well-approximated) as another deep neural network architecture, and we will denote these functions as $f_w(z)$, where $w$ are the weights of the neural net.
\paragraph{Adversarial GMM Networks (AGMM)} Thus we can apply our general approach presented in (ref) (simplified for the case when $U=B=1$, $\lambda = \delta^2$, $\mu = 2\delta^2 (4L^2 + 27)$, where $L$ is a bound on the lipschitzness of the operator $T$ with respect to the two function space norms and $\delta$ is a bound on the critical radius of the function spaces ${\mathcal F}_{3}$ and $\hat{{\mathcal G}}_{1,L^2}$):
for some constant $c>1$ that depends on the lipschitzness of the operator $T$. The AGMM criterion for training neural networks is closely related to the work of bennett2019deep. However, the regularization presented in bennett2019deep is not a simple second moment penalization. Here we show that such re-weighting is not required if one simply wants fast projected MSE rates (in (ref) we provide further discussion). Moreover, in (ref), we show how to derive intuition from our RKHS analysis to develop an architecture for the test function network that under conditions is guaranteed to contain the set of functions of the form $Th$. This leads to an MMD-GAN style adversarial GMM approach, where we consider test functions of the form: $f(z) = \frac{1}{s}\sum_{i=1}^s \beta_i K(c_i, g_w(z))$, where $c_i$ are parameters that could also be trained via gradient descent. The latter essentially corresponds to adding what is known as an RBF layer at the end of the adversary neural net (denoted as KLayerTrained in experiments). Finally, in (ref), we provide heuristic methods for solving the non-convex/non-concave zero-sum game, using first order dynamics.
We will show that we can reduce the problem presented in (ref) to a regression oracle over the function space ${\mathcal F}$ and a classification oracle over the function space ${\mathcal H}$ (see (ref) for more details). We will assume that we have a regression oracle that solves the square loss problem over ${\mathcal F}$: for any set of labels and features $z_{1:n}, u_{1:n}$ it returns
Moreover, we assume that we have a classification oracle that solves the weighted binary classification problem over ${\mathcal H}$ w.r.t. the accuracy criterion: for any set of sample weights $w_{1:n}$, binary labels $v_{1:n}$ in $\{0, 1\}$ and features $x_{1:n}$:
In practice, we will consider a random forest regression method as the oracle over ${\mathcal F}$ and a binary decision tree classification method as the oracle for ${\mathcal H}$ (which we will refer to as RFIV). Prior work on random forests for causal inference has focused primarily on learning forests that capture the heterogeneity of the treatment effect of a treatment, but did not account for non-linear relationships between the treatment and the outcome variable. The method proposed in this section makes this possible. Observe that the convexity of the set $A$ is violated by the random forest function class with a bounded set of trees. Albeit in practice this non-convexity can be alleviated by growing a large set of trees on bootstrap sub-samples or using gradient boosted forests as oracles for ${\mathcal F}$. Moreover, observe that we solely addressed the optimization problem and postpone the statistical part of random forests (e.g. critical radius) to future work (see also Appendix (ref)).
In the appendix we also provide further applications of our main theorems. In (ref) we show how our theorems apply to the case where ${\mathcal H}$ and ${\mathcal F}$ are growing linear sieves, which is a typical approach to non-parametric estimation in the econometric literature (see e.g. chen2012estimation). In (ref) we analyze the case where ${\mathcal H}$ and ${\mathcal F}$ are function classes defined via shape constraints. We analyze the case of total variation bound constraints and convexity constraints. This applications provides analogues of the convex regression and the isotonic regression to the endogenous regression setting and draws connections to recent works in econometrics on estimation subject to monotonicity constraints Chetverikov2017.
\paragraph{Experimental Design.} We consider the following data generating processes: for $n_x=1$ and $n_z\geq 1$
While, when $n_x=n_z>1$, then we consider the following modified treatment equation:
We consider several functional forms for $h_0$ including absolute value, sigmoid and sin functions (more details in (ref)) and several ranges of the number of samples $n$, number of treatments $n_x$, number of instruments $n_z$ and instrument strength $\gamma$. We consider as classic benchmarks 2SLS with a polynomial features of degree $3$ (2SLS) and a regularized version of 2SLS where ElasticNetCV is used in both stages (Reg2SLS).
In addition to these regimes, we consider high-dimensional experiments with images, following the scenarios proposed in bennett2019deep where either the instrument $z$ or treatment $x$ or both are images from the MNIST dataset consisting of grayscale images of $28\times28$ pixels. We compare the performance of our approaches to that of bennett2019deep, using their code. A full description of the DGP is given in the supplementary material.
\paragraph{Results.} The main findings are: i) for small number of treatments, the RKHS method with a Nystrom approximation (NystromRKHS), outperforms all methods (Figure (ref)), ii) for moderate number of instruments and treatments, Random Forest IV (RFIV) significantly outperforms most methods, with second best being neural networks (AGMM, KLayerTrained) (Figure (ref)), iii) the estimator for sparse linear hypotheses can handle an ultra-high dimensional regime (Figure (ref)), iv) neural network methods (AGMM, KLayerTrained) outperform the state of the art in prior work bennett2019deep for tasks that involve images (Figure (ref)). The figures below present the average MSE across $100$ experiments ($10$ experiments for Figure (ref)) and two times the standard error of the average MSE.