Automatic Double Machine Learning for Continuous Treatment Effects

Automatic Double Machine Learning for Continuous Treatment Effects

\paragraph

Introduction

In this paper we propose a new nonparametric estimator of continuous treatment effects and prove asymptotic normality. Continuous treatment effects are of significant importance in applied economics. For example, labor economists study how the number of hours in a job training program impacts worker earnings flores2007estimating, and political economists seek to understand how distance to polling locations impacts propensity to vote cantoni2020precinct\footnote{Empirical examples from additional recent applied economics papers are given in Appendix (ref)}.

There are many interesting statistical objects associated with continuous treatment effects. The object we specifically focus on estimating in this paper is the expected value of the outcome variable at a given level of the treatment variable. This object is also known as an average dose-response function or the average structural function blundell2001endogeneity. Under conditional unconfoundedness assumptions, our parameter has a causal interpretation, and if conditional unconfoundedness does not hold, the parameter is still of descriptive interest.

Current applied work often estimates continuous treatment effects by imposing linear functional form assumptions. This may lead to treatment effect estimates that are biased and hard to interpret. To avoid parametric functional form assumptions, we suggest a nonparametric estimator. Our nonparametric estimator incorporates machine learning (ML) based on Lasso.

Though ML methods like Lasso produce accurate predictions, they induce bias in the parameter estimates. This “regularization” bias is due to the bias-variance trade-off inherent in minimizing the mean square error by the ML algorithms. Given the high variance of the ML methods, to keep mean square error low the algorithms must accept bias increases. This impacts not only the mean of the estimator, but estimation of its asymptotic distribution as well. Therefore simply estimating standard objects of interest by plugging in ML methods (called “plug-in” estimators) can lead to obstacles in statistical inference. For instance, using an off-the-shelf Lasso estimate of the regression function and plugging into the first step of the estimation procedure will not produce valid confidence intervals chernozhukov2016locally. There are also inherent model selection problems that arise when using ML.

This paper aims to avoid these inference problems by adapting techniques for incorporating ML methods into estimators while preserving desirable statistical properties of these estimators. Specifically we build on work in chernozhukov2016locally and chernozhukov2018double that propose solutions to “debias” machine learning methods (DML). One important step in debiasing estimators is through sample splitting: averaging over observations of the data different than those used to estimate the models. Another important debiasing step is adding a debiasing term to the equation that defines the parameter of interest.

In the initial DML literature the structure of debiasing terms is found analytically and then estimated by ML methods. However, in some settings the structure is not known, and even when it is, this estimation procedure might be undesirable. In these cases we can instead estimate these debiasing terms through a new procedure called automatic double debiased machine learning (Auto-DML) that estimates the debiasing term directly chernozhukov2018automatic. For the estimands of interest in this paper --- points on the dose-response curve --- directly estimating the debiasing term is undesirable, for reasons described below, and we will instead use Auto-DML. However, tools from the current Auto-DML framework cannot be immediately applied to our continuous treatment effect setting given the localization around treatment values that estimation requires.

The contribution of this paper is to adapt the tools of Auto-DML toward the problem of estimating continuous treatment effects. The papers that have studied the DML framework for continuous treatment effects include su2019non and colangelo2020double. This paper is most closely related to the work of colangelo2020double, which presents and proves statistical properties of the DML version of the estimator we study here. A key difference between our estimators is how the debiasing term is constructed. Different versions of debiased continuous treatment effect estimators are given by kennedy_non-parametric_2017 and kallus2018policy.

In our continuous treatment setting, the structure of the debiasing term is known, but there are several reasons we would still want to estimate continuous treatment effects with Auto-DML rather than DML. Specifically, the debiasing term includes the multiplicative inverse of the generalized propensity score (MIGPS), which is the multiplicative inverse of the probability density of the treatment variable given the covariates ($\frac{1}{f(t|X)}$)\footnote{Readers may be more familiar with the binary treatment case where the equivalent is Inverse Probability-of-Treatment Weighting (IPW) or Horvitz Thompson weights $\frac{1}{e(X)}$ where $e(x) = {\rm I}\kern-0.18em{\rm P}(T = 1| X)$. }. Current DML estimators involve plugging in an estimated MIGPS and then inverting it, which can be numerically unstable. Our Auto-DML approach estimates the MIGPS directly instead. We incorporate this estimated term additively into the bias correction term rather than inversely.

