# Nonparametric / 'deep' instrumental variables

Here’s a pretty fun paper on deep instrumental variables. Deep learning is actually totally ancillary to the core idea of the paper, which is that you can approach nonparametric instrumental variables as two supervised learning problems. Here, I’ll do two things:

- I. background:
- what nonparametric instrumental variables is
- some theory relating to identification/estimation
- how the DeepIV paper approaches the problem

- II. re-implementation:
- re-implement the core idea of breaking nonparametric IV into two supervised learning problems, but without any deep learning
- build a simple but plausible-looking simulation environment to apply this re-implementation
- discuss / illustrate some potential challenges to this core idea

### I. Background / theory

#### I.1. Instrumental variables (IV)

People are often interested in how things affect other things:

- If you sell stuff, you might care about how prices affect number of sales.
- If you run a social media company, you might be interested in how the number of ads you show affect how much people use your platform.
- If you work in education policy, you might care about how education affects peoples’ lifetime earnings.

This general situation can be described as follows:

- , where and are observable but is not
- you’re trying to learn what is from your data.

To see how each of the examples above maps to the formula above:

- is how many units something sold, is the price.
- is the number of hours someone spends on your social media site, is the fraction of their news feed that’s ads.
- is somebody’s total lifetime earnings, is their education level.

So, the most natural thing you can do to figure out how affects is to just… look at your data. As in, for a given , you can look at all the various s associated with that, and then just average them to get . It’s clear that if , then for all , so that as you get more and more data you’ll eventually learn the function .

This method of requires that . Unfortunately, this assumption is often violated. For each of the above examples:

- The higher priced items are also likely to be items that people can’t buy anywhere else, and as such may sell more, so that the ‘substitutability’ of an item drives a part of , so that is correlated with price .
- Your social media platform’s machine learning team has already optimized ad targeting a bit, so that people that get shown more ads are naturally more patient, and more patient people spend more time on your platform, so that patience drives a part of , so that is correlated with number of ads the person sees .
- People who are more intelligent tend to get more education, but intelligence itself also affects how much you earn, so that intelligence drives a part of , so that is correlated with education level .

To deal with this problem, people often turn to something called ‘instrumental variables’. An instrumental variable is something such that:

- affects
- doesn’t affect via any other channels, only through its impact on
- is independent of

In the pricing example, you might run an experiment where you randomly discount some products by some percentage. The discount amount has an obvious affect on the price , is independent of quality or anything else affecting (because it’s random), and (probably) doesn’t affect purchases except through price. So the discount percentage would be your instrument .

With your instrument, you can basically do something like this:

- look at variation in
- this variation in affects , which then affects
- so, the change in when varies must be due purely to changes , since by assumption:
- has no effect on , and
- has no other impact on

- this allows you to deduce the impact of on .

We’ll be more formal about this in the next section.

#### I.2. Nonparametric IV, more formally

(Note: you can see Newey Powell 2003 for more details. This is called ‘nonparametric’ because we’re treating the function as a general function, as opposed to being defined by a fixed set of parameters.)

Alright, let’s formally define our setup: \begin{equation} y=g(x)+\epsilon,\quad \mathbb{E}[e\mid z]=0, \quad \text{possibly } \mathbb{E}[e\mid x]\neq 0 \tag{1} \label{ivdef} \end{equation}

- are all observable, is not
- is 1-dimensional, and can be multidimensional
- is your outcome variable, are covariates, are instruments

So, the formal version of the argument given at the end of the previous section is basically integrating the expression in \eqref{ivdef}: Some points:

- drops out because .
- is observable, since it’s just the distribution of given .
- is observable, since it’s just the average of given .
- So, the problem of figuring out what is amounts to solving for given this equation.

Alright, let’s do some examples of equation \eqref{fredholm} to fix intuition:

