ddml: Double/debiased machine learning in Stata

ddml: Double/debiased machine learning in Stata

\graphicspath{{fig/}}

\makeatletter \makeatother

\hypersetup{hidelinks}

\inserttype[st0001]{article} { Achim Ahrens\\ETH Z\"urich \\[email removed] \and Christian B. Hansen\\ University of Chicago\\[email removed] \and Mark E. Schaffer\\ Heriot-Watt University\\Edinburgh, United Kingdom \\ [email removed] \and Thomas Wiemann\\ University of Chicago\\ [email removed] }

Introduction

Identification of causal effects frequently relies on an unconfoundedness assumption, requiring that treatment or instrument assignment is sufficiently random given observed control covariates. Estimation of causal effects in these settings then involves conditioning on the controls. Unfortunately, estimators of causal effects that are insufficiently flexible to capture the effect of confounds generally do not produce consistent estimates of causal effects even when unconfoundedness holds. For example, blandhol2022tsls highlight that TSLS estimands obtained after controlling linearly for confounds do not generally correspond to weakly causal effects even when instruments are valid conditional on controls. Even in the ideal scenario where theory provides a small number of relevant controls, theory rarely specifies the exact nature of confounding. Thus, applied empirical researchers wishing to exploit unconfoundedness assumptions to learn causal effects face a nonparametric estimation problem.

Traditional nonparametric estimators suffer greatly under the curse of dimensionality and are quickly impractical in the frequently encountered setting with multiple observed covariates.\footnote{For example, the number of coefficients in polynomial series regression with interaction terms increases exponentially in the number of covariates.} These difficulties leave traditional nonparametric estimators essentially inapplicable in the presence of increasingly large and complex data sets, e.g.\ textual confounders as in roberts2020a or digital trace data hangartner2021a. Tools from supervised machine learning have been put forward as alternative estimators. These approaches are often more robust to the curse of dimensionality via the exploitation of regularization assumptions. A prominent example of a machine learning-based causal effects estimator is Post-Double Selection Lasso (PDS-Lasso) of Belloni2014a, which fits auxiliary lasso regressions of the outcome and treatment(s), respectively, against a menu of transformed controls. Under an approximate sparsity assumption, which posits that the DGP can be approximated well by a relatively small number of terms included in the menu, this approach allows for precise treatment effect estimation. The lasso can also be used for approximating optimal instruments Belloni2012. Lasso-based approaches for estimation of causal effects have become a popular strategy in applied econometrics gilchrist2016something,dhar2022reshaping, partially facilitated by the availability of software programs in Stata (\stcmd{pdslasso}\xspace, Ahrens2018; statacorp2019) and R (hdm, Chernozhukov2016b).

Although approximate sparsity is a weaker regularization assumption than assuming a linear functional form that depends on a known low-dimensional set of variables, it may not be suitable in a wide range of applications. For example, giannone2021economic argue that approximate sparsity may provide a poor description in several economic examples. There is thus a potential benefit to expanding the set of regularization assumptions and correspondingly considering a larger set of machine learners including, for example, random forests, gradient boosting, and neural networks. While the theoretical properties of these estimators are an active research topic (see, e.g., Athey2019a, farrell2021deep), machine learning methods are widely adopted in industry and practice for their empirical performance. To facilitate their application for causal inference in common econometric models, Chernozhukov2018 propose Double/Debiased Machine Learning (DDML), which exploits Neyman orthogonality of estimating equations and cross-fitting to formally establish asymptotic normality of estimators of causal parameters under relatively mild convergence rate conditions on nonparametric estimators.

DDML increases the set of machine learners that researchers can leverage for estimation of causal effects. Deciding which learner is most suitable for a particular application is difficult, however, since researchers are rarely certain about the structure of the underlying data generating process. A practical solution is to construct combinations of a diverse set of machine learners using stacking Wolpert1992,Breiman1996a. Stacking is a meta learner given by a weighted sum of individual machine learners (the “base learners”). When the weights corresponding to the base learners are chosen to maximize out-of-sample predictive accuracy, this approach hedges against the risk of relying on any particular poorly suited or ill-tuned machine learner.

In this article, we introduce the Stata package \stcmd{ddml}\xspace, which implements DDML for Stata.\footnote{This article refers to version v1.4.2 of \stcmd{ddml}\xspace.} \stcmd{ddml}\xspace adds to a small number of programs for causal machine learning in Stata (Ahrens2018, statacorp2019, statacorp2021). We briefly summarize the four main features of the program:

The article proceeds as follows. Section (ref) outlines DDML for the Partially Linear and Interactive Models under conditional unconfoundedness assumptions. Section (ref) outlines DDML for Instrumental Variables (IV) models. Section (ref) discusses how stacking can be combined with DDML and provides evidence from Monte Carlo simulations illustrating the advantages of DDML with stacking. Section (ref) explains the features, syntax and options of the program. Section (ref) demonstrates the program's usage with two applications.

DDML with Conditional Unconfoundedness

This section discusses DDML for the Partially Linear Model and the Interactive Model in turn. The exposition follows Chernozhukov2018. Both models are special cases of the general causal model

where $f_0$ is a structural function, $Y$ is the outcome, $D$ is the variable of interest, $\bm{X}$ are observed covariates, and $U$ are all unobserved determinants of $Y$ (i.e., other than $D$ and $\bm{X}$).\footnote{Since in (ref), $(D, \bm{X}, U)$ jointly determine $Y$, the model is also dubbed the “all causes model” (see, e.g., heckman2007econometric). Note that the model can equivalently be put into potential outcome notation with potential outcomes defined as $Y(d)\equiv f_0(d, \bm{X}, U)$.} The key difference between the Partially Linear Model and the Interactive Model is their position in the trade-off between functional form restrictions on $f_0$ and restrictions on the joint distribution of observables $(D,\bm{X})$ and unobservables $U$. For both models, we highlight key parameters of interest, state sufficient identifying assumptions, and outline the corresponding DDML estimator. A random sample $\{(Y_i, D_i, \bm{X}_i)\}_{i=1}^n$ from $(Y, D, \bm{X})$ is considered throughout.

The Partially Linear Model (partial)

The Partially Linear Model imposes the estimation model

where $\theta_0$ is a fixed unknown parameter. The key feature of the model is that the controls $\bm{X}$ enter through the unknown and potentially nonlinear function $g_0$. Note that $D$ is not restricted to be binary and may be discrete, continuous or mixed. For simplicity, we assume that $D$ is a scalar, although \stcmd{ddml}\xspace allows for multiple treatment variables in the Partially Linear Model.

The parameter of interest is $\theta_0$, the causal effect of $D$ on $Y$.\footnote{The interpretation of $\theta_0$ can be generalized. For example, the results of angrist1999empirical imply that in the general causal model (ref), $\theta_0$ is a positively weighted average of causal effects (e.g., conditional average treatment effects) under stronger identifying assumptions. The basic structure can also be used to obtain valid inference on objects of interest, such as projection coefficients, in the presence of high-dimensional data or nonparametric estimation without requiring a causal interpretation.} The key identifying assumption is given in Assumption (ref).\footnote{Discussions of Partially Linear Model typically show identification under the stronger assumption that $E[U\vert D,\bm{X}]=0$. We differentiate here to highlight differences between the Partially Linear Model and Interactive Model.}

To show identification of $\theta_0$, consider the score

where $\bm{W} \equiv (Y,D,\bm{X})$, and $\ell$ and $m$ are nuisance functions. Letting $m_0(\bm{X})\equiv E[D\vert \bm{X}] $ and $\ell_0(\bm{X})\equiv E[Y\vert \bm{X}]$, note that

by Assumption (ref). When in addition $E[Var(D\vert \bm{X})] \neq 0$, we get

Equation (ref) is constructive in that it motivates estimation of $\theta_0$ via a simple two-step procedure: First, estimate the conditional expectation of $Y$ given $\bm{X}$ (i.e., $\ell_0$) and of $D$ given $\bm{X}$ (i.e., $m_0$) using appropriate nonparametric estimators (e.g., machine learners). Second, residualize $Y$ and $D$ by subtracting their respective conditional expectation function (CEF) estimates, and regress the resulting CEF residuals of $Y$ on the CEF residuals of $D$. This approach is fruitful when the estimation error of the first step does not propagate excessively to the second step. DDML leverages two key ingredients to control the impact of the first step estimation error on the second step estimate: 1) second step estimation based on Neyman orthogonal scores and 2) cross-fitting. As shown in Chernozhukov2018, this combination facilitates the use of any nonparametric estimator that converges sufficiently quickly in the first and potentially opens the door for the use of many machine learners.

