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.
61,765 characters · 16 sections · 59 citation commands
Causal Gradient Boosting: Boosted Instrumental Variable Regression
Gradient boosting method is considered one of the leading machine learning (ML) algorithms for supervised learning with structured data. There is a large body of evidence showing that gradient boosting dominates in a significant number of ML competitions conducted on Kaggle\footnote{For reference see \url{https://www.kaggle.com/dansbecker/xgboost}}. However, recent literature hartford2017 has shown that traditional supervised machine learning methods do not perform well in the presence of endogeneity in the explanatory variables.
A common approach to correct for the endogeneity bias is to use instrumental variables (IVs). Nonparametric instrumental variables (NPIV) methods have regained popularity among applied researchers over the last decade as they do not require imposing (possibly) implausible parametric assumptions on the target function. On the other hand, existing nonparameteric estimation techniques require the researcher to specify a target function approximation (ideally driven by some ex-ante understanding of the data generating process), e.g. a sieve space, which in turn drives the choice of unconditional moment restrictions (or simply put, the choice of IV basis functions). If the approximation is bad, it will lead to misspecification issues, and if the IVs are “weak”\footnote{By “weak” we mean that the IV basis functions are weakly correlated with the target function basis functions. It is hard to define a notion of weak instruments in the NPIV set-up of newey_powell2003, since there is no explicit reduced form. In the triangular simultaneous equations models, where the explicit reduced form exists, han2014 defined weak IVs as a sequence of reduced-form functions where the associated rank shrinks to zero.}, most likely the standard NPIV asymptotic techniques will no longer be valid. Moreover, the complexity of both modelling and estimation explodes when there are more than a handful of inputs.
In this paper, we introduce an algorithm that allows to learn the target function in the presence of endogenous explanatory variables in a data driven way, meaning that the researcher does not have to make a stance on neither the form of the target function approximation nor the choice of instruments. We build on gradient boosting algorithm to transform the standard NPIV problem into a learning problem that accounts for endogeneity in explanatory variable, and thus, we call our algorithm boostIV.
We also consider a couple extensions to the boostIV algorithm that might improve its finite sample performance. First, we show how to incorporate optimal IVs, i.e. IVs that achieve the lowest asymptotic variance chamberlain1987asymptotic. Second, we augment the boostIV algorithm with a post-processing step where we re-estimate the weights on the learnt basis functions, we call this algorithm post-boostIV. The idea is based on friedman_popescu2003 who propose to learn an ensemble of basis functions and then apply lasso to perform basis function selection.
To avoid potentially severe finite sample bias due to the double use of data, we resort to the cross-fitting idea of chernozhukov2017dml. For the boostIV algorithm we split the data to learn instruments and basis functions on different data folds. We add an additional layer of cross-fitting to the post-boostIV algorithm to update the weights on the learnt basis functions.
Our method has a number of advantages over the standard NPIV approach. First, our approach allows the researcher to be completely agnostic to the choice of basis functions and IVs. Both basis functions and instruments are learnt in a data driven way which picks up the underlying data structure. Second, the method becomes even more attractive when the dimensionality of the problem grows, as the standard NPIV methods suffer greatly from the curse of dimensionality. Intuitively, learning via boosting should be able to construct basis functions that approximately represent the underlying low dimensional data features. However, our approach does not work in purely high-dimensional settings where the number of regressors exceeds the number of observations.
We study the performance of boostIV and post-boostIV algorithms in a series of Monte Carlo experiments. We compare the performance of our algorithms to both the standard sieve NPIV estimator and a variety of modern ML estimators. Our results demonstrate that boostIV performs at worst on par with the state of the art ML estimators. Moreover, we find no empirical evidence that post-boostIV achieves superior performance compared to boostIV and vice versa. However, adding the post-processing step reduces the amount of boosting iterations needed for the algorithm to converge rendering it (potentially) computationally more efficient\footnote{To be more precise, there is a trade-off at play. One boostIV iteration takes less time than one post-boostIV iteration as the latter algorithm includes an additional estimation step plus one more layer of cross-fitting. As a result, if adding the post-processing step reduces the amount of boosting iterations significantly, then we achieve computational gains. It might not be the case otherwise.}.
This paper brings together two strands of literature. First, our approach contributes to the literature on nonparametric instrumental variables modeling. newey_powell2003 propose to replace the linear relationships in standard linear IV regression with linear projections on a series of basis functions (also see blundell2007_SNPIV for an application to Engel-curve estimation). darolles2011nonparametric and hall2005nonparametric suggest to nonparametrically estimate the conditional distribution of endogenous regressors given the instruments, $F(x|z)$, using kernel density estimators. However, despite their simplicity and flexibility, both approaches are subject to the curse of dimensionality. Machine learning literature has recently also contributed to the nonparametric IV literature. hartford2017 propose a DeepIV estimator which first estimates $F(x|z)$ with a mixture of deep generative models on which then the structural function is learned with another deep neural network. Kernel IV estimator of singh2019kiv exploits conditional mean embedding of $F(x|z)$, which is then used in the second stage kernel ridge regression. muandet2019 avoid the traditional two stage procedure by focusing on the dual problem and fitting just a single kernel ridge regression.
Second, we exploit insights from the boosting literature. Originally boosting came out as an ensemble method for classification in the computational learning theory schapire1990strength,freund1995boosting, freund1997AdaBoost. Later on friedman2000additive draw connections between boosting and statistical learning theory by viewing boosting as an approximation to additive modeling. A different perspective on boosting as a gradient descent algorithm in a function space that connects boosting to the more common optimization view of statistical inference breiman1998arcing,breiman1999prediction,friedman2001greedy. $L_{2}$-boosting introduced by buhlmann2003L2boost provides a powerful tool to learning regression functions. A comprehensive boosting review can be found in buhlmann2007boosting.
The remainder of the paper is organized as follows. Section (ref) briefly introduces the NPIV framework. Section (ref) describes the standard boosting procedure. We present boostIV and post-boostIV in Section (ref). Section (ref) talks about hyperparameter tuning. Section (ref) discusses consistency. We illustrate the numerical performance of our algorithms in Section (ref). Section (ref) concludes. All the proofs and mathematical details are left for the Appendix.
Consider the standard conditional mean model of newey_powell2003
where $y$ is a scalar random variable, $g$ is an unknown structural function of interest, $x$ is a $d_{x} \times 1$ vector of (potentially) endogenous explanatory variables, $z$ is a $d_{z} \times 1$ vector of instrumental variables, and $\varepsilon$ is an error term\footnote{The approach can easily be extended to cases where only some of the regressors are endogenous. Suppose $x = (x_{1},\,x_{2})$ where $x_{1}$ consists of endogenous regressors and $x_{2}$ is a vector of exogenous regressors. Let $w$ be a vector of excluded instruments and set $z = (w,\,x_{2})$. This perfectly fits into the model described by (ref).}. Suppose that the model is identified and the completion condition holds, i.e. for all measurable real functions $\delta$ with finite expectation,
Intuitively this condition implies that there is enough variation in the instruments to explain the variation in $x$.
The conditional expectation of (ref) yields the integral equation
where $F$ denotes the conditional cdf of $x$ given $z$. Solving for $g$ directly is an ill-posed problem as it involves inverting linear compact operators (see e.g. kress1989). Note that the model in (ref) does not have an explicit reduced form, i.e. a functional relationship between endogenous and exogenous variables, however, it is implicitly embedded in $F$. Thus, from the estimation perspective we have two objects to estimate: (i) the conditional cdf F(x|z) and (ii) the structural function $g$.
A common approach in applied work is to assume that the relationships between $y$ and $x$ as well as $x$ and $z$ are linear, which leads to a standard 2SLS estimator. However, it can be a very restrictive assumption in practice, which can result in misspecification bias. A lot of more flexible non-parametric extension to 2SLS have been developed in the econometrics literature. The standard approach is ti use the series estimator of newey_powell2003 who propose to replace the linear relationships with a linear projections on a series of basis functions.
To illustrate the approach let us approximate $g$ with a series expansion
where $p^{L}(x) = (p_{1}(x),\dots,p_{L}(x))$ is a series of basis functions. It allows us to rewrite the conditional expectation of $y$ given $z$ as
Let $q^{K}(z) = (q_{1}(z),\dots,q_{K}(z))$ be a series of IV basis functions. This implies a 2SLS type estimator of $\gamma$
where $\hat{\mathbb{E}}[p^{L}(x)|z] = q^{K}(z)\left(q^{K}(z)'q^{K}(z)\right)^{-}q^{K}(z)'p^{L}(x)$. Given $L,K \rightarrow \infty$ as $n \rightarrow \infty$, asymptotically one can recover the true structural function. However, in finite samples one has to truncate the sieve at some value. Despite that, the performance of the estimator hinges crucially on the choice of the approximating space, especially in high dimensions. Moreover, NPIV estimators suffer greatly from the curse of dimensionality which renders them inapplicable even in many applications. Alternatively, we propose a data-driven approach, which is agnostic to the choice of sieve/approximating functions.
Boosting is a greedy algorithm to learn additive basis function models of the form
where $\varphi_{m}$ are generated by a simple algorithm called a weak learner or base learner. The weak learner can be any classification or regression algorithm, such as a regression tree, a random forest, a simple single-layer neural network, etc. One could boost the performance (on the training set) of any weak learner arbitrarily high, provided the weak learner could always perform slightly better than chance\footnote{This is relevant when applied to classification problems. For regression problems any simple method such as least squares regression, regression stump, or one or two-layered neural network will work.} schapire1990strength, freund1996experiments. It is a very nice feature, since the only thing we need to make a stance on is the form of the weak learner, which is much less restrictive than choosing a sieve.
The goal of boosting is to solve the following optimization problem
where $L(y,y')$ is a loss function and $f$ is defined by (ref). Since the boosting estimator depends on the choice of the loss function, the algorithm to solve (ref) should be adjusted for a particular choice. Instead, one can use a generic version called gradient boosting friedman2001greedy, mason2000boosting, which works for an arbitrary loss function.
breiman1998arcing showed that boosting can be interpreted as a form of the gradient descent algorithm in function space. This idea then was further extended by friedman2001greedy who presented the following functional gradient descent or gradient boosting algorithm:
The key point is that we do not go back and adjust earlier parameters. The resulting basis functions learnt from the data are $\phi(x) = (\phi(x;\theta_{1}),\dots,\phi(x;\theta_{M}))$. The number of iterations $M$ is a tuning parameter, which can be optimally tuned via cross-validation or some model selection criterion (see Section (ref) for more details).
The main complication in the NPIV set-up is that $x$ is potentially endogenous, otherwise learning the structural function via boosting would be straightforward. Moreover, we cannot learn basis functions in the first step and then construct IVs in the second. Dependence of the basis functions for the structural equation on the instruments and vice versa suggests an iterative algorithm.
Before we introduce the algorithm, we need to set up the boosting IV framework first. Combining (ref) and (ref) gives
Hence, the conditional expectation of $y$ given $z$ becomes
Note that (ref) and (ref) closely resemble their standard NPIV counterparts (ref) and (ref). The only difference is that the form of basis functions for boosting must be estimated, while for the standard NPIV it has to be ex-ante specified. Unlike standard boosting, where the goal is to learn $\mathbb{E}[y|x]$, in the presence of endogeneity, we want to learn $\mathbb{E}[y|x]$, implying that in each boosting iteration we have to learn the conditional expectation of the weak learner given the IVs.
To keep things clear and simple, we focus on $L_{2}$-boosting which assumes the squared loss function. buhlmann2003L2boost show that $L_{2}$-boosting is equivalent to iterative fitting of residuals. In the IV context, it means that at step $m$ the loss has the form
where $r_{m} \equiv y - g_{m-1}$ is the current residual. Thus, at step $m$ the optimal parameters minimize the loss between the residuals and the conditional expectation of the weak learner given the instruments,
However, the conditional expectation $\mathbb{E}[\varphi(x;\theta)|z]$ is unknown and has to be estimated.
A simple way to estimate the conditional expectation in (ref) is to project\footnote{In general we do not have to use a projection, we can use a more complex model to estimate the conditional expectation.} the weak learner on the space spanned by IVs
where $P_{A} = A(A'A)^{-1}A'$ is a projection matrix. The exogeneity condition in (ref) implies that any function of $z$ can serve as an instrument. However, we do not need any function, we need such a transformation of $z$ that will give us strong instruments, i.e. instruments that explain the majority of the variation in the endogenous variables. We follow gandhi2019cross_learner and introduce an additional step on which we learn the instruments. Let $\mathcal{H}(\cdot;\eta)$ be a class of IV functions parameterized by $\eta$. This formulation allows us to use various off-the-shelf algorithms such as Neural Networks, Random Forests, etc. to learn $\mathcal{H}(\cdot;\eta)$. Given the learnt IV transformation $\mathcal{H}(\cdot;\eta)$, we can rewrite (ref) as
Since the basis function parameters $(\alpha,\,\theta)$ depend on the IV transformation parameters $\eta$ and vice verse, we propose an algorithm that iterates between two steps. At the first step we learn instruments, i.e. $\eta_{m}$, given the basis functions parameter estimates from the previous iteration $(\alpha_{m-1},\,\theta_{m-1})$, then at the second step we learn new parameter estimates $(\alpha_{m},\,\theta_{m})$ given the instruments from the first step. We can draw an analogy with the canonical two-stage least squares, where we estimate the reduced form in the first stage, and the structural equation in the second. The details are provided in Algorithm (ref).
We call this algorithm the Naive boostIV, since we use the same data to learn both the instruments and the basis functions. Asymptotically this will not affect the properties of the estimator, however, in finite samples biases from the first stage will propagate to the second. This issue can be especially severe if we use regularized estimators in the first stage as the regularization bias will heavily affect the second stage estimates. To get around this issue we resort to cross-fitting.
Let $D = \{y_{i},x_{i},z_{i}\}_{i = 1}^{n}$ be our data set, where $D_{i}$ are iid. Split the data set into a $K$-fold partition, such that each partition $D_{k}$ has size $\left\lfloor \frac{n}{K} \right\rfloor$, and let $D_{k}^{c}$ be the excluded data. The boostIV procedure with cross-fitting is described in Algorithm (ref).
Our boostIV algorithm also allows to incorporate optimal instruments in the sense of chamberlain1987asymptotic, i.e. instruments that achieve the smallest asymptotic variance. Assuming conditional homoskedasticity, the optimal instrument vector of chamberlain1987asymptotic at step $m$ is
where
is the conditional expectation of the derivative of the conditional moment restriction with respect to the boosting parameters, and $\sigma^{2}_{m} = \mathbb{E}[r_{m}^{2}|z]$ is the conditional variance of the error term at step $m$. Thus, the IV transformation parameters $\eta_{m}$ are implicitly embedded in a particular approximation used to estimate $D_{m}(z)$.
The main complication with using optimal IVs is that they are generally unknown, hence, the common approach is to consider approximations. The parametrization in (ref)-(ref) allows us to use any off-the-shelf statistical/ML method to estimate the optimal functional form for the instruments. Moreover, the iterative nature of the algorithm allows us to use the estimates from step $m-1$ as proxies.
An important feature of the forward stage-wise additive modeling is that we do not go back and adjust earlier parameters. However, we might want to revisit the weights on the learnt basis functions to achieve a better fit. This can be seen as a way of post-processing our boostIV procedure.
The whole procedure can be broken down into two stages:
Note that the basis functions $(\hat{\varphi}_{1}(x),\,\dots,\,\hat{\varphi}_{M}(x))$ are causal in the sense that they are constructed using estimated parameters $\theta$ that identify a causal relationship between $x$ and $y$.
Boosting is an example of an ensemble method which combines various predictions with appropriate weights to get a better prediction. In the context of boostIV it works in the following way. We exploit the variation in IVs to get causal parameters $\theta$. Given the estimated parameters, we can treat each learnt basis function $\varphi(x;\theta_{m})$ as a separate prediction obtained by fitting a base learner. Then, the post-processing step in (ref) can be simply seen as model averaging.
We can estimate optimal weights $\hat{\beta}$ by simply running a least squares regression as in (ref) or use any other method such as Random Forrests, Neural Networks, boosting, etc. To avoid carrying over any biases from the estimation of $(\hat{\varphi}_{1}(x),\,\dots,\,\hat{\varphi}_{M}(x))$ into the choice of $\beta$, we use cross-fitting once again, which is a generalization of the stacking idea of wolpert1992.
Boosting performance crucially depends on the number of boosting iterations, in other words, $M$ is a tuning parameter. A common way to tune any ML algorithm is cross-validation (CV). The most popular type of CV is $k$-fold CV. The idea behind $k$-fold CV is to create a number of partitions (validation datasets) from the training dataset and fit the model to the training dataset (sans the validation data). The model is then evaluated against each validation dataset and the results are averaged to obtain the cross-validation error. In application to boosting, we can estimate the CV error for a grid of candidate tuning parameters (number of iterations) and pick $M^{*}$ that minimizes the CV error. Alternatively, buhlmann2007boosting show how to apply AIC and BIC criteria to boosting in the exogenous case. However, it is not clear how to adjust those criteria for the presence of endogeneity.
Both standard $k$-fold cross validation and the model selection criteria considered in buhlmann2007boosting can be computationally costly as it is necessary to compute all boosting iterations under consideration for the training data. To surpass this issue, we apply early stopping to $k$-fold CV. The idea behind early stopping is to monitor the behavior of the CV error and stop as soon as the performance starts decreasing, i.e. CV error goes up.
Algorithm (ref) provides implementation details for the $k$-fold CV with early stopping for either boostIV or post-boostIV procedure. The early stopping criterion compares the CV error evaluated for the model based on $M_j$ boosting iterations to the CV error evaluated for the model based on $M_i$, $M_i < M_j$. If $CV^{err}(M_j) > CV^{err}(M_i) + \epsilon$, where $\epsilon > 0$ but close to zero is a numerical error tolerance level, then we stop and set $M^{*} = M_i$, otherwise, continue the search. If the criterion is not met for any of the candidate tuning parameters, we pick the largest value $M^{*} = \bar{M}$.
An alternative solution would be to use a slice of the dataset as the validation sample and tune the number of iterations using the observations from the validation sample. We actually use this approach in our simulations since it significantly reduces the computational burden.
In this section, we show that under mild conditions boostIV is consistent. Theoretical properties of post-boostIV are beyond the scope of the paper and are left for future research.
We borrow the main idea from zhang_yu2005boostingES and modify it accordingly to apply it to the GMM criterion. Let $g(W_{i}, f) = (y_{i} - f(x_i))z_{i}$ denote a $k \times 1$ moment function, then $g_{0}(f) = \mathbb{E}[g(W_{i}, f)]$ is the population moment function and $\hat{g}(f) = n^{-1}\sum_{i=1}^{n}g(W_{i}, f)$ is its sample analog. Also let $\Omega$ denote a $k \times k$ positive semi-definite weight matrix and $\hat{\Omega}$ be its sample analog. Thus, the population GMM criterion and its sample analog are
The form of the GMM criterion in (ref) corresponds to the form of the empirical objective function in zhang_yu2005boostingES with the loss function replaced by the moment function.
We follow zhang_yu2005boostingES and replace the functional gradient decent step (ref) leading to the 2SLS fitting procedure on every iteration with an approximate minimization involving a GMM criterion. We can do that since the 2SLS solution is a special case of a GMM solution with an appropriate weighting matrix.
As zhang_yu2005boostingES show, the consistency of the boosting procedure consists of two parts: (i) numerical convergence of the procedure itself, i.e. the algorithm achieves the true minimum of the objective function, and (ii) statistical convergence that insures the uniform convergence of the sample criterion to its population analog. We will treat these two steps separately in the following subsections, and then combine them to demonstrate consistency of the boostIV.
To demonstrate numerical convergence, we first have to verify that the sample GMM criterion in (ref) satisfies Assumption 3.1 from zhang_yu2005boostingES.
Following zhang_yu2005boostingES, I introduce some additional notation. Let $\mathcal{S}$ be a set of real-valued functions and define
which forms a linear function space. Also, for all $f \in \text{span}(S)$ define the 1-norm with respect to the basis $S$ as
Note that Assumption (ref)(a) restricts the step size $\alpha_{m}$. friedman2001greedy argues that restricting the step size is always preferable in practice, thus, we will restrict our attention to this case\footnote{zhang_yu2005boostingES provide a short discussion on how to deal with the unrestricted step size, however, the argument relies on the exact minimization which greatly complicates the analysis.}. Moreover, $\Lambda_m$ is allowed to depend on the previous steps of the algorithm. Assumption (ref)(b) requires the step size $h_j$ to be small ($\sum_{j=0}^{\infty}h_{j}^{2}<\infty$) preventing large oscillation, but not too small ($\sum_{j=0}^{\infty}h_{j} = \infty$) ensuring that $f_{m}$ can cover the whole $\text{span}(\mathcal{S})$. The following theorem establishes the main numerical convergence result.
We need to show that the sample GMM criterion uniformly converges to its population analog, then under proper regularity conditions we will be able to ensure consistency of boostIV.
To show that the sample GMM criterion converges uniformly to its population analog, we will first bound the moment function and then we will show that it is sufficient to put a bound on the criterion function.
To bound the rate of uniform convergence of the moment function, we appeal to the concept of Rademacher complexity. Let $\mathcal{H} = {h(w)}$ be a set of real-valued functions. Let $\{\zeta_{i}\}_{i=1}^{n}$ be a sequence of binary random variables such that $\zeta_{i}$ takes values in $\{-1,1\}$ with equal probabilities. Then the sample or empirical Rademacher complexity of class $\mathcal{H}$ is given by
We also denote $R(\mathcal{H}) = \mathbb{E}_{W}\hat{R}(\mathcal{H})$ to be the expected Rademacher complexity, where $\mathbb{E}_{W}$ is the expectation with respect to the sample $W = (W_{1},\,\dots,\,W_{n})$. Note that the definition in (ref) differs from the standard definition of Rademacher complexity where there is an absolute value under the supremum sign vaart_wellner1996weak. The current version of Rademacher complexity has the merit that it vanishes for function classes consisting of single constant function, and is always dominated by the standard Rademacher complexity. Both definitions agree for function classes which are closed under negation meir_zhang2003generalization.
For many classes the Rademacher complexity can be calculated directly, however, to obtain a more general result we need to bound $R(\mathcal{S})$. Using the results from Section 4.3 in zhang_yu2005boostingES we can bound the expected Rademacher complexity of the weak learner class by
where $C_{\mathcal{S}}$ is a constant that solely depends on $\mathcal{S}$. zhang_yu2005boostingES also show that popular weak learners such as two-level neural networks and trees basis functions satisfy the requirements. However, zhang_yu2005boostingES point out that in general the bound may be slower than root-n. In Appendix (ref) we derive an alternative bound on $R(\mathcal{S})$ that works for any class with finite VC dimension. The derived VC bound is slower by the factor of $\log(n)$ that appears in a lot of ML algorithms.
Condition (ref) allows us to bound the moment function which leads to a bound on the rate of uniform convergence of the GMM criterion. The formal statements of the results are presented below.
In this section we put together the arguments for numerical and statistical convergence presented in the previous subsections to prove consistency of the boostIV algorithm. We start with a general decomposition illustrating the proof strategy and highlighting where exactly numerical and statistical convergence step in.
Suppose that we run the boostIV algorithm and stop at an early stopping point $\hat{m}$ that satisfies $\mathbb{P}(||\hat{f}_{\hat{m}}||_{1} \leq \beta_{n}) = 1$ for some sample-independent $\beta_{n} \geq 0$. Let $f^{*}$ be a unique minimizer of the population criterion, i.e. $Q(f^{*}) = \inf_{f \in \text{span}(\mathcal{S})}Q(f)$. By the triangle inequality, we get the following decomposition
We can bound the first term using the uniform bound on the sample GMM criterion in Theorem (ref), this is the statistical convergence argument. In order to bound the second term, we have to appeal to the numerical convergence argument in Theorem (ref). As a result, since $Q(\hat{f}_{\hat{m}}) \rightarrow Q(f^{*})$ as $n\rightarrow\infty$, it follows that $\hat{f}_{\hat{m}} \overset{p}{\to} f^{*}$. The following theorem formalizes the result.
To begin with, we consider a simple low-dimensional scenario with one endogenous variable and two instruments.
where instruments $z_{j} \sim U[-3,\,3]$ for $j = 1,\,2$, $e \sim \mathcal{N}(0,\,1)$ is the confounder, $\delta,\,\gamma \sim \mathcal{N}(0,\,0.1)$ are additional noise components, and $\rho$ is the parameter measuring the degree of endogeneity, which we set to 0.5 in the simulations. We focus on four specifications of the structural function:
We compare the performance of boostIV and post-boostIV with the standard NPIV estimator using the cubic polynomial basis, Kernel IV (KIV) regression of singh2019kiv\footnote{Code: \url{https://github.com/r4hu1-5in9h/KIV}}, DeepIV estimator of hartford2017\footnote{We use the latest implementation of the econML package: \url{https://github.com/microsoft/EconML}} and DeepGMM estimator of bennett2019deepGMM\footnote{Code: \url{https://github.com/CausalML/DeepGMM}}. We use $1,000$ observations for both train and test sets and $500$ observations for the validation set. Our results are based on $200$ simulations for each scenario.
We plot our results in Figures (ref) which shows the average out of sample fit across simulations (orange line) compared to the true target function (black line). Table (ref) presents the out-of-sample MSE across simulations. First thing to notice is that NPIV fails to capture different functional form subtleties. Second, DeepIV's performance does not improve upon the one of NPIV. Moreover, even though DeepGMM estimates have lower bias than the ones of NPIV and DeepIV (except for the log function), they are quite volatile across simulations leading to higher MSE. BoostIV performs on par with KIV both in terms of the bias term as they are able to recover the underlying structural relation, and in terms of the variance leading to low MSE. Finally, the post-processing step helps to further improve upon boostIV's performance by reducing bias. On top of that, post-boostIV requires less iterations to converge. We use $5,000$ iterations for boostIV, while post-boostIV uses on average $50$ iterations.\footnote{In this experiment we do not tune boostIV, we just pick a large enough number of iterations for it to converge. However, we do tune post-boostIV.}
Consider the following data generating process:
where $y_{i} \in \mathbb{R}$ is the response variable, $x_{i} \in \mathbb{R}^{d_{x}}$ is the vector of potentially endogenous variables, $z_{i} \in \mathbb{R}^{d_{z}}$ is the vector of instruments, $\varepsilon_{i} \in \mathbb{R}$ is the structural error term, and $v_{i} \in \mathbb{R}^{d_{x}}$ is the vector of the reduced form errors. Function $h(\cdot)$ is the structural function of interest, and function $g(\cdot)$ governs the reduced form relationship between the endogenous regressors and instrumental variables.
Instruments are drawn from a multivariate normal distribution, $z_{i} \sim \mathcal{N}(0, \Sigma_{z})$, where $\Sigma_{z}$ is just an identity matrix. The error terms are described by the following relationship:
where $\rho$ is the correlation between $\varepsilon$ and the elements of $v$, which controls the degree of endogeneity.
We consider two structural function specifications:
We also consider two different choices of the reduced form function $g(\cdot)$:
We use $1,000$ observations for the train set and $500$ observations for both the validation and test sets. We run $200$ simulations for each scenario. The results are summarized in Tables (ref) and (ref).
In this section we apply our algorithms to a more economically driven example of demand estimation. Demand estimation is a cornerstone of modern industrial organization and marketing research. Besides its practical importance, it poses a challenging estimation problem which modern econometric and statistical tools can be applied to.
We consider a nonparametric demand estimation framework of gandhi2020flexible (hereafter, GNT). GNT is a flexible framework that combines the nonparametric identification arguments of berry2014identification with the dimensionality reduction techniques of gandhi_houde2019.
In market $t$, $t=1,\dots,T$, there is a continuum of consumers choosing from a set of products $\mathcal{J} = \{0,1,\dots,J\}$ which includes the outside option, e.g. not buying any good. The choice set in market $t$ is characterized by a set of product characteristics $\chi_{t}$ partitioned as follows:
where $x_{t} \equiv (x_{1t}, \dots, x_{Jt})$ is a vector of exogenous observable characteristics (e.g. exogenous product characteristics or market-level income), $p_{t} \equiv (p_{1t}, \dots, p_{Jt})$ are observable endogenous characteristics (typically, market prices) and $\xi_{t} \equiv (\xi_{1t}, \dots, \xi_{Jt})$ represent unobservables potentially correlated with $p_{t}$ (e.g. unobserved product quality). Let $\mathcal{X}$ denote the support of $\chi_{t}$. Then the structural demand system is given by
where $\Delta^{J}$ is a unit $J$-simplex. The function $\sigma$ gives, for every market $t$, the vector $s_{t}$ of shares for the $J$ goods.
Following berry2014identification, we partition the exogenous characteristics as $x_{t} = \left(x_{t}^{(1)}, x_{t}^{(2)} \right)$, where $x_{t}^{(1)} \equiv \left(x_{1t}^{(1)}, \dots, x_{Jt}^{(1)}\right)$, $x_{jt} \in \mathbb{R}$ for $j \in \mathcal{J} \backslash \{0\}$, and define the linear indices
and let $\delta_{t} \equiv (\delta_{1t}, \dots, \delta_{Jt})$. Without loss of generality, we can normalize $\beta_{j} = 1$ for all $j$ (see berry2014identification for more details). Given the definition of the demand system, for every market $t$
Following berry2013connected and berry2014identification, we can show that there exists at most one vector $\delta_{t}$ such that $s_{t} = \sigma\left(\delta_{t}, p_{t}, x_{t}^{(2)}\right)$, meaning that we can write
We can rewrite (ref) in a more convenient form to get the following estimation equation
Note that in (ref) the inverse demand is indexed by $j$, meaning that we have to estimate $J$ inverse demand functions. To circumvent this problem, gandhi_houde2019 suggest transforming the input vector space under the linear utility specification to get rid of the $j$ subscript. GNT follow this idea and show that Equation (ref) can be rewritten as
where $g$ is such that
and $\omega_{jt} \equiv \left(s_{jt}, \{s_{kt}, d_{jkt}\}_{j \neq k}\right)$, where $d_{jkt} = \tilde{x}_{jt} - \tilde{x}_{kt}$ and $\tilde{x}_{t} \equiv \left(p_t, x_{t}^{(2)}\right)$.
Let $y_{jt} \equiv \log(s_{jt} / s_{0t}) - x_{jt}^{(1)}$, then we can rewrite equation (ref) in a more convenient form
Thus, (ref) is our structural equation where $\omega_{jt}$ contains endogenous variables. Assume we have a cost-shifter $w_{jt}$ that is exogenous, then given $\mathbb{E}[\xi_{jt}|x_{t}, w_{t}] = 0$ for $j \in \mathcal{J}\backslash\{0\}$\footnote{We also need the completeness condition to be satisfied, see berry2014identification for more details.} we can estimate the inverse demand function $g$. To constructs instruments, GNT transform the input space $(x_{t},\,w_{t})$ similarly to $\omega_{jt}$. Let $\zeta_{jt} \equiv \{\Delta_{jkt}\}_{j\not=k}$, where $\Delta_{jkt} = z_{jt} - z_{kt}$ and $z_{t} = (x_{t},\,w_{t})$. Thus, we can perform estimation based on $\mathbb{E}[\xi_{jt}|\zeta_{jt}] = 0$.
In our model design there are $T=100$ markets with $J \in \{10,\,20\}$ products and $K \in \{10,\, 20\}$ nonlinear characteristics besides the price. We compare the performance of boostIV and post-boostIV with KIV and DeepGMM. We drop NPIV since in our design it fails due to the curse of dimensionality. We also drop DeepIV as it suffers from the exploding gradient problem.
Table (ref) summarizes the results. The first thing to notice is that KIV performs the worst, while in the previous experiments it was one of the best performing estimators. It has both high bias and variance. DeepGMM has smaller bias, but the variance is still big. Our algorithms clearly dominate KIV and GMM, with post-boostIV delivering the best MSE results, while having slightly higher bias compared to boostIV.
In this paper we have introduced a new boosting algorithm called boostIV that allows to learn the target function in the presence of endogenous regressors. The algorithm is very intuitive as it resembles an iterative version of the standard 2SLS regression. We also study several extensions including the use of optimal instruments and the post-processing step.
We show that boostIV is consistent and demonstrates an outstanding finite sample performance in the series of Monte Carlo experiments. It performs especially well in the nonparameteric demand estimation example which is characterized by a complex nonlinear relationship between the target function and explanatory features.
Despite all the advantages of boostIV, the algorithm does not allow for high-dimensional settings where the number of regressors and/or instruments exceeds the number of iterations. We also believe it is possible to extend our algorithm in the spirit similar to XGBoost chen2016xgboost that could decrease the computation time taken by the algorithm. These would be interesting directions for future research.