- Suppose .
- In that case, is just the delta distribution, supported purely on .
- So, \eqref{fredholm} reduces to , which is just .
- This makes sense, since if then there’s no need for an instrument. Your is independent of . So you can just literally look at to get .

- Suppose has no impact on , as in for all .
- In that case, \eqref{fredholm} becomes which implies that .
- Thus, you’re not going to be able to learn in this case (unless is constant).
- This is also intuitive: if your instrument doesn’t have any impact on your , then it’s not possible for you to deduce the impact of on via this instrument because changing this instrument doesn’t change your at all.

- Suppose , and also that for some vector and where has rank .
- \eqref{fredholm} becomes which implies .
- You can estimate since and you observe and .
- You can estimate since you observe and and spans the same space as as you vary .
- This is just your typical linear instrumental variables setup, where first you estimate how your affects your (which is the ), and then plug in these estimates of given (which is ) to estimate .

#### I.3. Identification? Estimation?

Given that we want to back out using \eqref{fredholm}, it would be nice if there were some identification conditions telling us whether/when there’s only a single function that satisfies \eqref{fredholm} out of some class of functions . Identification holds IFF a ‘completeness’ condition holds. Roughly this completeness condition amounts to a requirement that any random variable for must somehow not exploit any information that’s present in but not in :

- Note that any function that satisfies must have
- So .
- By definition is identified exactly when implies . This last thing is the ‘completeness’ condition.

Let’s look again at the various examples above to understand this completeness condition:

- If , then clearly fixing exactly pins down , so this condtion holds, and as thus is identified.
- If has no impact on , then any that averages out to is going to satisfy , so that we will not be able to identify any that’s not constant.
- If , and also that with full rank, then may contain more information than , but the full rank condition guarantees that varying induces variation in in as many directions as itself, so that the spans the same space as . Then, because is linear, the extra information that has that’s not in is irrelevant, as variation in suffices to recover , so that is identified.

There are some theoretical results on this completeness condition, though maybe not as extensive as we’d like.

- There are fairly simple situations where completeness fails (see Severini Tripathi 2006).
- But completeness is ‘generic’ in some sense so long as you have as many instruments as regressors (see Andrews 2011).
- The sense in which completeness is generic is an extension of the standard ‘measure-1 wrt Lebesgue’ definition extended to infinite dimensional spaces.
- The genericness is relative to a fairly expansive looking class of distributions (as in, consider the set of all distributions in this fairly expansive class, then the complete ones are generic).
- This seems like a pretty useful result that suggests we shouldn’t worry too much, though I’m not entirely convinced I trust my intuition about these kind of unnatural infinite-dimensional objects.

- Also, it’s impossible to perform hypothesis testing for whether completeness holds, based on observations of (see Canay Santos Shaikh 2013).
- Though this untestability is due to the presence of incomplete distributions that are arbitrarily close to any complete distribution, and if you kind of lump these arbitrarily-close-to-being-complete distributions with the complete ones, testing becomes possible (see Freyberger 2017).

In practice, identification is good, but if we want to estimate we need some sort of consistency result. On this front, we have the problem that the solution to \eqref{fredholm} may not be a continuous function of and . This discontinuity is a problem, because typically we’ll estimate stuff by getting a bunch of data so that our empirical versions of and converge towards the truth, and then we invoke an argument by continuity that our resulting estimate of will also converge towards the truth as we get more data. If we have no guarantee of this continuity then there’s no guarantee that that our estimate of will converge to the truth even if our estimates of and become arbitrarily good. You can often circumvent this kind of issue by regularizing things appropriately, but in that case you’ll have to argue consistency on a case by case basis. For example, Newey Powell 2003 provides some consistency results assuming that is some linear combination of some basis functions + regularization.

Overall, general theory for identification/estimation of nonparametric IV is somewhat incomplete, so that we don’t have the kind of totally solid guarantees that maybe we’d like.

#### I.4. Deep IV