In addition, directly estimating the MIGPS has desirable balancing properties. Dealing with balancing and trimming is one of the key components of applied work, and Auto-DML is tailor-made for balancing in comparison to DML. Current DML approaches for continuous treatment effects must conduct implicit trimming by using bounded kernels to achieve numerical stability; however, with trimming, we are no longer conducting inference on the original population, but instead on a sub-population of some form. There is limited theoretical justification for ad hoc trimming and censoring of this type crump2009dealing. In comparison, Auto-DML does not require trimming through a kernel, and has been shown to work well without much additional trimming; e.g., in the binary IV treatment case singh2019biased. Empirically we find through Monte Carlo simulations that when we avoid implicit trimming our estimator decreases root mean square error (RMSE) --- up to 50% in certain specifications --- compared to current methods.

The issue of balancing and overlap is especially important in the continuous case. The number of people close to a fixed value $t$ with particular values of the covariates may be very small. Therefore, it has been noted that researchers must be careful when using inverse propensity weighting for continuous treatments because the effects may be exquisitely sensitive to the specification of conditional density hernan2010causal. A motivation for this paper is to see if the advantages of Auto-DML over DML in the binary case transfer over to the continuous case.

The paper will continue in the following way. Section (ref) sets up the framework of the paper and introduces the parameter of interest. The assumptions needed to identify this parameter of interest are given in Section (ref). In Section (ref) we describe the learning problem for estimating our parameter, which includes our novel MIGPS estimation procedure. Section (ref) provides the theoretical results. Simulation results are described in Section (ref). Section (ref) concludes.

Setup

Notation

We assume we have independent and identically distributed data $(W_1, \cdots, W_n)$ where the $W_i = (X_i, T_i, Y_i)$ are copies of a random variable $W$ with support $\mathcal{W} = \mathcal{X} \times \mathcal{T} \times \mathcal{Y}$, with a cumulative distribution function (cdf) $F_{YTX}(Y,T,X)$. We use capital letters to denote random variables and lowercase letters to denote their possible values. For each unit in a large population $X_i \in {\rm I}\kern-0.18em{\rm R}^{p}$ denotes a vector of covariates, with $p$ potentially large, and $T_i \in {\rm I}\kern-0.18em{\rm R}$ as the continuous treatment. We use the potential outcome framework rubin1974estimating and $Y_i(t)$ denotes the potential outcome that would have been observed for individual $i$ under treatment level $t$.

To simplify presentation let $\gamma_0(t, x)$ denote the true conditional expectation function ${\rm I}\kern-0.18em{\rm E}[Y| T = t, X = x]$. We also use $\alpha_0(t, X)$ to denote the true multiplicative inverse of the generalized propensity score (MIGPS) $\frac{1}{f_{T|X}(t|X)}$. As will be explained later, these two functions\footnote{they are functions and we stick to the language of the literature and call them nuisance parameters. } will be our nuisance parameters. They are called nuisance parameters because though we need them in order to estimate of parameter of interest, they themselves are not inherently of interest.

We use the $|\cdot|_q$ as the $\ell_q$ norm of a vector, and denote by $\|\cdot\|$ the $\mathcal{L}_2$ norm of a random variable i.e. $\|X_i\| = \sqrt{{\rm I}\kern-0.18em{\rm E}[X_i]^2} $ . Latter when have random matrices

We denote the kernel function by $K_h(x) = \frac{1}{h}K(\frac{X - x}{h})$\footnote{We primarily use the gaussian kernel, and in the Section (ref) we discuses the epanechnikov kernel as well}. Here $h$ is the bin-width of the kernel. Let the roughness of the kernel $K$ we are using be denoted by $R(K) := \int_{-\infty}^{\infty} K(u) du$ as in section 2.2 of hansen2009lecture.

When we use sample splitting on the data, the full sample will be split in $L$ different folds. Let $W_{\ell}$ for $\ell \in 1:L$ denote the data that is the $\ell$-th fold. Let $W_{\ell}^c$ denote the data that is the complement of the $\ell$-th fold \footnote{so if $L = 5$, then $W_1^c = W_2 \cup W_3 \cup W_4 \cup W_5 $}.

Parameter of Interest

For a fixed value $t$ the goal is to estimate

This object is also known as an average dose-response function. It is also known as the average structural function blundell2001endogeneity.