Neyman orthogonality refers to a property of score functions $\psi$ that ensures local robustness to estimation errors in the first step. Formally, it requires that the Gateaux derivative with respect to the nuisance functions evaluated at the true values is mean-zero. In the context of the partially linear model, this condition is satisfied for the moment condition (ref):

where the derivative is with respect to the scalar $r$ and evaluated at $r=0$. Heuristically, we can see that this condition alleviates the impact of noisy estimation of nuisance functions as local deviations of the nuisance functions away from their true values leave the moment condition unchanged. We refer to Chernozhukov2018 for a detailed discussion but highlight that all score functions discussed in this article are Neyman orthogonal.

Cross-fitting ensures independence between the estimation error from the first step and the regression residual in the second step. To implement cross-fitting, we randomly split the sample into $K$ evenly-sized folds, denoted as $I_1,\ldots, I_K$. For each fold $k$, the conditional expectations $\ell_0$ and $m_0$ are estimated using only observations not in the $k$th fold -- i.e., in $I^c_k\equiv I \setminus I_k$ -- resulting in $\hat{\ell}_{I^c_{k}}$ and $\hat{m}_{I^c_{k}}$, respectively, where the subscript ${I^c_{k}}$ indicates the subsample used for estimation. The out-of-sample predictions for an observation $i$ in the $k$th fold are then computed via $\hat{\ell}_{I^c_{k}}(\bm{X}_i)$ and $\hat{m}_{I^c_{k}}(\bm{X}_i)$. Repeating this procedure for all $K$ folds then allows for computation of the DDML estimator for $\theta_0$:

where $k_i$ denotes the fold of the $i$th observation.\footnote{We here omit the constant from the estimation stage. Since the residualized outcome and treatment may not be exactly mean-zero in finite samples, \stcmd{ddml}\xspace includes the constant by default in the estimation stage of partially linear models.}

We summarize the DDML algorithm for the Partially Linear Model in Algorithm 1:\footnote{Algorithm 1 corresponds to the `DML2' algorithm in Chernozhukov2018. Chernozhukov2018 recommend `DML2' over the alternative `DML1' algorithm, which fits the final estimator by fold.}

Chernozhukov2018 give conditions on the joint distribution of the data, in particular on $g_0$ and $m_0$, and properties of the nonparametric estimators used for CEF estimation, such that $\hat{\theta}_n$ is consistent and asymptotically normal. Standard errors are equivalent to the conventional linear regression standard errors of $Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)$ on $D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)$. ddml computes the DDML estimator for the Partially Linear Model using Stata's regress. All standard errors available for linear regression in Stata are also available in \stcmd{ddml}\xspace, including different heteroskedasticity and cluster-robust standard errors.\footnote{See help regress\#\#vcetype for available options.}

\paragraph{Remark 1: Number of folds.} The number of cross-fitting folds is a necessary tuning choice. Theoretically, any finite value is admissable. Chernozhukov2018 report that four or five folds perform better than only using $K=2$. Based on our simulation experience, we find that more folds tends to lead to better performance as more data is used for estimation of conditional expectation functions, especially when the sample size is small. We believe that more work on setting the number of folds would be useful, but believe that setting $K = 5$ provides is likely a good baseline in many settings.

\paragraph{Remark 2: Cross-fitting repetitions.} DDML relies on randomly splitting the sample into $K$ folds. We recommend running the cross-fitting procedure more than once using different random folds to assess randomness introduced via the sample splitting. \stcmd{ddml}\xspace facilitates this using the rep(integer) options, which automatically estimates the same model multiple times and combines the resulting estimates to obtain the final estimate. By default, \stcmd{ddml}\xspace reports the median over cross-fitting repetitions. \stcmd{ddml}\xspace also supports the average of estimates. Specifically, let $\hat\theta_n^{(r)}$ denote the DDML estimate from the $r$th cross-fit repetition and $\hat{s}_n^{(r)}$ its associated standard error estimate with $r=1,\ldots,R$. The aggregate median point estimate and associated standard error are defined as \[ \breve{\hat{\theta}}_n = {\textrm{median}}{\left(\left(\hat\theta_n^{(r)}\right)_{r=1}^R\right)} \quad\textrm{and}\quad \breve{\hat{s}}_n = \sqrt{\textrm{median}{\left(\left((\hat{s}_n^{(r)})^2 + (\hat\theta_n^{(r)} - \breve{\hat{\theta}}_n)^2\right)_{r=1}^R\right)}}. \] The aggregate mean point estimates and associated standard error are calculated as \[ \bar{\hat{\theta}}_n = \frac{1}{R}\sum_{r=1}^R \hat\theta_n^{(r)} \quad\textrm{and}\quad \bar{\hat{s}}_n = \sqrt{\textrm{hmean}{\left(\left((\hat{s}_n^{(r)})^2 + (\hat\theta_n^{(r)} - \bar{\hat{\theta}}_n)^2\right)_{r=1}^R\right)}}, \] where $\textrm{hmean}{()}$ is the harmonic mean.\footnote{The harmonic mean of $x_1,\ldots,x_n$ is defined as $\textrm{hmean}{(x_1,\ldots,x_n)}=n\left(\sum_{i=1}^n \frac{1}{x_i} \right)^{-1}$. We use the harmonic mean as it is less sensitive to outlier values.}

\paragraph{Remark 3: Cluster-dependence and folds.} Under cluster-dependence, we recommend randomly assigning folds by cluster; see fcluster(varname).

The Interactive Model (interactive)

The Interactive Model is given by

where $D$ takes values in $\{0,1\}$. The key deviations from the Partially Linear Model are that $D$ must be a scalar binary variable and that $D$ is not required to be additively separable from the controls $\bm{X}$. In this setting, the parameters of interest we consider are