Thankfully, we’re now at a point in history where it’s not too insane to just throw function approximators at various problems and have it work pretty well even if it’s not entirely clear why. So, even with a lack of theoretical guarantees about identification / estimation, it’s not a bad idea to just go ahead and just try to do something reasonable, and hope that things work out.

The core idea of the Deep IV paper is to do this very reasonable thing:

- Do standard supervised learning to estimate . This is ‘stage1’.
- Use your as a proxy for the truth , find so that is approximately satisfied in your data. This is ‘stage2’

This is quite intuitive: if we have lots of data, we can get a pretty good estimate of , and also of . We’re going to hope that we live in a world where is identified and not too ill-posed. If that’s the case, then as we get a lot of data, we’ll be able to learn .

This strategy is also quite nice, because it breaks the problem of estimating into two supervised leanring problems, both of which are fairly easy to throw function approximators at.

- For the stage1, both and are observable. The complication here is that, instead of trying to predict as a function of , you’re trying to predict the entire distribution of given . Distributions are generally much more difficult to predict, particularly if is multidimensional, but there’s various things you can do. Two immediate examples:
- Assume there’s some parametric form for , so you’re only predicting a finite number of parameters. E.g. assume normal / mixture of normals (this is what the Deep IV paper does).
- If the distribution of is just 1-dimensional, then you can just estimate quantiles of .
- Discretize your -space and estimate a categorical distribution on the grid.

- For stage2, this is also a supervised learning problem: You have observations , and you’re trying to find some such that \eqref{fredholm} holds approximately in your data. That is, you want to pickg that makes and very close to each other in your data. So, you can think about minimizing the squared loss:

where we’ve replaced the internal integral with a sum, because the integral will need to be implemented by sampling a bunch of and then taking a sum.

- You can parameterize your as a function of some stuff, in which case you can just go ahead and do gradient descent or something like that. The Deep IV paper uses a neural network here.
- There’s some amount of complication here, in that there’s a sum on the inside of your square that you need to evaluate every time you want to evaluate the objective function, and that’s kind of hard to deal with. So, one thing you can do is to just move the sum outside of the square. This amounts to using a different objective function (that upper bounds the previous one via Jensen):

With this, you can just randomly draw a bunch of s for each according to , treat this new set of as your data and then just do straight-up supervised learning on this instead. However, just swapping out one objective function for another is a bit questionable, and certainly if we had consistency before we would lose it by doing this (as we’ll see below, this fails in quite trivial cases). But at the end of the day, everything that we’re doing here is a bit suspicious from a theoretical perspective, and this swapping out kind of feels not that unreasonable. So, worth trying out, especially because it’s easy.

Ok, so now regarding the ‘deep’ in DeepIV: for both stage1 and stage2, the paper proposes to use neural networks. The stage1 model uses a neural network to learn a mixture of gaussians to approximate , and stage2 uses a neural network to optimize one of the two objectives \eqref{trueobj} or \eqref{upperobj}.

### II. Re-implementation

#### II.1. Implementation details

Deep learning is great, but we can just as well use other function approximators in its stead. Especially for applications domains other than image/text, more classical methods may often be preferable for ease/performance/interpretability reasons.

So I went ahead and wrote up an alternative implementation based on tree boosting and conditional quantiles rather than neural networks and Gaussian mixtures.

The main class of interest is `NonparametricIV`

, implemented in `npiv/nonparametric_iv.py`

, which encapsulates the two stages described above.

- the stage1 model approximates via some number of predicted quantiles
- this is done through LightGBM’s quantile regression functionality.
- the default set of quantiles is defined here in
`NonparametricIV._init_stage1()`

- by default, 10 quantiles are used, from .05 to .95
- the training parameters for each quantile are all the same by default

- the training of the stage1 models is implemented here in
`NonparametricIV.train_stage1()`

- this literally just trains a LightGBM model for each quantile using the defined parameters

- these models are then used to predict the conditional quantiles for each observation in the training data in
`NonparametricIV._generate_stage2_data()`