This parameter of interest $\beta_{t}$ is implicitly defined by the moment function $g$

Here $\gamma_0(t,X)$ is the true conditional expectation regression function, and it is a nuisance parameter that must be estimated in order to estimate the parameter of interest.

\paragraph{Example}

Now we define a $\beta_{t}$ for an example data generating process (DGP) in order to help the reader better understand our object of interest. Let's say we have a covariate $X_1 \sim N(1,1)$ and noise $\epsilon \sim N(0, 1)$ and treatment $T \sim N(2,1)$ and our outcome $Y = e^{T} + 3 X_1 + \epsilon$. Let's say we take $t = 0$, and so we are interested in $\beta_{0}$. In this case $\gamma_0(0, X) = e^{0} + 3X_1$ and so ${\rm I}\kern-0.18em{\rm E}[\gamma_0(0, X)] = e^{0} + 3 {\rm I}\kern-0.18em{\rm E}[X_1] = 1 + 3 = 4$. Hence in this toy example $\beta_0 = 4$

Identification

Our parameter of interest $\beta_{t}$ is a function of potential outcomes which are not observed in the data directly. Therefore, we need to make assumptions which enable us to use the observational data to do inference about ${\rm I}\kern-0.18em{\rm E}[Y(t)]$. The following characterization of $\beta_{t}$ will be useful,

The assumptions outlined in this section are the ones required for identification of causal effects. They must hold for every $t \in \mathcal{T}$ that a researchers wants to do inference on. As mentioned in kennedy_non-parametric_2017, when a researcher has randomized experimental data, this assumptions hold. With observation data these assumptions are harder to justify, and generally impossible to impossible to test, however the estimate is still useful. The estimator still gives us an adjusted measure of association that is interesting in its own right.

Conditional unconfoundedness assumption is also known as the “ignorability” assumption. It is also related to the exogeneity assumptions made in applied economics work. The assumptions intuitively means that controlling for covariates $X$, the treatment level $t$ is effectively random. See imbens_role_2000 for weaker form of this.

Overlap is also known as the “positivity” assumption. We are assuming that the propensity score is uniformly bounded away from 0 for all values in the support of the pre-treatment variables. As discussed in imbens_role_2000, when one has a continuous treatment this may be harder to satisfy than the more commonly studied binary treatment case. This condition also becomes harder to satisfy when the dimensionality of the covariates is large.

Consistency is a causal assumption that is not always explicitly stated, but is basically always assumed in some form. Consistency means that the observed outcome for individuals with treatment level $t$ equals her outcome if she had received treatment $t$ hernan2010causal.

Learning Problem

Debiased Moment

Now we are ready to explain our proposed estimation procedure. Above in equation (ref) we gave the moment function that defined our parameter of interest $g(W, \beta, \gamma_{0}) = \gamma_0(t, X) - \beta $.

If a researcher wanted to estimate the parameter with ML, a first pass at the problem could be a “plug-in” approach. In such an approach, we would use ML to fit the model for nuisance parameter $\gamma(t,X)$, and then predictions of the model would be used to create our parameter $\beta_t$ according to the moment function.

However this plug-in approach leads to large bias in our estimate of $\beta_t$ as explained in the introduction. The DML approach gives us a way to create a new debiased moment function which enables us to avoid the bias of the plug-in approach. We denote this new debiased moment function by $\psi$.

We call $\phi(w, \beta, \gamma, \alpha)$ our “debiasing term”. As with our original moment function (ref), it is also a function of $\gamma(t, X)$, but now we introduce a new second nuisance parameter $\alpha(t, X)$ corresponding to the multiplicative inverse of the propensity score (MIGPS). In the continuous treatment effect setting of our paper we know that the correct debiasing term is $\alpha_0(t, X) = \frac{1}{f(t|X)}$. Our paper will be following the Auto-DML literature and estimating $\alpha$ directly chernozhukov2018automatic. This is in contrast to the current literature, in which $\hat{f}(t|X)$ is estimated as a function of $X$ and then $\hat{\alpha}(t,X) = \frac{1}{\hat{f}(t|X)}$.

Note that our debiasing term $\phi(w, \beta, \gamma, \alpha)$ is a function of a kernel $K_h(T_i - t)$, where we are localizing around the specific treatment level we are interested in $t$. As explained in colangelo2020double, as the bin-width of the kernel $h \rightarrow 0$ we have Neyman orthogonality as defined in chernozhukov2016locally neyman1959optimal.