which correspond to the average treatment effect (ATE) and average treatment effect on the treated (ATET), respectively.

Assumptions (ref) and (ref) below are sufficient for identification of the ATE and ATET. Note the conditional mean independence condition stated here is stronger than the conditional orthogonality assumption sufficient for identification of $\theta_0$ in the Partially Linear Model.

Under assumptions (ref) and (ref), we have

so that identification of the ATE and ATET immediately follows from their definition.\footnote{In the defined Interactive Model under Assumption (ref), the heterogeneity in treatment effects that the ATE and ATET average over is fully observed since $U$ is additively separable. Under stronger identifying assumptions, the DDML ATE and ATET estimators outlined here also apply to the ATE and ATET in the general causal model (ref) that average over both observed and unobserved heterogeneity. See, e.g., Belloni2017.}

In contrast to Section (ref), second-step estimators are not directly based on the moment conditions used for identification. Additional care is needed to ensure local robustness to first-stage estimation errors (i.e., Neyman orthogonality). In particular, the Neyman orthogonal score for the ATE that Chernozhukov2018 consider is the efficient influence function of Hahn1998

where $\bm{W}\equiv (Y, D, \bm{X})$. Similarly for the ATET,

Importantly, for $g_0(D,\bm{X})\equiv E[Y\vert D, \bm{X}]$, $m_0(\bm{X})\equiv E[D\vert \bm{X}]$, and $p_0 \equiv E[D]$, Assumptions (ref) and (ref) imply

and we also have that the Gateaux derivative of each condition with respect to the nuisance parameters $(g_0,m_0,p_0)$ is zero.

As before, the DDML estimators for the ATE and ATET leverage cross-fitting. The DDML estimators of the ATE and ATET based on $\psi^{\textrm{ATE}}$ and $\psi^{\textrm{ATET}}$ are

where $\hat{g}_{I^c_{k}}$ and $\hat{m}_{I^c_{k}}$ are cross-fitted estimators for $g_0$ and $m_0$ as defined in Section (ref). Since $D$ is binary, the cross-fitted values $\hat{g}_{I^c_{k}}(1,\bm{X})$ and $\hat{g}_{I^c_{k}}(0,\bm{X})$ are computed by only using treated and untreated observations, respectively. $\hat{p}_{I^c_{k}}$ is a cross-fitted estimator of the unconditional treatment probability.

\stcmd{ddml}\xspace supports heteroskedasticity and cluster-robust standard errors for $\hat{\theta}_n^{\textrm{ATE}}$ and $\hat{\theta}_n^{\textrm{ATET}}$. The algorithm for estimating the ATE and ATET are conceptually similar to Algorithm 1. We delegate the detailed outline to Algorithm A.1 in the Appendix. Mean and median aggregation over cross-fitting repetitions are implemented as outlined in Remark 2.

DDML with Instrumental Variables

This section outlines the Partially Linear IV Model, the Flexible Partially Linear IV Model, and the Interactive IV Model. The discussion is again based on Chernozhukov2018. As in the previous section, each model is a special case of the general causal model (ref). The discussion in this section differs from the preceding section in that identifying assumptions leverage instrumental variables $\bm{Z}$. The two partially linear IV models assume strong additive separability as in (ref), while the Interactive IV Model allows for arbitrary interactions between the treatment $D$ and the controls $\bm{X}$ as in (ref). The Flexible Partially Linear IV Model allows for approximation of optimal instruments\footnote{We only accommodate approximation of optimal instruments under homoskedasticity. The instruments are valid more generally but are not optimal under heteroskedasticity. Obtaining optimal instruments under heteroskedasticity would require estimating conditional variance functions.} as in Belloni2012 and Chernozhukov2015a, but relies on a stronger independence assumption than the Partially Linear IV Model. Throughout this discussion, we consider a random sample $\{(Y_i, D_i, \bm{X}_i,\bm{Z}_i)\}_{i=1}^n$ from $(Y, D, \bm{X}, \bm{Z})$.

Partially Linear IV Model (iv)

The Partially Linear IV Model considers the same functional form restriction on the causal model as the Partially Linear Model in Section (ref). Specifically, the Partially Linear IV Model maintains

where $\theta_0$ is the unknown parameter of interest.\footnote{As in Section (ref), the interpretation of $\theta_0$ can be generalized under stronger identifying assumptions. See angrist2000interpretation.}

The key deviation from the Partially Linear Model is that the identifying assumptions leverage instrumental variables $Z$, instead of directly restricting the dependence of $D$ and $U$. For ease of exposition, we focus on scalar-valued instruments in this section but we emphasize that \stcmd{ddml}\xspace for Partially Linear IV supports multiple instrumental variables and multiple treatment variables.

Assumptions (ref) and (ref) below are sufficient orthogonality and relevance conditions, respectively, for identification of $\theta_0$.

To show identification, consider the score function

where $\bm{W}\equiv (Y, D, \bm{X}, Z)$. Note that for $\ell_0(\bm{X}) \equiv E[Y\vert \bm{X}]$, $m_0(\bm{X})\equiv E[D\vert \bm{X}]$, and $r_0(\bm{X})\equiv E[Z\vert \bm{X}]$, Assumption (ref) implies $E[\psi(\bm{W}; \theta_0, \ell_0,m_0,r_0)]=0.$ We will also have that the Gateux derivative of $E[\psi(\bm{W}; \theta_0, \ell_0,m_0,r_0)]$ with respect to the nuisance functions $(\ell_0,m_0,r_0)$ will be zero. Rewriting $E[\psi(\bm{W}; \theta_0, \ell_0,m_0,r_0)]=0$ then results in a Wald expression given by

where Assumption (ref) is used to ensure a non-zero denominator.

The DDML estimator based on Equation (ref) is given by

where $\hat{\ell}_{I^c_{k}}$, $\hat{m}_{I^c_{k}}$, and $\hat{r}_{I^c_{k}}$ are appropriate cross-fitted CEF estimators.

Standard errors corresponding to $\hat{\theta}_{n}$ are equivalent to the IV standard errors where $Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)$ is the outcome, $D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)$ is the endogenous variable, and $Z_i-\hat{r}_{I^c_{k_i}}(\bm{X}_i)$ is the instrument. \stcmd{ddml}\xspace supports conventional standard errors available for linear instrumental variable regression in Stata, including heteroskedasticity and cluster-robust standard errors. Mean and median aggregation over cross-fitting repetitions are implemented as outlined in Remark 2. In the case where we have multiple instruments or endogenous regressors, we adjust the algorithm by residualizing each instrument and endogenous variable as above and applying two-stage least squares with the residualized outcome, endogenous variables, and instruments.

Flexible Partially Linear IV Model (fiv)

The Flexible Partially Linear IV Model considers the same parameter of interest as the Partially Linear IV Model. The key difference here is that identification is based on a stronger independence assumption which allows for approximating optimal instruments using nonparametric estimation, including machine learning, akin to Belloni2012 and Chernozhukov2015a. In particular, the Flexible Partially Linear IV Model leverages a conditional mean independence assumption rather than an orthogonality assumption as in Section (ref). As in Section (ref), we state everything in the case of a scalar $D$.