- that is, the implementaiton here uses these predicted quantiles to approximate the conditional distribution of rather than estimating some parametric distribution and then sampling from that
- this method duplicates the data a lot, so that we have a long dataframe with a row for each (observation, quantile) pair, leading to 10x the number of rows if using default settings
- this
`stage2_data`

is then used for training the stage2 models

- the stage2 model targets either the true objective \eqref{trueobj} or the upper bound objective \eqref{upperobj} via several possible function approximators
- there are two function approximators available here:
- LightGBM, which is the default one
- and linear, which exists mostly for sanity checking purposes (or I guess if you have strong parametric beliefs about the shape of the stage2 model)
- the
`stage2_model_type`

initialization argument defines whether to use LightGBM

- the
`stage2_objective`

initialization argument defines whether to use the true objective \eqref{trueobj} or the upper bound objective \eqref{upperobj} - this gives us 4 different possible configurations of the stage 2 model, which are implemeted in
`NonparametricIV.train_stage2()`

- the two cases with upper bound objective are relatively straightforward:
- upper bound objective, linear approximator:
- this is probably the easiest case to implement, as it’s literally just linear regression

- upper bound objective, LightGBM approximator:
- this is also fairly straightforward, as it’s just straight-up regression objective using LightGBM

- upper bound objective, linear approximator:
- the two cases with the true objective are more involved, as the expression in \eqref{trueobj} has a sum inside the square
- this requires some custom objectives which we implment in
`npiv/custom_objectives.py`

)- these custom objective assume that there is a
`grp_col`

variable that indicates which rows of the`stage2_data`

correspond to the same observation (recall that we generated e.g. 10 predicted quantiles for each observation, so that we have e.g. 10 rows for each observation) `grouped_sse_loss()`

is what it sounds like: it computes the sum of squared errors, except that we sum the predicted and true y-values within all observations with a single group before taking the squared error, which is what the true objective \eqref{trueobj} requires`grouped_sse_linear()`

just uses a linear function to generated the predicted y-value, then computed the grouped sse loss`grouped_sse_loss_grad_hess()`

analytically computes the first and second derivatives of`grouped_sse_loss()`

with respect to each predicted y-value

- these custom objective assume that there is a
- true objective, linear approximator:
- this involves optimizing the custom objective
`grouped_sse_linear()`

- we do so by relying on some gradient-free optimization method from scipy (by default Nelder-Mead)

- this involves optimizing the custom objective
- true objective, LightGBM approximator
- LightGBM just requires that we provide the gradient and hessian of the objective function wrt the predicted y-values, which is why we have this
`grouped_sse_loss_grad_hess()`

function - we just call LightGBM with this custom objective, using
`grouped_sse_loss()`

as the corresponding custom metric

- LightGBM just requires that we provide the gradient and hessian of the objective function wrt the predicted y-values, which is why we have this

- this requires some custom objectives which we implment in

- the two cases with upper bound objective are relatively straightforward:

- there are two function approximators available here:

In addition, in order to assess how good our trained models are, we’ll need some functions for introspecting the behavior of these models

- for linear models, you can just… look the coefficients
- for nonparametric models, like the LightGBM models, we can do something similar, where we look at the marginal effects of each variable
- this is implemented in
`ModelWrapper.marginal_effect_plots()`

in`npiv/model_wrapper.py`

- it literally takes an input dataframe, perturbs the x-column in question a tiny bit, records the corresponding change in y, and then returns the distribution of that across all observations in the dataframe
- note that this just produces the coefficients for linear models (which forms the basis for one of the unit tests)
- as a result, if the true model is linear, we would prefer that our estimated marginal effects all be clustered very tightly around the corresponding coefficient of the true linear model

- this is implemented in

#### II.2. Synthetic experiment setting