Estimation

We use the empirical analog of the debiased moment function (ref) as our estimator.

To construct this estimate of $\beta_t$ in practice there are two different estimation stages: stage 1 for estimating nuisance parameters $\hat{\gamma}$ and $\hat{\alpha}$, and stage 2 for estimating the parameter of interest $\hat{\beta}_t$. Now we provide a brief outline of the process because going into the details of each step

Now we go into details of the construction.

\paragraph{Stage 1}

Estimating $\hat{\alpha}(t, X)$ and $\hat{\gamma}(t, X)$.

We approximate the value of the MIGPS at $t$ by a suitable linear combination of basis functions of the observed covariates, with coefficients localized at $t$. Let $b(x)$ be a $p \times 1$ dictionary of functions.\footnote{For example, when we simulate the estimator in section (ref) we set $b(X)$ to be a fifth order polynomial set of the covariate variables.}

Our goal is to find a vector of coefficients $\hat{\rho}_t$ \footnote{Note that $\hat{\rho}_t$ has a $t$ subscript because the coefficient is for a specific level of the continuous treatment $t$} for our dictionary such that

$$\widehat{\frac{1}{f_{T|X}(t|X)} } := \hat{\alpha}(t, X) = b(X)\hat{\rho}_t.$$

As we'll discuss below, we can find the appropriate $\hat{\rho}_t$ by solving the following Lasso problem,

where $\hat{M}$ and $\hat{Q}$ are defined as

Here $r_L$ is the Lasso regularization parameter. Often, when running Lasso, $r_L$ is selected via cross-validation. However, given that we don't know what true treatment effect, we are not able to use cross validation to pick $r_L$ in a traditional way in our case. Therefore, we follow the iterative procedure in chernozhukov2018automatic to determine its value in practice.

Before continuing, we give some intuition for the structure of estimation problem.

The goal is to find $\hat{\rho}_t$ such that

If we actually observed $\frac{1}{f_{T|X}(t|X)}$ we could find the $\hat{\rho}_t$ that enables us to approximate it best by solving a weighted least squares problem of the following form,

approximate with

We formalize this intuition in the next Lemma.

Note from the first order conditions of the Lasso we get this following balancing equation

From this expression we see that $b(X_i)' \rho_{L}$ serve to approximately balancing the overall sample average with the sample average of the group with treatment value $T$ close to our $t$ of interest. This type of balancing condition for the binary treatment case is given in chernozhukov2018automatic, athey2016approximate, zubizarreta2015stable.

Finally, to estimate $\hat{\gamma}(t,X)$ (the CEF), we similarly project $\gamma_0(t,X)$ onto a $p$-dimensional dictionary. We focus on using Lasso in this paper. Given that we make predictions $\hat{\gamma}(t, X)$ at a fixed point $t$, we require uniform rates, which are available for Lasso.

\paragraph{Stage 2}

From stage 1 we have functions $\hat{\alpha}_{\ell}(t,X)$ and $\hat{\gamma}_{\ell}(t,X)$ for each of the folds $\ell$. Recall that $\hat{\alpha}_{\ell}(t,X)$ and $\hat{\gamma}(t,X)$ were fit using $W_{\ell}^C$, and now in stage two we use these functions to find fitted values for observations in $W_{\ell}$. Sum up all observations in all folds, and sum up all folds and divide by $n$ as in equation (ref)

Similarly, sum up over all observations at in equation (ref) to calculate the estimate of the variance.

Asymptotic Normality

Now we prove the asymptotic normality of the estimator we described in the previous section. Like estimation section, we start with discussing the first stage of the estimation process in Section (ref)

Stage 1

The normality of our estimator of $\beta_t$ depends on the estimation of the multiplicative inverse of the generalized propensity score (MIGPS) $\alpha(t,X)$ and the conditional expectation function $\gamma(t,X)$. We provide conditions for $\mathcal{L}_2$ convergence rates for our estimators of the nuisance parameters.

We adapt the assumptions from chernozhukov2016locally to our setting.

This assumption could also be weakened to allow for the bound on the basis functions to be an increasing function of the sample size $B_n$ rather than a constant $C$.