Assumption (ref) implies that for any function $\tilde{p}(\bm{Z},\bm{X})$, it holds that

where $\ell_0(\bm{X})=E[Y\vert \bm{X}]$ and $m_0(\bm{X})=E[D\vert \bm{X}]$. Identification based on (ref) requires that there exists some function $\tilde{p}$ such that

A sufficient assumption is that $D$ and $\bm{Z}$ are not mean independent conditional on $\bm{X}$. This condition allows setting $\tilde{p}(\bm{Z},\bm{X}) = E[D\vert \bm{Z},\bm{X}]$ which will then satisfy (ref).\footnote{The choice $\tilde{p}(\bm{Z},\bm{X}) = E[D\vert \bm{Z},\bm{X}]$ results in the optimal instrument, in the sense of semi-parametric efficiency, under homoskedasticity.} Assumption (ref) is a consequence of this non-mean independence.

Consider now the score function

where $\bm{W}\equiv (Y, D, \bm{X}, \bm{Z})$. Note that for $\ell_0(\bm{X}) \equiv E[Y\vert \bm{X}]$, $m_0(\bm{X})\equiv E[D\vert \bm{X}]$, and $p_0(\bm{Z},\bm{X}) \equiv E[D\vert \bm{Z},\bm{X}]$, Assumption (ref) and the law of iterated expectations imply $E[\psi(\bm{W}; \theta_0, \ell_0, m_0, p_0)]=0$ and the Gateaux differentiability condition holds. Rewriting then results in a Wald expression given by

where Assumption (ref) ensures a non-zero denominator.

The DDML estimator based on the moment solution (ref) is given by

where $\hat{\ell}_{I^c_{k}}$, $\hat{m}_{I^c_{k}}$, and $\hat{p}_{I^c_{k}}$ are appropriate cross-fitted CEF estimators.

In simulations, we find that the finite sample performance of the estimator in (ref) improves when the law of iterated expectations applied to $E[p_0(\bm{Z},\bm{X})] = m_0(\bm{X})$ is explicitly approximately enforced in estimation. As a result, we propose an intermediate step to the previously considered two-step DDML algorithm: Rather than estimating the conditional expectation of $D$ given $\bm{X}$ directly, we estimate it by projecting first-step estimates of the conditional expectation of $p_0(\bm{Z},\bm{X})$ onto $\bm{X}$ instead. Algorithm 2 outlines the LIE-compliant DDML algorithm for computation of (ref).

Standard errors corresponding to $\hat{\theta}_{n}$ in (ref) are the same as in Section (ref) where the instrument is now given by $\hat{p}_{I^c_{k_i}}(\bm{Z}_i, \bm{X}_i)-\hat{m}_{I^c_{k_i}}(\bm{X}_i)$. Mean and median aggregation over cross-fitting repetitions are as outlined in Remark 2.

Interactive IV Model (interactiveiv)

The Interactive IV Model considers the same causal model as in Section (ref); specifically

where $D$ takes values in $\{0,1\}$. The key difference from the Interactive Model is that this section considers identification via a binary instrument $Z$ representing assignment to treatment.

The parameter of interest we target is

where $p_0(Z,\bm{X})\equiv \Pr(D=1\vert Z, \bm{X})$. Here, $\theta_0$ is a local average treatment effect (LATE). Note that in contrast to the LATE developed in imbens1994identification, we follow the exposition in Chernozhukov2018 where “local” does not strictly refer to compliers but instead observations with a higher propensity score -- i.e., a higher probability of complying.\footnote{Identification of the conventional complier-focused LATE is achieved under stronger conditional independence and monotonicity assumptions not introduced in this paper. Under these stronger assumptions, the DDML LATE estimator outlined here targets the conventionally considered LATE parameter.}

Identification again leverages Assumptions (ref) and (ref) made in the context of the Flexible Partially Linear IV Model. In addition, we assume that the propensity score is weakly monotone with probability one, and that the support of the instrument is independent of the controls.

Assumptions (ref)-(ref) imply that

where $\ell_0(Z, \bm{X}) \equiv E[Y\vert Z, \bm{X}]$, verifying identification of the LATE $\theta_0$. Akin to Section (ref), however, estimators of $\theta_0$ should not directly be based on Equation (ref) because the estimating equations implicit in obtaining (ref) do not satisfy Neyman-orthogonality. Hence, a direct estimator of $\theta_0$ obtained by plugging nonparametric estimators in for nuisance functions in (ref) will potentially be highly sensitive to the first step nonparametric estimation error. Rather, we base estimation on the Neyman orthogonal score function

where $\bm{W} \equiv (Y, D, \bm{X}, Z)$ . Note that under Assumptions (ref)-(ref) and for $\ell_0(Z,\bm{X}) \equiv E[Y\vert Z,\bm{X}]$, $p_0(Z,\bm{X}) \equiv E[D|Z,\bm{X}]$, and $r_0(\bm{X})\equiv E[Z\vert \bm{X}]$, we have $E[\psi(\bm{W};\theta_0, \ell_0, p_0, r_0)]=0$ and can verify that its Gateaux derivative with respect to the nuisance functions local to their true values is also zero.

The DDML estimator based on the orthogonal score $\psi$ is then

where $\hat{\ell}_{I^c_{k}}$, $\hat{p}_{I^c_{k}}$, and $\hat{r}_{I^c_{k}}$ are appropriate cross-fitted CEF estimators. Since $Z$ is binary, the cross-fitted values $\hat{\ell}_{I^c_{k}}(1,\bm{X})$ and $\hat{p}_{I^c_{k}}(1,\bm{X})$, as well as $\hat{\ell}_{I^c_{k}}(0,\bm{X})$ and $\hat{p}_{I^c_{k}}(0,\bm{X})$ are computed by using only assigned and unassigned observations, respectively.

\stcmd{ddml}\xspace supports heteroskedasticity and cluster-robust standard errors for $\hat{\theta}_n$. Mean and median aggregation over cross-fitting repetitions are implemented as outlined in Remark 2.

The choice of machine learner

Chernozhukov2018 show that DDML estimators are asymptotically normal when used in combination with a general class of machine learners satisfying a relatively weak convergence rate requirement for estimating the CEFs. While asymptotic properties of common machine learners remain an highly active research area, recent advances provide convergence rates for special instances of many machine learners, including lasso bickel2009simultaneous, Belloni2012, random forests wager2015adaptive, wager2018, Athey2019a, neural networks schmidt2020nonparametric, farrell2021deep, and boosting luo2016high. It seems likely that many popular learners will fall under the umbrella of suitable learners as theoretical results are further developed. However, we note that currently known asymptotic properties do not cover a wide range of learners, such as very deep and wide neural networks and deep random forests, as they are currently implemented in practice.