We’ll build a synthetic setting that is (1) realistic-ish, and (2) fairly simple, to function as a minimal test for our implementation of nonparametric/deep iv. The **simulation setup is linear in both stage1 and stage2**, so as to be as trivial as possible.

Implementation is in the `IVSimulator`

class in `npiv/iv_simulator.py`

- the setting is e-commerce
- some exogenous product features
`exog_x_cols`

(intuitively, like colors or weight or something) - one endogenous feature, the
`log_price`

- endogeneity comes from the fact that product elasticity (which is unobserved) is correlated with price
- this correlation comes from the prices being chosen to be almost profit-maximizing given elasticity

- instrument is a random experimental price shock of up to is either the positive or negative direction
- y-variable
`log_sales`

- true mapping from
`exog_x_cols`

and`instrument`

to quantiles of`log_price`

is linear- the unobserved optimal prices are created here via
`generate_log_optimal_prices()`

and`generate_log_costs()`

that are both linear, and then some random noise and the instrument are added - we’ll try to recover these quantiles in stage1 of NonparametricIV

- the unobserved optimal prices are created here via
- true mapping from
`exog_x_cols`

and`log_price`

to`log_sales`

is linear, with coefficients in`log_sales_coefs`

- we’ll try to recover these true coefficients in stage2 of NonparametricIV

#### II.3. Some points, illustrated via synthetic data

We can now go about applying NonparametricIV to data generated from IVSimulator to illustrate some interesting points.

Some preliminaries:

- data:
- we’ll do a super simple example with 2 exogenous covariates
- a total of 50000 randomly generated observations for training:
- 20000 for stage1: 16000 for training, 4000 for early stopping
- 20000 for stage2: 16000 for training, 4000 for recording validation performance

- this is a lot of data for an extremely trivial setup

- preliminary analysis via linear regression and 2SLS:
- as expected, linear regression of
`log_sales`

on`exog_x_cols`

and`log_price`

does a poor job of recovering the true coefficients, due to endogeneity (cell [6]) - as expected, 2 stage least squares (e.g. linear IV) works quite well, as the simulation setting is in fact linear in both stage1 and stage2

- as expected, linear regression of

Now, some points:

**NonparametricIV works… kind of?**- section 1.a. : the stage1 models are extremely good, and predicts quantiles that are extremely concentrated around the truth
- this is certainly far better than what you would expect on real life problems

- section 1.b. : the stage2 model’s marginal effects don’t look too terrible? (we use the true objective, not the upper bound objective here, for reasons discussed below)
- the average marginal effect on the
`log_price`

is around -3, whereas the truth is -4 - for the
`exog_x_cols`

the marignal effects are around 0, whereas truth is around .7 - however, there’s fairly large variance around the average marginal effects around there mean, which is due to the trained stage2 lgb model being kind of very jagged (see cell 16 for the partial dependency plot) relative to the truth (which is linear and thus smooth)
- compared to ignoring endogeneity, this approach is cleary better, which is unsurprising because endogeneity is a pretty severe issue in this simulation
- compared to 2SLS it’s cleary worse, which is unsurprising because 2SLS makes very restrictive linearity assumptions that happen to be correct here

- the average marginal effect on the

- section 1.a. : the stage1 models are extremely good, and predicts quantiles that are extremely concentrated around the truth
**stage2 model is EXTREMELY sensitive to small imperfections in the stage1 quantiles**- in section 2, we do NonparametrIV as above, except that instead of predicting the conditional quantiles, we directly use the true ones
- the training data is the same 16000 observations

- section 2.c.IV. has the analogue of section 1 above, except with true quantiles rather than predicted
- performance is better than in section 1: marginal effect of
`log_price`

is centered at -3.6 vs -3.1 in section 1, and marigal effects on`exog_x_cols`

are also closer to the truth - this is a bit surprising, given how
*incrediby good*the stage1 quantile models in section 1 were, that going from those predictions to truth still produced large improvements in stage2 performance