The next two assumptions - Assumptions (ref) and (ref) control the complexity of the true function $\alpha_0(t,X)$. Intuitively, the less complex $\alpha_0(t,X)$, the faster we can estimate it. Therefore the complexity of $\alpha_0(t,X)$ governs the convergence rate $\|\hat{\alpha}_0(t,X) - {\alpha}(t,X)\|$. When using the estimator we can assume either Assumption (ref) or (ref) - whichever we find plausible. If we use Assumption (ref) we will get a slower rate, and if we use Assumption (ref) we will get a faster rate. We call the complexity in Assumption (ref) the dense regime and in Assumption (ref) the sparse regime.

We will prove the convergence rates under the two different regimes in Lemmas (ref) and (ref) respectively.

Assumption (ref) says that $\alpha_0$ can be approximated by our dictionary $b$. An example of when this assumption holds is when $\alpha_0$ is a linear combination of the elements of the $b$ dictionary. The specific $\epsilon_n$ rate assumed here follows from Lemma (ref).

Part 3 of assumption (ref) is a population version of the restricted eigenvalue condition of bickel2009simultaneous as adapted in chernozhukov2018automatic. A clear introduction of the restricted eigenvalue condition is given in tibshirani2016closer.

To give some intuition about this part of assumption (ref), we relate it to the “no perfect multicollinearity” assumption of classic linear regression. When doing a linear regression of an outcome variable $Y$ on a covariate matrix $X$, we have to make sure that we don't have muticolinarity in order to ensure that matrix ${\rm I}\kern-0.18em{\rm E}[X^{T}X]$ is invertible. The matrix will only be invertible if none of the eigenvalues of the matrix are zero. We do not need such a strong condition in our case for matrix $Q$. Given that we are assuming a form of sparsity, we only require the invertability for the sub-matrices of $Q$ that we are considering.

To satisfy this assumption we set $a_n = \ln(\ln(n))$, following chatterjee2015prediction. In practice, we pick a data-driven following chernozhukov2018automatic.

The proof for Lemma (ref) is given in Appendix (ref). We denote the rate proven here as $\epsilon_n = \sqrt{\frac{\ln(p)}{hn}} + h^2$ for the rest of the paper. This rate impacts the rates at which the nuisance parameters can estimated. In many examples in the Auto-DML literature - like the binary treatment case with IV - the rate $\epsilon = \sqrt{\frac{\ln(p)}{n}}$. In our case we have a slower rate with an additional $h^2$ because we are estimating at a fixed point $t$.

The assumption that $\epsilon_n = o(r_L)$ means that our regularization parameter $r_L$ must go to zero slightly slower than $\epsilon_n$, which in our case $\epsilon_n = O_p(\sqrt{\frac{\ln(p)}{hn}} + h^2)$. One way to impose that $\epsilon_n = o(r_L)$ is to set $r_L$ proportional to $\epsilon_n\ln(\ln(n))$ in large samples. The $\ln(\ln(n))$ term does not impact the rates.

Hence

The proof for Lemma (ref) is given in Appendix (ref)

The proof for Lemma (ref) is given in Appendix (ref)

These regime-specific bounds on $d_{\gamma}$ are sufficient conditions for the DML product condition given in Corollary (ref). This controls the interaction remainder

Corollary (ref) is proven in Appendix (ref). In this corollary we see the trade off in the error permitted in estimating our two nuisance parameters - the regression function $\gamma_0(t, X)$ and the MIGPS $\alpha_0(t,X)$.

Stage 2

To prove the normality of our estimator we follow the structure of Lemma 15 of chernozhukov2016locally.

Proofs following the structure of Lemma 15 from chernozhukov2016locally require three different types of assumptions in order to prove asymptotic normality require three types of assumptions. First are the mild mean square consistency conditions - which are satisfied given the rate conditions given in Lemmas (ref) and (ref) Assumption (ref). Second is an assumption that the controls the interaction of the nuisance parameters, that is controlled by Corollary (ref). Lastly, there is an assumption that controls that average of the double robustness term $\psi$ n our case as $h$ goes to zero this reminder goes to zero

The $B_t$ is a bias term due to the fact the density estimators have with a non-zero asymptotic bias when one uses an optimal bandwidth $h$ - as explained in section 2.14 of hansen2009lecture. We could make $B_t$ smaller by selecting a sub-optimal bandwidth, but then our estimator would have a slower convergence rate, so we will avoid this. We derive the asymptotic bias and variance terms in Lemmas (ref) and (ref) in the appendix.

Numerical examples