The relative robustness of DDML to the first-step learners leads to the question of which machine learner is the most appropriate for a given application. It is ex ante rarely obvious which learner will perform best. Further, rather than restricting ourselves to one learner, we might want to combine several learners into one final learner. This is the idea behind stacking generalization, or simply “stacking”, due to Wolpert1992 and Breiman1996a. Stacking allows one to accommodate a diverse set of base learners with varying tuning and hyper-tuning parameters. It thus provide a convenient framework for combining and identifying suitable learners, thereby reducing the risk of misspecification. Ahrens2023_applied introduce short-stacking which reduces the computational cost of pairing DDML and stacking drastically, as well as pooled stacking which enforces common weights across cross-fitting folds.

We discuss stacking approaches to DDML estimation in Section (ref). Section (ref) demonstrates the performance of DDML in combination with stacking approaches using a simulation.

DDML and stacking

Our discussion of stacking in the context of DDML focuses on the Partially Linear Model in (ref), but we highlight that DDML and stacking can be combined in the same way for all other models supported in \stcmd{ddml}\xspace. Suppose we consider $J$ machine learners, referred to as base learners, to estimate the CEFs $\ell_0(\bm{X}) \equiv E[Y|\bm{X}]$ and $m_0(\bm{X}) \equiv E[D|\bm{X}]$. The set of base learners could, for example, include cross-validated lasso and ridge with alternative sets of predictors, gradient boosted trees with varying tree depth and feed-forward neural nets with varying number of hidden layers and neurons. Generally, we recommend considering a relatively large and diverse set of base learners, and including some learners with alternative tuning parameters.

We randomly split the sample into $K$ cross-fitting folds, denoted as $I_1,\ldots,I_K$. In each cross-fitting step $k$, we define the training sample as $I_k^c\equiv T_k$, comprising all observations excluding the cross-fitting hold-out fold $k$. This training sample is further divided into $V$ cross-validation folds, denoted as $T_{k,1},\ldots,T_{k,V}$. The stacking regressor fits a final learner to the training sample $T_k$ using the cross-validated predicted values of each base learner as inputs. A typical choice for the final learner is constrained least squares (CLS) which restricts the weights to be positive and sum to one. The stacking objective function for estimating $\ell_0(\bm{X})$ using the training sample $T_k$ is then defined as:

where $w_{k,j}$ are referred to as stacking weights. We use $\hat{\ell}^{(j)}_{T_{k,v(i)}^c}(\bm{X}_i)$ to denote the cross-validated predicted value for observation $i$, which is obtained from fitting learner $j$ on the sub-sample $T_{k,v(i)}^c \equiv T_k \setminus T_{k,v(i)}$, i.e., the sub-sample excluding the fold $v(i)$ into which observation $i$ falls. The stacking predicted values are obtained as $\sum_j \hat{w}_{k,j} \hat{\ell}_k^{(j)}(\bm{X}_i)$ where each learner $j$ is fit on the step-$k$ training sample $T_k$. The objective function for estimating $m_0(\bm{X})$ is defined accordingly.

CLS frequently performs well in practice and facilitates the interpretation of stacking as a weighted average of base learners Hastie2009. It is, however, not the only sensible choice of combining base learners. For example, stacking could instead select the single learner with the lowest quadratic loss, i.e., by imposing the constraint $w_{k,j}\in\{0,1\}$ and $\sum_{k,j}w_{k,j}=1$. We refer to this choice as “single best” and include it in our simulation experiments. We implement stacking for DDML using \stcmd{pystacked}\xspace Ahrens2022.

\paragraph{Pooled stacking.} A variant of stacking specific to DDML is pooled stacking. Standard stacking fits the final learner $K$ separate times, once in each cross-fitting step, yielding $K$ separate sets of stacking weights $\hat{w}_{k,j}$ for the $J$ learners. With DDML pooled stacking, we can impose the additional constraint in (ref) that the weights are the same across all cross-fit folds, $\hat{w}_{k,j}=\hat{w}_{j}, \forall \: k$. By returning a single set of stacking weights, pooled stacking imposes an additional degree of regularization and facilitates interpretation but suffers from the same high computational cost as pairing DDML with (regular) stacking.

\paragraph{Short-stacking.} Stacking and pooled stacking rely on cross-validation. In the context of DDML we can also exploit the cross{\it -fitted} predicted values directly for stacking. That is, we can directly apply CLS to the cross-fitted predicted values for estimating $\ell_0(\bm{X})$ (and similarly $m_0(\bm{X})$): \[ \underset{w_1,\ldots,w_J}{\min} \sum_{i=1}^n \left( Y_i - \sum_{j=1}^J w_j\hat{\ell}^{(j)}_{I_{k(i)}^c}(\bm{X}_i) \right)^2, \qquad \textrm{s.t.}\ w_j\geq 0,\ \sum_{j=1}^J |w_j|=1\] We refer to this form of stacking that utilizes the cross-fitted predicted values as short-stacking as it takes a short-cut. This is to contrast it with regular stacking which estimates the stacking weights for each cross-fitting fold $k$. The main advantage of short-stacking relative to standard stacking is the lower computational cost, because short-stacking does not require the fitting of the $j$ learners on sub-samples to obtain the cross-validated predicted values $\hat{\ell}^{(j)}_{T_{k,v(i)}^c}(\bm{X}_i)$ needed for standard stacking. Furthermore, short-stacking (like pooled stacking) also produces a single set of weights for the entire sample, which facilitates interpretation and implies a higher degree of regularization. A potential disadvantage of short-stacking is that it is more susceptible to over-fitting issues since stacking weights and structural parameters are estimated using the same cross-fitted predicted values. We thus recommend only considering short-stacking in regular settings where the number of candidate learners is small relative to $N$ Ahrens2023_applied. Algorithm A.4 in the Appendix summarizes the short-stacking algorithm for the Partially Linear Model.\footnote{While short-stacking can be applied in a similar fashion to other conditional expectations, a complication arises in the Flexible Partially Linear IV Model where the cross-fitted predicted values of $E[D|\bm{X}]$ depend on $E[D|\bm{X},\bm{Z}]$. We describe the algorithm that accounts for this in the Appendix; see Algorithm A.5.}

Monte Carlo simulation

To illustrate the advantages of DDML with stacking, we generate artificial data based on the Partially Linear Model

where both $\varepsilon_i$ and $u_i$ are independently drawn from the standard normal distribution. We set the target parameter to $\theta_0=0.5$ and the sample size to either $n=100$ or $n=1000$. The controls $\bm{X}_i$ are drawn from the multivariate normal distribution with $N(0,\bm{\Sigma})$ where $\Sigma_{ij}=(0.5)^{|i-j|}$. The number of controls is set to $p=\dim(\bm{X}_i)=50$, except in DGP 5 where $p=7$. The constants $c_Y$ and $c_D$ are chosen such that the $R^2$ in (ref) and (ref) are approximately equal to 0.5. To induce heteroskedasticity, we set \[\sigma_D\left(\bm{X}_i\right)=\sqrt{\frac{\left(1+g(\bm{X}_i)\right)^2}{\frac{1}{n}\sum_i\left(1+g(\bm{X}_i)\right)^2}} \quad\text { and } \quad \sigma_Y\left(D_i, \bm{X}_i\right)=\sqrt{\frac{\left(1+\theta_0 D_i+g(\bm{X}_i)\right)^2}{\frac{1}{n}\sum_i\left(1+\theta_0 D_i+g(\bm{X}_i)\right)^2}}\] The nuisance function $g(\bm{X}_i)$ is generated using five exemplary DGPs, which cover linear and nonlinear processes with varying degrees of sparsity and varying number of observed covariates: \[

\] DGP 1 is a linear design involving many negligibly small parameters. While not exactly sparse, the design can be approximated well through a sparse representation. DGP 2 is linear in the parameters and exactly sparse, but includes interactions and second-order polynomials. DGPs 3-5 are also exactly sparse but involve complex nonlinear and interaction effects. DGP 4 and 5 are identical, except that DGP 5 does not add nuisance covariates that are unrelated to $Y$ and $D$.