- performance is better than in section 1: marginal effect of
- similarly if we take the true quantiles and add noise to them as in section 2.e., then performance predictably decreases
- this is thematically similar to a fairly key point from another economics/machine learning paper on estimating stuff in two stages, which requries L2 convergence of some stage1 estimate in order to get -consistent estimates of some stage2 parameter

- in section 2, we do NonparametrIV as above, except that instead of predicting the conditional quantiles, we directly use the true ones
**regularization is pretty hard to get right**- in both sections 1.b and 2.c.IV, the loss-minimizing number of rounds to train was very small (<1000), and the implied marginal effects computed at that point would have been much worse than the ones generated after training for 10000 rounds
- thus, performing early stopping via a validation set (as you would do in a typical supervised learning context) leads to significant bias
- however, the variance of the estimated marginal effects is also much smaller if we do early stopping
- this feels like a fairly standard example of bias-variance tradeoff

- in section 2.d. we increase a regularization parameter dramatically (
`min_gain_to_split`

), and then re-train the stage2 model- this leads to biased marginal effects, e.g. on average around -2 for
`log_price`

rather than the true value of -4 - the validation loss function is actually lower here relative to 2.c.IV

- this leads to biased marginal effects, e.g. on average around -2 for
- thus, it’s unclear if model selection via validation-set performance is a good approach here

- in both sections 1.b and 2.c.IV, the loss-minimizing number of rounds to train was very small (<1000), and the implied marginal effects computed at that point would have been much worse than the ones generated after training for 10000 rounds
**upper-bound objective is trivially inconsistent**- this is easy to see:
- in 2SLS, the stage2 just takes each observation and replaces the
`log_price`

and replaces it with`log_price`

`exog_x_cols`

,`instrument`

- in the stage2 of NonparametricIV, we’re producing multiple quantiles of (
`log_price`

`exog_x_cols`

,`instrument`

) instead of just the mean - this is identical to taking the specification in stage2 of 2SLS and adding noise to
`log_price`

`exog_x_cols`

,`instrument`

- this is errors in variables, which leads to attenuation bias

- in 2SLS, the stage2 just takes each observation and replaces the
- see sections 2.c.I for the linear version, and 2.c.III for the LGB version; both produce quite bad results

- this is easy to see:

#### II.4. Potential follow-up work

So, overall it looks like this core idea of breaking up nonparametric instrumental variables into two different supervised learning problems has some promise. Some points regarding future work:

- sharing information across quantiles
- right now, each quantile model is trained separately
- there may be some benefit to jointly estimating the conditional distribution of instead
- the DeepIV implementation does something like this, but at the cost of imposing a somewhat more structured parametric form for this distribution
- maybe some post-processing shrinkage of the predicted quantiles so that e.g. implied standard deviation for similar observations don’t vary too much

- regularization / early stopping seems hard to get right in stage2, so maybe consider using some algorithm in stage2 that doesn’t overfit so easily
- so, probably some averaging-based method like random forests, which is somewhat more robust to overfitting than gradient boosting or neural networks
- didn’t do this here because the LightGBM random forest doesn’t appear to currently support custom objectives…

- also, sample splitting is a pretty reliable way to robustify stuff, so maybe do that as well in training these models (e.g. build trees using one subset of the data, fill in the leaf values with another)
- again, something that would take some more effort to implement

- so, probably some averaging-based method like random forests, which is somewhat more robust to overfitting than gradient boosting or neural networks
- as a minimal example, we’ve only evaluated this algorithm on a fairly trivial example here, so we can probably do some more examples to build intuition
- probably try more data
- probably try more covariates
- probably try some nonlinearities

- the custom objectives used here are very slow, as they’re implemented using aggregation functions in pandas see
`npiv/custom_objectives.py`

)- this is fine for protyping, but taking 10 minutes to train a model on 16000 observations renders this fairly impractical for any real-life work
- so maybe modify LightGBM and directly implement these loss functions in C++