Simulation study

In this section we give simulation results for our estimator, and compare our estimator to the estimator proposed by colangelo2020double. In their paper colangelo2020double give simulation results to showcase the performance of their estimator under a specific DGP, and we use the same setup for our simulations.

DGP

Consider two independent standard normal noise variables $\nu$ and $\epsilon$

Create 100 covariates $X$, and which are distributed normally also with mean zero and standard deviation 1, and are correlated with another according to the covariance matrix $\Sigma$. The $diag(\Sigma) = 1$, and $(i,j)$-entry $\Sigma_{ij} = 0.5$ for $|i-j| = 1$ and $\Sigma_{ij} = 0$ for $|i-j| > 1$ for $i,j = 1,\cdots,100$

Treatment variable $T$ is a function of noise $\nu$, the covariates $X$, and a vector of parameters $\theta$. The $j-th$ element of $\theta$ is $\theta_j = 1/j^2$. $\Phi$ denotes the CDF of $N(0,1)$.

This means that our generalized propensity score is defined as

Lastly our outcome variable,

Thus the potential outcome $Y(t) = 1.2t+1.2X'\theta+ t^2 + tX_1 + \epsilon$. The parameter of interest in the simulations is the average dose response function at $t = 0$, i.e., $\beta_0 = E[Y(0)] = 0$.

We want to point out an important detail about the DGP. Because $\Phi$ is bounded between the values of 0 and 1, our $f(t| X = x)$ can only vary between 0.22 and 0.53. Therefore we expect estimators to perform well, since there are no small propensity scores that have to be estimated precisely in order to prevent large bias.

The estimator from colangelo2020double is

In comparison to ours

A very key component to both of these estimators is the kernel $K$. We present results using both epanechnikov and gaussian kernels. We used the code from colangelo2020double papers to produce Table 1.

The first column gives the samples size of the simulation $n$. The second column gives the number of folds that were used in cross fitting $L$. Third column gives $c_h$ the size of the bin-width parameter used for kernel. The bin-width formula is given in equation (ref). In this equation $\sigma_T$ is the standard deviation of the treatment variable.

We can see that the table using the epanechnikov kernel in comparison to the gaussian kernel leads to lower bias and root mean square error (RMSE), and higher coverage. However, since the epanechnikov kernel is bounded, it is doing trimming implicitly, which as we discussed in the introduction is something we want to avoid as we could possibly lead to inference on a different population. The results of the colangelo2020double estimator with a gaussian kernel are given in Table (ref). We see that with high $c_h$ the estimator with the gaussian kernel does particularly badly - with the estimator losing proper coverage.

The simulation results for our estimator are given in Table (ref). We use the gaussian kernel in order to avoid implicit trimming. The $c_h$ is the kernel bin-width parameter like above. The $c_l$ is a parameter for the lasso regularization parameter $r_L$; the formula for $r_L$ is given by

We see that our results are pretty sensitive to the choice of $c_l$, and we shall work to reduce this sensitivity in future versions. For the reasonable $c_l$ values in most cases our estimator improves upon the Colangelo and Lee estimator with a gaussian kernel. This is especially true for the simulations with the smaller sample size of n = 500 --- there the RMSE is cut in half. With the larger sample size of n=1000 our estimators can decrease RMSE by about a third.

It is important to note the both estimators perform badly with large bin-widths --- at $c_h$ = 1.5 we see that our estimators no longer have valid coverage. Therefore the choice of the bin-width is also very important for the performance of estimators of this type.

Conclusion

In this paper, we presented a new estimator of continuous treatment effects and proved that it is asymptotically normal. Our estimator uses a new debiasing method that draws from both the DML and ADML literatures. We estimate the MIGPS – which is the critical component of the debiasing term – directly, as opposed to current methods that estimate the propensity of treatment and then invert it. Our direct estimation theoretically gives us improved numerical stability and nice automatic balancing properties. Empirically in Monte Carlo simulations, we find that our estimator decreases RMSE – up to 50% in certain specifications – compared to current methods. In future versions of our paper we plan to make three key improvements. First, on the theoretical side, we plan to explain more intuition about our estimator's balancing properties. Second, on the empirical side, we want to improve the procedure that selects the optimal regularization parameter. Lastly, we also want to apply the estimator to a real world empirical example - specifically we plant to estimate the spatial gender wage gap as was studied in liu_geography_2020.