We consider DDML with the following supervised machine learners for cross-fitting the CEFs:\footnote{All base learners have been implemented using \stcmd{pystacked}\xspace. We use the defaults of \stcmd{pystacked}\xspace for parameter values and settings not mentioned here.}

We use the above set of learners as base learners for DDML with stacking approaches. Specifically, we estimate DDML using stacking, short-stacking and pooled stacking which we combine with CLS and the single-best learner. We set the number of folds to $K=20$ if $n=100$, and $K=5$ if $n=1000$. That is, we adapt the number of folds $K$ to the total sample size $n$ to ensure that the CEF estimators are trained on sufficiently large training samples.

For comparison, we report results for OLS and PDS-lasso with base controls, PDS-Lasso with Poly 5, PDS-Lasso with Poly 2 + Interactions, and an oracle estimator using the full sample.\footnote{The PDS-Lasso estimators set tuning parameters using the default in \stcmd{pdslasso}\xspace.} The oracle estimator presumes knowledge of the function $g(\bm{X})$, and obtains estimates by regressing $Y$ on the two variables $D$ and $g(\bm{X})$.

We report simulation median absolute bias (MAB) and coverage rates of 95% confidence intervals (CR) for DGPs 1-3 in Table (ref). We delegate results for DGPs 4 and 5, including a brief discussion, to Appendix B. DDML estimators leveraging stacking approaches perform favorably in comparison to individual base learners in terms of bias and coverage. The relative performance of stacking approaches seems to improve as the sample size increases, likely reflecting that the stacking weights are more precisely estimated in larger small samples. For $n=1000$, the bias of stacking with CLS is at least as low as the bias of the best-performing individual learner under DGP 1-2, while only gradient boosting and neural net yield a lower bias than stacking under DGP 3.

Results for coverage are similar with stacking-based estimates being comparable with the best performing feasible estimates and the oracle when $n = 1000.$ With $n = 100$, coverage of confidence intervals for stacking-based estimators are inferior to coverages for a small number of the individual learners but still competitive and superior than most learners. Looking across all results, we see that stacking provides robustness to potentially very bad performance that could be obtained from using a single poorly performing learner.

There are overall little performance differences among the six stacking estimators considered., suggesting that short-stacking has a substantial practical advantage due to its lower computational cost. Ahrens2023_applied report that short-stacking reduces the compute time by a factor of $1/V$ where $V$ is the number of cross-validation folds. There is some evidence that the single-best selector outperforms CLS in very small sample sizes in DGPs 2-3, but not in DGP 1 (and also not in DGPs 4-5, see Table B.1). We suspect that the single-best selector works better in scenarios where there is one base learner that clearly dominates.

The mean-squared prediction errors (MSPE) and the average stacking weights, which we report in Tables B.2 and B.3 in the Appendix, provide further insights into how stacking with CLS functions. CLS assigns large stacking weights to base learners with a low MSPE, which in turn are associated with a low bias. Importantly, stacking assigns zero or close-to-zero weights to poorly specified base learners such as the highly regularized random forest, which in all three DGPs ranks among the individual learners with highest MSPE and highest bias. The robustness to misspecified and ill-chosen machine learners, which could lead to misleading inference, is indeed one of our main motivations for advocating stacking approaches to DDML.

DDML with stacking approaches also compares favorably to conventional full-sample estimators. In the relatively simple linear DGP 1, DDML with stacking performs similarly to OLS and the infeasible oracle estimator---both in terms of bias and coverage---for $n=100$ and $n=1000$. In the more challenging DGPs 2 and 3, the bias of DDML with stacking is substantially lower than the biases of OLS and the PDS-Lasso estimators. While the bias and size distortions of DDML with stacking are still considerable in comparison to the infeasible oracle for $n=100$, they are close to the oracle for $n=1000$. The results overall highlight the flexibility of DDML with stacking to flexibly approximate a wide range of DGPs provided a diverse set of base learners is chosen.

The program

In this section, we provide an overview of the \stcmd{ddml}\xspace package. We introduce the syntax and workflow for the main programs in Section (ref). Section (ref) lists the options. Section (ref) covers the simplified one-line program \stcmd{qddml}\xspace. We provide an overview of supported machine learning programs in Section (ref). Finally, Section (ref) adds a note on how to ensure replication with \stcmd{ddml}\xspace.

Syntax: \stcmd{ddml}\xspace

The ddml estimation proceeds in four steps.

Step 1: Initialize \stcmd{ddml}\xspace and select model.

where {\it model} selects between the Partially Linear Model (partial), the Interactive Model (interactive), the Partially Linear IV Model (iv), the Flexible Partially Linear IV Model (fiv), and the Interactive IV Model (interactiveiv). This step creates a persistent Mata object with the name provided by mname(\emph{name}\texttt{)} in which model specifications and estimation results will be stored. The default name is {\it m0}.

At this stage, the user-specified folds for cross-fitting can be set via integer-valued Stata variables (see {\tt foldvar({\it varlist})}). By default, observations are randomly assigned to folds and {\tt kfolds({\it integer})} determines the number of folds (the default is 5). Cluster-randomized fold splitting is supported (see {\tt fcluster({\it varname})}). The user can also select the number of times to fully repeat the cross-fitting procedure (see {\tt rep({\it integer})}).

Step 2: Add supervised machine learners for estimating conditional expectations.

In this second step, we select the machine learning programs for estimating CEFs.

where {\it cond_exp} selects the conditional expectation to be estimated by the machine learning program command. At least one learner is required for each conditional expectation. Table (ref) provides an overview of which conditional expectations are required by each model. The program command is a supervised machine learning program such as cvlasso or \stcmd{pystacked}\xspace (see compatible programs in Section (ref)). The options {\it cmdopt} are specific to that program.

Step 3: Perform cross-fitting.

This step implements the cross-fitting algorithm. Each learner is fit iteratively on training folds and out-of-sample predicted values are obtained. Cross-fitting is the most time-consuming step, as it involves fitting the selected machine learners repeatedly.

Step 4: Estimate causal effects.

In the last step, we estimate the parameter of interest for all combination of learners added in Step 2.

To report and post selected results, we can use {\tt ddml estimate} with the {\tt replay} option:

Utilities

{\tt ddml describe} provides information about the model setup and or results:

{\tt ddml} stores many internal results on associate arrays, notably the various stacking weights. These can be retrieved using {\tt ddml extract:}

{\tt ddml export} saves the estimated conditional expectations and other variables to a CSV file:

Options

\paragraph{Step 1 options: Initialization.}

\paragraph{Step 2 options: Adding learners.}

\paragraph{Step 3 options: Cross-fitting.}

\paragraph{Step 4 options: Estimation.}

Re-fitting the final learner using {\tt ddml estimate} with stacking options is generally very fast because it does not require re-cross-fitting.

Short syntax: \stcmd{qddml}\xspace

The \stcmd{ddml}\xspace package includes the wrapper program \stcmd{qddml}\xspace which provides a one-line syntax for estimating a \stcmd{ddml}\xspace model. The one-line syntax follows the syntax of pdslasso and ivlasso Ahrens2018. The main restriction of \stcmd{qddml}\xspace compared to the more flexible multi-line syntax is that \stcmd{qddml}\xspace only allows for one user-specified machine learner (in addition to regress, which is added by default).

\stcmd{qddml}\xspace has integrated support for \stcmd{pystacked}\xspace, and \stcmd{pystacked}\xspace, is the default learner in all equations. The syntax for \stcmd{qddml}\xspace options differs depending on whether \stcmd{pystacked}\xspace is used as the learner in each equation.

Syntax when used with \stcmd{pystacked}\xspace

The \stcmd{pystacked}\xspace option sets the options for all the conditional expectations estimated by \stcmd{pystacked}\xspace; the _y, _d and _z variants control the options sent to the corresponding conditional expectation estimations. Other options are as in \stcmd{ddml}\xspace.

Syntax when used with other learners

The cmd option sets the options for all the conditional expectations estimated by \stcmd{pystacked}\xspace; the ycmd, dcmd and zcmd variants control the options sent to the corresponding conditional expectation estimations. The cmdopt option can be used either to set the options for all equations, or, by replacing the asterisk with y, d or z, set the options for the corresponding conditional expectation estimation. Other options are as in \stcmd{ddml}\xspace.

Supported machine learning programs

\stcmd{ddml}\xspace is compatible with any supervised ML program in Stata that supports the typical “reg\,\,y\,\,x” syntax, comes with a post-estimation \stcmd{predict} and supports \stcmd{if} statements. We have tested \stcmd{ddml}\xspace with the following programs:

Stata programs that are currently not supported can be added relatively easily using wrapper programs (see parsnip2 for an example).

Inspecting results and replication

In this section we discuss how to ensure replicability when using \stcmd{ddml}\xspace. We also discuss some tools available for tracing replication failures. First, however, we briefly describe how \stcmd{ddml}\xspace stores results.

\stcmd{ddml}\xspace stores estimated conditional expectations in Stata's memory as Stata variables. These variables can be inspected, graphed and summarized in the usual way. Fold ID variables are also stored as Stata variables (by default named \verb+m0_fid_+r, where \verb+m0+ is the default model name and r is the cross-fitting repetition). \stcmd{ddml}\xspace models are stored on Mata structs and using Mata's associative arrays. Specifically, the \stcmd{ddml}\xspace model created by \stcmd{ddml}\xspace init is an mStruct, and information relating to the estimation of conditional expectations are stored in eStructs. Results relating to the overall model estimation are stored in associative arrays that live in the \texttt{mStruct}, and results relating to the estimation of conditional expectations are stored in associative arrays that live in the corresponding \texttt{eStruct}s.

Replication tips:

Applications

We demonstrate the ddml workflow using two applications. In Section (ref), we apply the DDML estimator to estimate the effect of 401(k) eligibility on financial wealth following poterba1995. We focus on the Partially Linear Model for the sake of brevity, but provide code that demonstrates the use of \stcmd{ddml}\xspace with the Interactive Model, Partially Linear IV Model and Interactive IV Model using the same application in Appendix C. Additional examples can also be found in the help file. Based on berry1995automobile, we show in Section (ref) how to employ \stcmd{ddml}\xspace for the estimation of the Flexible Partially Linear IV Model which allows both for flexibly controlling for confounding factors using high-dimensional function approximation of confounding factors and for estimation of optimal instrumental variables.

401(k) and financial wealth

The data consists of $n=9915$ households from the 1991 SIPP. The application is originally due to poterba1995, but has been revisited by Belloni2017, Chernozhukov2018, and wuthrich2021, among others. Following previous studies, we include the control variables age, income, years of education, family size, as well as indicators for martial status, two-earner status, benefit pension status, IRA participation, and home ownership. The outcome is net financial assets and the treatment is eligibility to enroll for the 401(k) pension plan.

We load the data and define three globals for outcome, treatment and control variables. We then proceed in the four steps outlined in Section (ref).

Step 1: Initialize ddml model.

We initialize the \stcmd{ddml}\xspace model and select the Partially Linear Model in (ref). Before initialization, we set the seed to ensure replication. This should always be done before ddml init, which executes the random fold assignment. In this example, we opt for four folds to ensure the readability of some of the output shown below, although we recommend considering a larger number of folds in practice.

Step 2: Add supervised machine learners for estimating conditional expectations.

In this step, we specify which machine learning programs should be used for the estimation of the conditional expectations $E[Y\vert \bm{X}]$ and $E[D\vert \bm{X}]$. For each conditional expectation, at least one learner is required. For illustrative purposes, we consider regress for linear regression, pystacked with the m(lassocv) option for cross-validated lasso (as an example of how to use pystacked to estimate a single learner) and rforest for random forests. When using rforest, we need to add the option \texttt{vtype(none)} since the post-estimation \texttt{predict} command of \texttt{rforest} does not support variable types.

The flexible \stcmd{ddml}\xspace syntax allows specification of different sets of covariates for different learners. This flexibility can be useful as, for example, linear learners such as the lasso might perform better if interactions are provided as inputs, whereas tree-based methods such as random forests may detect certain interactions in a data-driven way. Here, we use interactions and second-order polynomials for the cross-validated lasso, but not for the other learners.

This application has only one treatment variable, but \stcmd{ddml}\xspace does support multiple treatment variables. To add a second treatment variable, we would simply add a statement such as ddml E[D|X]: reg D2 \$X where D2 would be the name of the second treatment variable. An example with two treatments is provided in the help file.

The auxiliary ddml sub-command describe allows us to verify that the learners were correctly registered:

Step 3: Perform cross-fitting.

The third step performs cross-fitting, which is the most time-intensive process. The shortstack option enables the short-stacking algorithm of Section (ref).

Six variables are created and stored in memory which correspond to the six learners specified in the previous step. These variables are called \verb+Y1_reg_1+, \verb+Y2_pystacked_1_1+, \verb+Y3_rforest_1+, \verb+D1_reg_1+, \verb+D2_pystacked_1_1+ and \verb+D3_rforest_1+. Y and D indicate outcome and treatment variable. The index 1 to 3 is a learner counter. reg, pystacked and \texttt{rforest} correspond to the names of the commands used. The \verb+_1+ suffix indicates the cross-fitting repetition. The additional \verb+_1+ in the case of \verb+D2_pystacked_1_1+ indicates the learner number (there is only a single \texttt{pystacked} learner).

After cross-fitting, we can inspect the mean-squared prediction errors by fold and learner:

Step 4: Estimate causal effects.

In this final step, we obtain the causal effect estimates. Since we requested short-stacking in Step 3, ddml shows the short-stacking result which relies on the cross-fitted values of each base learner. In addition, the specification that corresponds to the minimum-MSE learners is listed at the beginning of the output.

Since we have specified three learners per conditional expectation, there are in total 9 specifications relying on the base learners (since we can combine \verb+Y1_reg_1+, \verb+Y2_pystacked_1+ and \verb+Y3_rforest_1+ with \verb+D1_reg_1+, \verb+D2_pystacked_1+ and \verb+D3_rforest_1+). To get all results, we add the allcombos option:

We can use the spec(string) option to select among the listed specifications. {\it string} is either the specification number, ss, st or ps, to get the short-stacking, standard stacking or pooled stacking specification, respectively, or \texttt{mse} for the specification corresponding to the minimal MSPE. In the example above, \texttt{spec(1)} reports in full the specification using \texttt{regress} for estimating both $E[Y\vert \bm{X}]$ and $E[D\vert \bm{X}]$. The \texttt{spec(}\emph{string}\texttt{)} option can be provided either in combination with \texttt{allcombos}, or after estimation in combination with the {\tt replay} option, for example:

Manual final estimation.

In the background, ddml estimate regresses \verb+Y1_reg_1+ against \verb+D1_reg_1+ with a constant. We can verify this manually:

Manual estimation using { regress} allows the use of { regress}'s postestimation tools.

Stacking

We next demonstrate DDML with stacking. To this end, we exploit the stacking regressor implemented in \stcmd{pystacked}\xspace. \stcmd{pystacked}\xspace allows combining multiple base learners with learner-specific settings and covariates into a final meta learner. The learners are separated by \verb+||+. method(name) selects the learner, xvars(varlist) specifies learner-specific covariates (overwritting the default covariates \verb+$X+) and \texttt{opt(}\emph{string}\texttt{)} passes options to the learners. In this example, we use OLS, cross-validated lasso and ridge, random forests and gradient boosting. We furthermore use parallelization with 5 cores. A detailed explanation of the \stcmd{pystacked}\xspace syntax can be found in Ahrens2022.

After cross-fitting, we retrieve the cross-fitted MSPE using ddml extract with show(mse) or examine the stacking weights using stweights:

Finally, in the estimation stage, we retrieve the results of DDML with stacking:

The DDML-specific stacking approaches of short-stacking and pooled stacking can be requested at either the cross-fitting or estimation steps. Re-fitting the final learner at the estimation step allows us to avoid repeating the computationally-intensive cross-fitting. Here, we request short-stacking and pooled stacking but using the single-best base learner.

Short-stacking is computationally substantially faster than regular stacking (or pooled stacking) because it avoids the cross-validation within cross-fit folds. Below, we disable regular stacking with the nostdstack option in the cross-fitting stage. In this example, where we use parallelization with 5 cores, the run time is only 50.7 seconds DDML with short-stakcing compared to 93.0 seconds for DDML with regular stacking.

One-line syntax.

\stcmd{qddml}\xspace provides a simplified and convenient one-line syntax. The main constraint of \stcmd{qddml}\xspace is that it only allows for one user-specified learner. The default learner is pystacked, which by default uses OLS, cross-validated lasso and gradient boosting as default learners. The pystacked base learners as well as non-pystacked commands can be be modified via various options. Below we use qddml with pystacked's default base learners. We omit the output for the sake of brevity.

The market for automobiles

For this demonstration, we follow chernozhukov2015 who estimate a stylized demand model using instrumental variables based on the data from berry1995automobile. The authors of the original study estimate the effect of prices on the market share of automobile models in a given year ($n=2217$). The controls are product characteristics (a constant, air conditioning dummy, horsepower divided by weight, miles per dollar, vehicle size). To account for endogenous prices, berry1995automobile suggest exploiting other products' characteristics as instruments. Following chernozhukov2015, we define the baseline set of instruments as the sum over all other products' characteristics, calculated separately for own-firm and other-firm products, which yields 10 baseline instruments. chernozhukov2015 also construct an augmented set of instruments, including first-order interactions, squared and cubic terms. In the analysis below, we extend chernozhukov2015 by applying DDML with stacking and a diverse set of learners including OLS, lasso, ridge, random forest and gradient boosted trees. We use the augmented set of controls for all base learners and OLS, which we include for reference.

We load and prepare the data:

\paragraph{Step 1: Initialize ddml model.}

Note that in the ddml init step, we include the option reps(5) which will result in running the full cross-fitting procedure five times, each with a different random split of the data. Replicating the procedure multiple times allows us to gauge the impact of randomness due to the random splitting of the data into subsamples.

\paragraph{Step 2: Add supervised machine learners for estimating conditional expectations.}

Estimation of a fiv model requires us to add learners for $E[Y|\bm{X}]$, $E[D|\bm{X},\bm{Z}]$ and $E[D|\bm{X}]$. Compared to the other models supported by \stcmd{ddml}\xspace, there is one complication that arises because, in order to estimate $E[D|\bm{X}]$, we exploit fitted values of $E[D|\bm{X},\bm{Z}]$ to impose LIE-compliance. Since these fitted values have not yet been generated, we use the placeholder \verb+{D}+ that in the cross-fitting stage will be internally replaced with estimates of $E[D|\bm{X},\bm{Z}]$. We use the learner(string) option to match one learner for $E[D|\bm{X}]$ with a learner for $E[D|\bm{X},\bm{Z}]$, and vname(varname\texttt{)} to indicate the name of the treatment variable.

\paragraph{Steps 3-4: Perform cross-fitting (output omitted) and estimate causal effects.}

\paragraph{Manual final estimation.} We can obtain the final estimate manually. To this end, we construct the instrument as $\hat{E}[D|X,Z]-\hat{E}[D|X]$ and the residualized endogenous regressor as $D-\hat{E}[D|X]$. The residualized dependent variable is saved in memory. Here we obtain the estimate from the first cross-fitting replication. We could obtain the estimate for replication $r$ by changing the “_1” to “_$r$”.

Conclusion

This article introduces the Stata program ddml for Double/Debiased Machine Learning. The program allows for flexible estimation of structural parameters in five econometric models, leveraging a wide range of supervised machine learners. While ddml is compatible with many existing machine learning programs in Stata, it is specifically designed to be used with pystacked, which allows combining several learners into a meta learner via stacking. We see several avenues for extensions: First, ddml primarily focuses on cross-sectional models. Some panel models are readily implementable in ddml, but expanding its capabilities for seamless use with a wide range of panel models would increase its practical relevance. Second, researchers and policymakers are frequently interested in learning treatment effects for specific sub-populations sharing observable characteristics. The estimation of conditional average treatment effects would be a natural extension to the existing ddml program. Third, \texttt{ddml} currently lacks under-identification diagnostics and weak-identification robust inference for instrumental variable regressions, which we hope to add in future releases.

Acknowledgments

We thank users who tested earlier versions of the program. We also thank Jan Ditzen, Ben Jann, Eroll Kuhn, Di Liu and Moritz Marbach for helpful comments, as well as participants at Stata User Conferences in Germany (2021), Italy (2022), Switzerland (2022) and the UK (2023). All remaining errors are our own.