\section{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, \citet{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 \citet{roberts2020a} or digital trace data \citep{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 \citet{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 \citep{Belloni2012}. Lasso-based approaches for estimation of causal effects have become a popular strategy in applied econometrics \citep[e.g.][]{gilchrist2016something,dhar2022reshaping}, partially facilitated by the availability of software programs in Stata (\stcmd{pdslasso}\xspace, \citealp{Ahrens2018}; \citealp{statacorp2019}) and R (\texttt{hdm}, \citealp{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, \citet{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., \citealp{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, \citet{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 \citep{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 1.2 of \stcmd{ddml}\xspace.} \stcmd{ddml}\xspace adds to a small number of programs for causal machine learning in Stata (\citealp{Ahrens2018}, \citealp{statacorp2019}, \citealp{statacorp2021}). We briefly summarize the four main features of the program:
\begin{enumerate}

\item \stcmd{ddml}\xspace supports flexible estimators of causal parameters in five econometric models: (1) the Partially Linear Model, (2) the Interactive Model (for binary treatment), (3) the Partially Linear IV Model, (4) the Flexible Partially Linear IV Model, and (5) the Interactive IV Model (for binary treatment and instrument).

\item \stcmd{ddml}\xspace supports data-driven combinations of multiple machine learners via stacking by leveraging \texttt{pystacked} \citep{Ahrens2022}, our complementary Stata frontend relying on the Python library \emph{scikit-learn} \citep{scikit-learn,sklearn_api}.

\item Aside from \texttt{pystacked}, \stcmd{ddml}\xspace can be used in combination with many other existing supervised machine learning programs available in or via Stata. \stcmd{ddml}\xspace has been tested with \texttt{lassopack} \citep{Ahrens2019}, \texttt{rforest} \citep{Schonlau2020}, \texttt{svmachines} \citep{Guenther2018}, and \texttt{parsnip} \citep{Huntington2021}. Indeed, the requirements for compatibility with \stcmd{ddml}\xspace are minimal: Any \texttt{eclass} program with the Stata-typical ``\texttt{reg\,\,y\,\,x}'' syntax, support for \texttt{if} conditions and post-estimation \texttt{predict} is compatible with \stcmd{ddml}\xspace.

\item \stcmd{ddml}\xspace provides flexible multi-line syntax and short one-line syntax. The multi-line syntax offers a wide range of options, guides the user through the DDML algorithm step-by-step, and includes auxiliary programs for storing, loading and displaying additional information. We also provide a complementary one-line version called \stcmd{qddml}\xspace (`quick' \stcmd{ddml}\xspace), which uses a similar syntax as \stcmd{pdslasso}\xspace and \texttt{ivlasso} \citep{Ahrens2018}.
\end{enumerate}

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

\section{DDML with Conditional Unconfoundedness}\label{sec:ddml_confounding}

This section discusses DDML for the Partially Linear Model and the Interactive Model in turn. Both models are special cases of the general causal model \begin{align}\label{eq:all_causes_1}
    Y = f_0(D,\bm{X},U),
\end{align}
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 \eqref{eq:all_causes_1}, $(D, \bm{X}, U)$ jointly determine $Y$, the model is also dubbed the ``all causes model'' (see, e.g., \citealp{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.

\subsection{The Partially Linear Model (\texttt{partial})}\label{sec:plm}
The Partially Linear Model imposes the estimation model
\begin{align}\label{eq:plm}
\begin{aligned}
    Y & = \theta_0 D + g_0(\bm{X}) + U
    \end{aligned}
\end{align}
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 \citet{angrist1999empirical} imply that in the general causal model \eqref{eq:all_causes_1}, $\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{assumption:plm_orthogonal}.\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.}
\begin{assumption}[Conditional Orthogonality]\label{assumption:plm_orthogonal}
$E[Cov(U, D\vert \bm{X})]= 0$ .
\end{assumption}

To show identification of $\theta_0$, consider the score \begin{align}\label{eq:plm_score}
    \psi(\bm{W}; \theta, m,\ell) = \Big(Y - \ell(\bm{X}) - \theta\big(D - m(\bm{X})\big)\Big)\big(D - m(\bm{X})\big),
\end{align}
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
\begin{align*}
    E\left[\psi(\bm{W}; \theta_0, m_0, \ell_0)\right] = 0
\end{align*}
by Assumption \ref{assumption:plm_orthogonal}. When in addition $E[Var(D\vert \bm{X})] \neq 0$, we get
\begin{align}\label{eq:plm_alpha0}
    \theta_0 = \frac{E\left[\big(Y - \ell_0(\bm{X})\big)\big(D - m_0(\bm{X})\big)\right]}{E\left[(D - m_0(\bm{X}))^2\right]}.
\end{align}

Equation \eqref{eq:plm_alpha0} 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 \citet{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 \eqref{eq:plm_score}:
\begin{align*}
    \partial_r \left\{E\left[\psi(W; \theta_0, m_0 + r(m - m_0),\ell_0 + r(\ell - \ell_0))\right]\right\}\vert_{r=0} = 0,
\end{align*}
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 \citet{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 \textit{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$:\begin{align}\label{eq:plm_estimator}
        \hat{\theta}_n = \frac{\frac{1}{n}\sum_{i=1}^n \big(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\big)\big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)}{\frac{1}{n}\sum_{i=i}^n \big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)^2},
\end{align}
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 \citet{Chernozhukov2018}. \citet[Remark~3.1]{Chernozhukov2018} recommend `DML2' over the alternative `DML1' algorithm, which fits the final estimator by fold.}

\noindent\begin{minipage}{\linewidth}
\begin{sttech}[Algorithm 1. DDML for the Partially Linear Model.]
Split the sample $\{(Y_i, D_i, \bm{X}_i)\}_{i=1}^n$ randomly in $K$ folds of approximately equal size. Denote $I_k$ the set of observations included in fold $k$ and $I_k^c$ its complement.
\begin{enumerate}[nosep]
    \item For each $k\in\{1,\ldots, K\}$:
    \begin{enumerate}
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $Y_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{\ell}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $D_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{m}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
    \end{enumerate}
    \item Compute \eqref{eq:plm_estimator}.
\end{enumerate}
\end{sttech}
\end{minipage}

\citet{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)$. \texttt{ddml} computes the DDML estimator for the Partially Linear Model using Stata's \texttt{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 \texttt{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. \citet[][Remark~3.1]{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 \stopt{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 \stopt{fcluster}{varname}.

\subsection{The Interactive Model (\texttt{interactive})}\label{sec:interactive}

The Interactive Model is given by
\begin{align}\label{eq:interactive}
\begin{aligned}
    Y &= g_0(D,\bm{X}) + U
\end{aligned}
\end{align}
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\begin{align*}
    \theta^{\textrm{ATE}}_0 &\equiv E[g_0(1, \bm{X}) - g_0(0,\bm{X})]\\
    \theta^{\textrm{ATET}}_0 & \equiv E[g_0(1, \bm{X}) - g_0(0,\bm{X})\vert D = 1],
\end{align*}
which correspond to the average treatment effect (ATE) and average treatment effect on the treated (ATET), respectively.

Assumptions \ref{assumption:interactive_cmi} and \ref{assumption:interactive_overlap} 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.
\begin{assumption}[Conditional Mean Independence]\label{assumption:interactive_cmi}
$E[U\vert D, \bm{X}]=0$.
\end{assumption}
\begin{assumption}[Overlap]\label{assumption:interactive_overlap}
$\Pr(D=1 \vert \bm{X})\in (0,1)$ with probability 1.
\end{assumption}

Under assumptions \ref{assumption:interactive_cmi} and \ref{assumption:interactive_overlap}, we have\begin{align*}
    E[Y\vert D, \bm{X}] = E[g_0(D,\bm{X})\vert D, \bm{X}] + E[U\vert D, \bm{X}] = g_0(D,\bm{X}),
\end{align*}
so that identification of the ATE and ATET immediately follows from their definition.\footnote{In the defined Interactive Model under Assumption \ref{assumption:interactive_cmi}, 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 \eqref{eq:all_causes_1} that average over both observed and unobserved heterogeneity. See, e.g., \cite{Belloni2017}.}

In contrast to Section \ref{sec:plm}, 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 \citet{Chernozhukov2018} consider is the efficient influence function of \citet{Hahn1998}
\begin{align*}
    \psi^{\textrm{ATE}}(\bm{W};\theta, g, m) =  \frac{D(Y-g(1, \bm{X}))}{m(\bm{X})} {-} \frac{(1-D)(Y - g(0, \bm{X}))}{1 - m(\bm{X})} {+} g(1, \bm{X}) {-} g(0, \bm{X}) {-} \theta,
\end{align*}
where $\bm{W}\equiv (Y, D, \bm{X})$. Similarly for the ATET,
\begin{align*}
    \psi^{\textrm{ATET}}(\bm{W};\theta, g, m, p) =  \frac{D(Y - g(0,\bm{X}))}{p} - \frac{m(\bm{X})(1-D)(Y-g(0,\bm{X}))}{p(1-m(\bm{X}))} - \theta.
\end{align*}
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{assumption:interactive_cmi} and \ref{assumption:interactive_overlap} imply
\begin{align*}E[\psi^{\textrm{ATE}}(\bm{W};\theta_0^{\textrm{ATE}}, g_0, m_0)] &= 0 \\  E[\psi^{\textrm{ATET}}(\bm{W};\theta_0^{\textrm{ATET}}, g_0, m_0, p_0)] &= 0;\end{align*}
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
\begin{align}
\begin{split}\label{eq:ddml_ate}
   \hat{\theta}^{\textrm{ATE}}_n &=\frac{1}{n}\sum_{i=1}^n\Bigg(\frac{D_i(Y_i-\hat{g}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{m}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-D_i)(Y_i-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{m}_{I^c_{k_i}}(\bm{X}_i)} \\ &\hspace{.5\linewidth} +\hat{g}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i)\Bigg),
   \end{split}
   \\
   \begin{split}
      \hat{\theta}_n^{\textrm{ATET}} &=\frac{1}{n}\sum_{i=1}^n\left(\frac{D_i(Y_i-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i))}{\hat{p}} - \frac{\hat{g}_{I^c_{k_i}}(0,\bm{X}_i)(1-D_i)(Y_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i))}{\hat{p}(1-\hat{m}_{I^c_{k_i}}(\bm{X}_i))}\right)  ,\label{eq:ddml_atet}
      \end{split}
\end{align}
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{sec:plm}. 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} \equiv \frac{1}{n}\sum_{i=1}^nD_i$ is the sample share of treated observations.

\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 Appendix \ref{app:algorithms}. Mean and median aggregation over cross-fitting repetitions are implemented as outlined in Remark~2.

\section{DDML with Instrumental Variables}\label{sec:ddml_iv}
This section outlines the Partially Linear IV Model, the Flexible Partially Linear IV Model, and the Interactive IV Model. As in the previous section, each model is a special case of the general causal model \eqref{eq:all_causes_1}. 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 \eqref{eq:plm}, while the Interactive IV Model allows for arbitrary interactions between the treatment $D$ and the controls $\bm{X}$ as in \eqref{eq:interactive}. 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 \citet{Belloni2012} and \citet{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})$.

\subsection{Partially Linear IV Model (\texttt{iv})}\label{sec:pliv}
The Partially Linear IV Model considers the same functional form restriction on the causal model as the Partially Linear Model in Section \ref{sec:plm}. Specifically, the Partially Linear IV Model maintains
\begin{align*}
\begin{aligned}
        Y&= \theta_0 D  + g_0(\bm{X}) + U, \\
\end{aligned}
\end{align*}
where $\theta_0$ is the unknown parameter of interest.\footnote{As in Section \ref{sec:plm}, the interpretation of $\theta_0$ can be generalized under stronger identifying assumptions. See \citet{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{assumption:pliv_orthogonal} and \ref{assumption:pliv_relevance} below are sufficient orthogonality and relevance conditions, respectively, for identification of $\theta_0$.
\begin{assumption}[Conditional IV Orthogonality]\label{assumption:pliv_orthogonal}
$E[Cov(U,Z\vert \bm{X})]= 0$.
\end{assumption}
\begin{assumption}[Conditional Linear IV Relevance]\label{assumption:pliv_relevance}
$E[Cov(D,Z\vert \bm{X})]\neq 0$.
\end{assumption}

To show identification, consider the score function \begin{align*}
    \psi(\bm{W}; \theta, \ell,m,r) = \Big(Y-\ell(\bm{X}) - \theta(D - m(\bm{X}))\Big)\Big(Z - r(\bm{X})\Big),
\end{align*}
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{assumption:pliv_orthogonal} 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 \begin{align}\label{eq:pliv_alpha0}
    \theta_0 = \frac{E\left[\big(Y-\ell_0(\bm{X})\big)\big(Z - r_0(\bm{X})\big)\right]}{E\left[\big(D - m_0(\bm{X})\big)\big(Z - r_0(\bm{X})\big)\right]},
\end{align}
where Assumption \ref{assumption:pliv_relevance} is used to ensure a non-zero denominator.

The DDML estimator based on Equation \eqref{eq:pliv_alpha0} is given by \begin{align}\label{eq:pliv_ddml_estimator}
    \hat{\theta}_{n} = \frac{\frac{1}{n}\sum_{i=1}^n \left(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\right)\left(Z_i-\hat{r}_{I^c_{k_i}}(\bm{X}_i)\right)}{\frac{1}{n}\sum_{i=i}^n \left(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\right)\left(Z_i-\hat{r}_{I^c_{k_i}}(\bm{X}_i)\right)},
\end{align}
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.

\subsection{Flexible Partially Linear IV Model (\texttt{fiv})}\label{sec:hdiv}
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 \citet{Belloni2012} and \citet{Chernozhukov2015a}.
In particular, the Flexible Partially Linear IV Model leverages a conditional mean independence assumption rather than an orthogonality assumption as in Section \ref{sec:pliv}. As in Section \ref{sec:pliv}, we state everything in the case of a scalar $D$.

\begin{assumption}[Conditional IV Mean Independence]\label{assumption:hdiv_iv_meanindependence}
$E[U\vert \bm{Z},\bm{X}] = 0$.
\end{assumption}
Assumption \ref{assumption:hdiv_iv_meanindependence} implies that for any function $\tilde{p}(\bm{Z},\bm{X})$, it holds that  \begin{align}\label{eq:hdiv_moment}
    E\left[\Big(Y - \ell_0(X) - \theta\big(D - m_0(\bm{X})\big)\Big)\Big(\tilde{p}(\bm{Z},\bm{X}) - E\left[\tilde{p}(\bm{Z},\bm{X})\vert \bm{X}\right]\Big)\right] = 0,
\end{align}
where $\ell_0(\bm{X})=E[Y\vert \bm{X}]$ and $m_0(\bm{X})=E[D\vert \bm{X}]$. Identification based on \eqref{eq:hdiv_moment} requires that there exists some function $\tilde{p}$ such that \begin{align}\label{eq:hdiv_ptld_relevance}
    E\left[Cov(D, \tilde{p}(\bm{Z},\bm{X})\vert \bm{X})\right] \neq 0.
\end{align}
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 \eqref{eq:hdiv_ptld_relevance}.\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{assumption:hdiv_relevance} is a consequence of this non-mean independence.
\begin{assumption}[Conditional IV Relevance]\label{assumption:hdiv_relevance}
$E[Var(E[D\vert \bm{Z},\bm{X}]\vert \bm{X})] \neq 0.$
\end{assumption}

Consider now the score function \begin{align*}
    \psi(\bm{W}; \theta, \ell, m, p) = \Big(Y-\ell(\bm{X}) - \theta(D - m(\bm{X}))\Big)\Big(p(\bm{Z}, \bm{X}) - m(\bm{X})\Big),
\end{align*}
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{assumption:hdiv_iv_meanindependence} 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 \begin{align}\label{eq:hdiv_alpha0}
    \theta_0 = \frac{E\left[\big(Y-\ell_0(\bm{X})\big)\big(p_0(\bm{Z}, \bm{X}) - m_0(\bm{X})\big)\right]}{E[\big(D-m_0(\bm{X})\big)\big(p_0(\bm{Z}, \bm{X}) - m_0(\bm{X})\big)]},
\end{align}
where Assumption \ref{assumption:hdiv_relevance} ensures a non-zero denominator.

The DDML estimator based on the moment solution \eqref{eq:hdiv_alpha0} is given by \begin{align}\label{eq:hdiv_ddml_estimator}
        \hat{\theta}_{n} = \frac{\frac{1}{n}\sum_{i=1}^n \big(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\big)\big(\hat{p}_{I^c_{k_i}}(\bm{Z}_i, \bm{X}_i)-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)}{\frac{1}{n}\sum_{i=i}^n \big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)\big(\hat{p}_{I^c_{k_i}}(\bm{Z}_i, \bm{X}_i)-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)},
\end{align}
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 \eqref{eq:hdiv_ddml_estimator} 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 \eqref{eq:hdiv_ddml_estimator}.

\noindent\begin{minipage}{\linewidth}
\begin{sttech}[Algorithm 2. LIE-compliant DDML for the Flexible Partially Linear IV Model.]
 Split the sample $\{(Y_i, D_i, \bm{X}_i, \bm{Z}_i)\}_{i=1}^n$ randomly in $K$ folds of approximately equal size. Denote $I_k$ the set of observations included in fold $k$ and $I_k^c$ its complement.
\begin{enumerate}[nosep]
    \item For each $k\in\{1,\ldots, K\}$:
    \begin{enumerate}
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $Y_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{\ell}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $D_i$ as the outcome and $(\bm{Z}_i,\bm{X}_i)$ as predictors. Obtain the out-of-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ for $i\in I_k$ and in-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ for $i\in I^c_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using the in-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ as the outcome and $X_i$ as predictors. Obtain the out-of-sample predicted values $\hat{m}_{I^c_k}(\bm{X}_i)$ for $i\in I_k.$
    \end{enumerate}
    \item Compute \eqref{eq:hdiv_ddml_estimator}.
\end{enumerate}
\end{sttech}
\end{minipage}

Standard errors corresponding to $\hat{\theta}_{n}$ in  \eqref{eq:hdiv_ddml_estimator} are the same as in Section \ref{sec:pliv} 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.

\subsection{Interactive IV Model (\texttt{interactiveiv})}\label{sec:interactiveiv}

The Interactive IV Model considers the same causal model as in Section \ref{sec:interactive}; specifically
\begin{align*}
\begin{aligned}
    Y &= g_0(D, \bm{X}) + U \\
\end{aligned}
\end{align*}
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$.

The parameter of interest we target is\begin{align}\label{interactiveiv:fuzzy_late}
\begin{aligned}
   \theta_0 & =  E\left[g_0(1, \bm{X}) - g_0(0, \bm{X})\vert \: p_0(1,\bm{X}) > p_0(0,\bm{X})\right],
\end{aligned}
\end{align}
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 \citet{imbens1994identification}, ``local'' here 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. Under these stronger assumptions, the DDML LATE estimator outlined here targets the conventionally considered LATE parameter.}

Identification again leverages Assumptions \ref{assumption:hdiv_iv_meanindependence} and \ref{assumption:hdiv_relevance} 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.
\begin{assumption}[Monotonicity]\label{assumption:iv_monotonicty}
$p_0(1, \bm{X})\geq p_0(0, \bm{X})$ with probability 1.
\end{assumption}
\begin{assumption}[IV Overlap]\label{assumption:iv_overlap}
$\Pr(Z=1\vert \bm{X})\in(0,1)$ with probability 1.
\end{assumption}

Assumptions \ref{assumption:hdiv_iv_meanindependence}-\ref{assumption:iv_overlap} imply that \begin{align}\label{eq:interactiveiv_flate}
    \theta_0 = \frac{E\left[\ell_0(1, \bm{X})-\ell_0(0, \bm{X})\right]}{E\left[p_0(1, \bm{X})-p_0(0, \bm{X})\right]},
\end{align}
where $\ell_0(Z, \bm{X}) \equiv E[Y\vert Z, \bm{X}]$, verifying identification of the LATE $\theta_0$. Akin to Section \ref{eq:interactive}, however, estimators of $\theta_0$ should not directly be based on Equation \eqref{eq:interactiveiv_flate} because the estimating equations implicit in obtaining \eqref{eq:interactiveiv_flate} do not satisfy Neyman-orthogonality. Hence, a direct estimator of $\theta_0$ obtained by plugging nonparametric estimators in for nuisance functions in  \eqref{eq:interactiveiv_flate} will potentially be highly sensitive to the first step nonparametric estimation error. Rather, we base estimation on the Neyman orthogonal score function \begin{align*}
    \psi(\bm{W};\theta, \ell, p, r) &= \frac{Z(Y - \ell(1, \bm{X}))}{r(\bm{X})} - \frac{(1-Z)(Y-\ell(0, \bm{X}))}{1-r(\bm{X})} + \ell(1, \bm{X})-\ell(0, \bm{X}) \\
    &  + \left[\frac{Z(D - p(1,\bm{X}))}{r(\bm{X})} - \frac{(1-Z)(D-p(0,\bm{X}))}{1-r(\bm{X})} + p(1,\bm{X})-p(0,\bm{X})\right] \times \theta
\end{align*}
where $\bm{W} \equiv (Y, D, \bm{X}, Z)$ . Note that under Assumptions \ref{assumption:hdiv_iv_meanindependence}-\ref{assumption:iv_overlap} 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 \begin{align}
\begin{split}
    \hat{\theta}_n &= \\
    &\frac{\frac{1}{n}\sum_{i}\left(\frac{Z_i(Y_i-\hat{\ell}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{r}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-Z_i)(Y_i-\hat{\ell}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{r}_{I^c_{k_i}}(\bm{X}_i)}+\hat{\ell}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{\ell}_{I^c_{k_i}}(0,\bm{X}_i)\right)}{\frac{1}{n}\sum_{i}\left(\frac{Z_i(D_i-\hat{p}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{r}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-Z_i)(D_i-\hat{p}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{r}_{I^c_{k_i}}(\bm{X}_i)}+\hat{p}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{p}_{I^c_{k_i}}(0,\bm{X}_i)\right)}, \label{eq:ddml_late}
\end{split}
\end{align}
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 only using treated and untreated 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.

\section{The choice of machine learner}\label{sec:stacking}
\citet{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 \citep{bickel2009simultaneous, Belloni2012}, random forests \citep{wager2015adaptive, wager2018, Athey2019a}, neural networks \citep{schmidt2020nonparametric, farrell2021deep}, and boosting \citep{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 \emph{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 \citet{Wolpert1992} and \citet{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.

We introduce stacking for DDML in Section~\ref{eq:stacking}. Section~\ref{eq:stacking_mc} demonstrates the performance of DDML in combination with stacking using a simulation.

\subsection{Stacking}\label{eq:stacking}
Our discussion of stacking in the context of DDML focuses on the Partially Linear Model in \eqref{eq:plm}, 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:
\[  \underset{w_{k,1},\ldots,w_{k,J}}{\min} \sum_{i\in T_k} \left( Y_i - \sum_{j=1}^J w_{k,j} \hat{\ell}^{(j)}_{T_{k,v(i)}^c}(\bm{X}_i)  \right)^2,   \qquad \textrm{s.t.}\ w_{k,j}\geq 0,\ \sum_{j=1}^J |w_{k,j}|=1,\]
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)$ which observation $i$ falls into. 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 \citep{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 \texttt{pystacked} \citep{Ahrens2022}.

\paragraph{Short-stacking.} Stacking relies 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 \emph{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 is the lower computational cost. Short-stacking also produces a single set of weights for the entire sample, rather than a different set of weights in each cross-fit fold and thus facilitates interpretation. 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.}

\subsection{Monte Carlo simulation}\label{eq:stacking_mc}
To illustrate the advantages of DDML with stacking, we generate artificial data based on the Partially Linear Model
\begin{align}
Y_i &= \theta_0D_i + c_Yg(\bm{X}_i) + \sigma_Y(D_i,\bm{X}_i)\varepsilon_i \label{eq:mc1}\\
D_i &= c_Dg(\bm{X}_i) + \sigma_D(\bm{X}_i)u_i \label{eq:mc2}
\end{align}
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 \eqref{eq:mc1} and \eqref{eq:mc2}
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:
\[
\begin{array}{ll}
\textrm{DGP 1:}\quad g(\bm{X}_i) &= \displaystyle\sum\nolimits_j 0.9^jX_{ij} \\
\textrm{DGP 2:}\quad g(\bm{X}_i) &= X_{i1}X_{i2}{+}X_{i3}^2{+}X_{i4}X_{i5}{+}X_{i6}X_{i7}{+}X_{i8}X_{i9}{+}X_{i10}{+} X_{i11}^2{+}X_{i12} X_{i13} \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 3:}\quad g(\bm{X}_i) &= \mathbbm{1}\{X_{i1}>0.3\}\mathbbm{1}\{X_{i2}>0\}\mathbbm{1}\{X_{i3}>-1\}  \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 4:}\quad g(\bm{X}_i) &= X_{i1} + \sqrt{|X_{i2}|} + \textrm{sin}(X_{i3}) + 0.3X_{i4}X_{i5} + X_{i6} + 0.3X_{i7}^2  \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 5:}\quad g(\bm{X}_i)&=\textnormal{same as DGP 4 with } p=7 \phantom{\displaystyle\sum\nolimits_j} \\
\end{array}\]
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 which 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 \texttt{pystacked}. We use the defaults of \texttt{pystacked} for parameter values and settings not mentioned here.}
\begin{enumerate}[noitemsep,topsep=0pt,itemsep=0pt,parsep=0pt,before=\vspace{0mm},after=\vspace{-0mm}]
    \item[1.-2.] Cross-validated lasso \& ridge with untransformed base controls
    \item[3.-4.] Cross-validated lasso \& ridge with 5th-order polynomials of base controls but no interactions (referred to as `Poly 5')
    \item[5.-6.] Cross-validated lasso \& ridge with second-order polynomials and all first order interaction terms (referred to as `Poly 2 + Inter.')
    \item[7.] Random forests (RF) with low regularization: base controls, maximum tree depth of 10, 500 trees and approximately $\sqrt{p}$ features considered at each split
    \item[8.] RF with medium regularization: same as 7., but with maximum tree depth of 6
    \item[9.] RF with high regularization: same as 7., but with maximum tree depth of 2
    \item[10.] Gradient boosted trees (GB) with low regularization: base controls, 1000 trees and a learning rate of 0.3. We enable early stopping which uses a 20\% validation sample to decide whether to stop the learning algorithm. Learning is terminated after 5 iterations with no meaningful improvement in the mean-squared loss of the validation sample.\footnote{We use a tolerance level of 0.01 to measure improvements.}
    \item[11.] GB with medium regularization: same as 10., but with learning rate of 0.1
    \item[12.] GB with high regularization: same as 10., but with learning rate of 0.01
    \item[13.] Feed-forward neural net with base controls and two layers of size 20
\end{enumerate}

\noindent We use the above set of learners as base learners for DDML with stacking approaches. Specifically, we estimate DDML using stacking regression with CLS, stacking with single-best learner, and short-stacking with CLS. We set the number of folds to $K=20$ if $n=100$, and $K=5$ if $n=1000$. For comparison, we report results for OLS using the base controls, PDS-Lasso with the 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{tab:mc1}. We delegate results for DGPs~4 and~5, including a brief discussion, to  Appendix~\ref{appendix:moremc}. DDML estimators leveraging stacking or short-stacking 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 imprecisely estimated in very 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 four stacking approaches considered. 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~\ref{tab:mc2}). 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~\ref{tab:stacking_weights_appendix} and \ref{tab:mspe_appendix} 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 compare 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.

\begin{sidewaystable}
    \centering\scriptsize\singlespacing
\begin{threeparttable}
    \caption{Bias and Coverage Rates in the Linear and Nonlinear DGPs}
    \label{tab:mc1}
        \begin{tabular}{rlcccccccccccccc}
    \toprule
    \midrule
          &       & \multicolumn{4}{c}{DGP 1} & & \multicolumn{4}{c}{DGP 2} & & \multicolumn{4}{c}{DGP 3} \\
    \cline{3-6} \cline{8-11} \cline{13-16} \\[-6pt]
     & & \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}    &&  \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}     &&  \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}       \\
    && MAB & CR.   &  MAB & CR   &&   MAB & CR   &  MAB & CR   &&  MAB & CR   & MAB & CR   \\ \midrule
          \multicolumn{2}{l}{Full sample:}\\
          \partialinput{1}{5}{Simul/sim_SJ/bias_1.tex}
           \multicolumn{2}{l}{DDML methods:}\\
           \multicolumn{2}{l}{\it ~~Base learners}\\
          \partialinput{5}{19}{Simul/sim_SJ/bias_1.tex}
          \multicolumn{2}{l}{\it ~~Meta learners}\\
          \partialinput{20}{23}{Simul/sim_SJ/bias_1.tex}
    \midrule
    \bottomrule
    \end{tabular}%
  \begin{tablenotes}[para,flushleft]
  \footnotesize
  \item \textit{Notes:} The table reports median absolute bias (MAB) and coverage rate of a  95\% confidence interval (CR). We employ standard errors robust to heteroskedasticity.
  For comparison, we report the following full sample estimators: infeasible Oracle, OLS, PDS-Lasso with base and two different expanded sets of covariates.
  DDML estimators use 20 folds for cross-fitting if $n=100$, and 5 folds if $n=1000$. Meta-learning approaches rely on all listed base learners.
  Results are based on 1,000 replications. Results for DGPs 4 and 5 can be found in Table~\ref{tab:mc2} in the Appendix.
  \end{tablenotes}
\end{threeparttable}
\end{sidewaystable}

\section{The program}\label{sec: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{sec:syntax_ddml}. Section~\ref{sec:options} lists the options. Section~\ref{sec:syntax_qddml} covers the simplified one-line program \stcmd{qddml}\xspace. We provide an overview of supported machine learning programs in Section~\ref{sec:supported_programs}. Finally, Section~\ref{sec:replication} adds a note on how to ensure replication with \stcmd{ddml}\xspace.

\subsection{Syntax: \stcmd{ddml}\xspace}\label{sec:syntax_ddml}
The \texttt{ddml} estimation proceeds in four steps.

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

\begin{stsyntax}
ddml init {\it model}
    \optional{,
    mname({\it name})
    vars({\it varlist})
    kfolds({\it integer})
    fcluster({\it varname})
    foldvar({\it varlist})
    reps({\it integer})
    tabfold
    vars({\it varlist})
    }
\end{stsyntax}
\noindent where {\it model} selects between the Partially Linear Model (\texttt{partial}), the Interactive Model (\texttt{interactive}), the Partially Linear IV Model (\texttt{iv}), the Flexible Partially Linear IV Model (\texttt{fiv}), and the Interactive IV Model (\texttt{interactiveiv}). This step creates a persistent Mata object with the name provided by \stopt{mname}{name} 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})}).

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

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

\begin{stsyntax}
ddml {\it cond\_exp}
    \optional{,
    mname({\it name})
    vname({\it varname})
    learner({\it name})
    vtype({\it string})
    predopt({\it string})
    }
    :
    {\it command} {\it depvar} {\it vars}
    \optional{,
    {\it cmdopt}
    }
\end{stsyntax}

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

\begin{table}[H]
\newcommand{\xmark}{}
    \centering
    \begin{tabular}{l|ccccc}\hline\hline
       {\it cond\_exp}      & \texttt{partial} &   \texttt{interactive} & \texttt{iv} & \texttt{fiv} &  \texttt{late} \\ \hline
       \verb+E[Y|X]+ &  \checkmark &  \xmark & \checkmark & \checkmark & \xmark\\
       \verb+E[Y|X,D]+ &\xmark & \checkmark &\xmark & \xmark & \xmark \\
       \verb+E[Y|X,Z]+ & \xmark & \xmark & \xmark & \xmark & \checkmark\\[.1cm]
       \verb+E[D|X]+ & \checkmark & \checkmark & \checkmark & \checkmark & \xmark \\
       \verb+E[D|Z,X]+ & \xmark & \xmark & \xmark & \checkmark & \checkmark  \\[.1cm]
       \verb+E[Z|X]+  & \xmark & \xmark & \checkmark & \xmark & \checkmark \\\hline\hline
    \end{tabular}
    \caption{The table lists the conditional expectations which need to be specified for each model.}
    \label{tab:cond_exp}
\end{table}

\subsubsection*{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.

\begin{stsyntax}
ddml crossfit
    \optional{,
    mname({\it name})
    {\ttfamily shortstack}
    }
\end{stsyntax}

\subsubsection*{Step 4: Estimate causal effects.}

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

\begin{stsyntax}
ddml estimate
    \optional{,
    mname({\it name})
    robust
    cluster({\it varname})
    vce({\it vcetype})
    att
    trim
    spec({\it string})
    rep({\it string})
    }
\end{stsyntax}

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

\begin{stsyntax}
ddml estimate
    \optional{,
    replay
    mname({\it name})
    spec({\it integer or string})
    rep({\it integer or string})
    fulltable
    notable
    allest
    }
\end{stsyntax}

\subsection{Options}\label{sec:options}

\paragraph{Step 1 options: Initialization.}
\begin{description}
    \item[\stopt{mname}{{\it name}}] name of the DDML model. Allows to run multiple DDML models simultaneously. Defaults to {\it m0}.
    \item[\stopt{kfolds}{{\it integer}}] number of cross-fitting folds. The default is 5.
    \item[\stopt{fcluster}{{\it varname}}] cluster identifiers for cluster randomization of folds.
    \item[\stopt{foldvar}{{\it varlist}}] integer variables to specify custom folds (one per cross-fitting repetition).
    \item[\stopt{reps}{{\it integer}}] number of cross-fitting repetitions, i.e., how often the cross-fitting procedure is repeated on randomly generated folds.
    \item[{{\tt tabfold}}]  prints a table with frequency of observations by fold.
\end{description}

\paragraph{Step 2 options: Adding learners.}
\begin{description}
    \item[\stopt{vname}{{\it varname}}] name of the dependent variable in the reduced form estimation.  This is usually inferred from the {\it command} line but is mandatory for the {\tt fiv} model.
    \item[\stopt{learner}{{\it varname}}] optional name of the variable to be created.
    \item[\stopt{vtype}{{\it string}}] optional variable type of the variable to be created. Defaults to {\it double}.  {\it none} can be used to leave the type field blank. (Setting \stopt{vtype}{{\it none}} is required when using \stcmd{ddml}\xspace with {\tt rforest}.)
    \item[\stopt{predopt}{{\it varname}}]  {\tt predict} option to be used to get predicted values.  Typical values could be {\tt xb} or {\tt pr}. Default is blank.
\end{description}

\paragraph{Step 3 options: Cross-fitting.}
\begin{description}
    \item[{\ttfamily shortstack}]  asks for short-stacking to be used.  Short-stacking runs constrained least squares on the cross-fitted predicted values to obtain a weighted average of several base learners.
\end{description}

\paragraph{Step 4 options: Estimation.}
\begin{description}
    \item[\stopt{spec}{{\it string}}] select specification. This can either be the specification number, {\it mse} for minimum-MSE specification (the default) or {\it ss} for short-stacking.
    \item[\stopt{rep}{{\it string}}] select cross-fitting repetitions. This can either be the cross-fit repetition number, {\it mn} for mean aggregation or {\it md} for median aggregation (the default). See Remark~2 for more information.
    \item[{\tt robust}] report SEs that are robust to the presence of arbitrary heteroskedasticity.
    \item[\stopt{cluster}{{\it varname}}] select cluster-robust variance-covariance estimator.
    \item[\stopt{vce}{{\it type}}] select variance-covariance estimator, e.g.\ \texttt{vce(hc3)} or \texttt{vce(cluster id)}. See \verb+help regress##vcetype+ for available options.
    \item[\stopt{trim}{{\it real}}] trimming of propensity scores for the Interactive and Interactive IV models. The default is 0.01 (that is, values below 0.01 and above 0.99 are set to 0.01 and 0.99, respectively).
    \item[{\tt atet}] report average treatment effect of the treated (default is ATE).
    \item[{\tt noconstant}] suppress constant term in the estimation stage (only relevant for partially linear models).
\end{description}

\subsection{Short syntax: \stcmd{qddml}\xspace}\label{sec:syntax_qddml}
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 \texttt{pdslasso} and \texttt{ivlasso} \citep{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 \texttt{regress}, which is added by default).

\subsubsection*{Syntax for Partially Linear and Interactive Model}
\begin{stsyntax}
qddml
    \depvar\
    {\it treatment\_vars}
    \textnormal{(}{\it controls}\textnormal{)}\verb+,+
    \verb+model(partial|interactive)+
    \optional{
    {\it options}
    }
\end{stsyntax}

\subsubsection*{Syntax for IV models}
\begin{stsyntax}
qddml
    \depvar\
    \textnormal{(}{\it controls}\textnormal{)}
    \textnormal{(}{\it treatment\_vars}=excluded\_instruments\textnormal{)}
    \verb+,+
    \verb+model(iv|late|fiv)+
    \optional{
    cmd({\it string})
    cmdopt({\it string})
    {\it ddml\_options}
    {\tt noreg}
    }
\end{stsyntax}

\noindent where {\it ddml\_options} options are listed in Section~\ref{sec:options} and are internally passed on to the \stcmd{ddml}\xspace sub-routines. \stopt{cmd}{string} selects the machine learning command to be used and \stopt{cmdopt}{string} allows to pass options to the estimation command. \stcmd{qddml}\xspace adds by default \stcmd{regress} to the user-specified learner; this behavior can be disabled with \stcmd{noreg}. See the \stcmd{qddml}\xspace help file for a full list of options.

\subsection{Supported machine learning programs}\label{sec:supported_programs}
\stcmd{ddml}\xspace is compatible with any supervised ML program in Stata that supports the typical ``\texttt{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:
\begin{itemize}[nosep]
    \item \texttt{lassopack} implements regularized regression, e.g.\ lasso, ridge, elastic net \citep{Ahrens2019}.
    \item \stcmd{pystacked}\xspace facilitates the stacking of a wide range of machine learners including regularized regression, random forests, support vector machines, gradient boosted trees and feed-forward neural nets using Python's \emph{scikit-learn} \citep{Ahrens2022,scikit-learn,sklearn_api}. In addition, \stcmd{pystacked}\xspace can also be used as a front-end to fit individual machine learners.

    \item \texttt{rforest} is a random forest wrapper for WEKA \citep{Schonlau2020,frank2009weka}.
    \item \texttt{svmachines} allows for the estimation of support vector machines using \texttt{libsvm} \citep{chang2011libsvm,Guenther2018}.
    \item The program \texttt{parsnip} of the package \texttt{mlrtime} provides access to R's \emph{parsnip} machine learning library through \texttt{rcall} \citep{Huntington2021,haghish2019seamless}. Using \texttt{parsnip} requires the installation of the supplementary wrapper program \texttt{parsnip2}.\footnote{Available from {\tt\url{https://github.com/aahrens1/parsnip2}.}}
\end{itemize}

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

\subsection{Inspecting results and replication}\label{sec: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_+\textit{r}, where \verb+m0+ is the default model name and \textit{r} is the cross-fitting repetition). \stcmd{ddml}\xspace models are stored on Mata \texttt{struct}s and using Mata's associative arrays. Specifically, the \stcmd{ddml}\xspace model created by \stcmd{ddml}\xspace \texttt{init} is an \texttt{mStruct}, and information relating to the estimation of conditional expectations are stored in \texttt{eStruct}s. 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:
\begin{itemize}
    \item Set the Stata seed before {\tt ddml init}. This ensure that the same random fold variable is used for a given data set.
    \item Using the same fold variable alone is usually not sufficient to ensure replication, since many machine learning algorithms involve randomization. That said, note that the fold variable is stored in memory and can be reused for subsequent estimations via the \stopt{foldvar}{varlist} option.
    \item Replication of \stcmd{ddml}\xspace results may require additional steps with some programs that rely on randomization in other software environments, e.g., R or Python. \stcmd{pystacked}\xspace uses a Python seed generated in Stata. Thus, when \stcmd{ddml}\xspace is used with \stcmd{pystacked}\xspace, setting the seed before {\tt ddml init} also guarantees that the same Python seed underlies the stacking estimation. Other programs relying on randomization outside of Stata might not behave in the same way. Thus, when using other programs, check the help files for options to set external random seeds. Try estimating each of the individual learners on the entire sample to see what settings need to be passed to them for their results to replicate.

    \item The \stcmd{ddml}\xspace \texttt{extract} utility can be used to retrieve and inspect a wide range of intermediate results and statistics from Mata structures. This can be useful for trying to track down estimation errors and replication anomalies.
    \item Beware of changing samples. Fold splits or learner idiosyncracies may mean that sample sizes vary slightly across learners, estimation samples and/or cross-fitting repetitions. Note that \stcmd{ddml}\xspace stores sample indicators for cross-fitting repetitions as Stata variables; \stcmd{ddml}\xspace \texttt{extract} with the \texttt{show(n)} option will report sample sizes by learner and fold. See the {\tt ddml extract} help file for more information.
    \item The \stcmd{ddml}\xspace \texttt{export} utility can be used to export the estimated conditional expectations, fold variables and sample indicators to a CSV format file for examination and comparison in other software environments.
\end{itemize}

\section{Applications}\label{sec:application}
We demonstrate the \texttt{ddml} workflow using two applications. In Section~\ref{sec:401k}, we apply the DDML estimator to estimate the effect of 401(k) eligibility on financial wealth following \citet{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~\ref{appendix:moreapplications}. Additional examples can also be found in the help file.
Based on \citet{berry1995automobile}, we show in Section~\ref{sec:ivhd_application} 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.

\subsection{401(k) and financial wealth}\label{sec:401k}
The data consists of $n=9915$ households from the 1991 SIPP. The application is originally due to \citet{poterba1995}, but has been revisited by \citet{Belloni2017}, \citet{Chernozhukov2018}, and \citet{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{sec:syntax_ddml}.

\begin{stlog}
\input{Applications/logs/partial_start.log.tex}
\end{stlog}

\subsubsection*{Step 1: Initialize \texttt{ddml} model.} We initialize the \stcmd{ddml}\xspace model and select the Partially Linear Model in \eqref{eq:plm}. Before initialization, we set the seed to ensure replication. This should always be done before \texttt{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.

\begin{stlog}
\input{Applications/logs/partial_init.log.tex}
\end{stlog}

\subsubsection*{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 \texttt{regress} for linear regression, \texttt{cvlasso} for cross-validated lasso and \texttt{rforest} for random forests. When using \texttt{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.

\begin{stlog}
\input{Applications/logs/partial_add.log.tex}
\end{stlog}

The flexible \stcmd{ddml}\xspace syntax allows specifying 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, for example, 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 \texttt{cvlasso}, 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 \texttt{ddml E[D|X]: reg D2 \$X} where \texttt{D2} would be the name of the second treatment variable. An example with two treatments is provided in the help file.

The auxiliary \texttt{ddml} sub-command \texttt{describe} allows to verify that the learners were correctly registered:

\begin{stlog}
\input{Applications/logs/partial_describe.log.tex}
\end{stlog}

\subsubsection*{Step 3: Perform cross-fitting.} The third step performs cross-fitting, which is the most time-intensive process. The \texttt{shortstack} option enables the short-stacking algorithm of Section~\ref{eq:stacking}.

\begin{stlog}
\input{Applications/logs/partial_crossfit.log.tex}
\end{stlog}

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_cvlasso_1+, \verb+Y3_rforest_1+, \verb+D1_reg_1+, \verb+D2_cvlasso_1+ and \verb+D3_rforest_1+. \texttt{Y} and \texttt{D} indicate outcome and treatment variable. The index \texttt{1} to \texttt{3} is a learner counter. \texttt{reg}, \texttt{cvlasso} and \texttt{rforest} correspond to the name of the commands used. The \verb+_1+ suffix  indicates the cross-fitting repetition.

After cross-fitting, we can inspect the mean-squared prediction errors by fold and learner:
\begin{stlog}
\input{Applications/logs/partial_extract.log.tex}
\end{stlog}

\subsubsection*{Step 4: Estimate causal effects.}

In this final step, we obtain the causal effect estimates. Since we requested short-stacking in Step~3, \texttt{ddml} shows 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 (denoted as \texttt{opt}).

\begin{stlog}
\input{Applications/logs/partial_estimate.log.tex}
\end{stlog}

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_cvlasso_1+ and \verb+Y3_rforest_1+ with \verb+D1_reg_1+, \verb+D2_cvlasso_1+ and \verb+D3_rforest_1+). To get all results, we add the \texttt{allcombos} option:

\begin{stlog}
\input{Applications/logs/partial_estimate_all.log.tex}
\end{stlog}

We can use the \stopt{spec}{string} option to select among the listed specifications. {\it string} is either the specification number, \texttt{ss} to get the short-stacking specification 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 \stopt{spec}{string} option can be provided either in combinations with \texttt{allcombos}, or after estimation in combination with the {\tt replay} option, for example:

\begin{stlog}
\input{Applications/logs/partial_spec.log.tex}
\end{stlog}

\subsubsection*{Manual final estimation.}
In the background, \texttt{ddml estimate} regresses \verb+Y1_reg_1+ against \verb+D1_reg_1+ with a constant. We can verify this manually:

\begin{stlog}
\input{Applications/logs/partial_reg.log.tex}
\end{stlog}

Manual estimation using {\ttfamily regress} allows the use of postestimation regression tools such as {\ttfamily avplot}.

\subsubsection*{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. For demonstration, we use \texttt{cvlasso}:

\begin{stlog}
\input{Applications/logs/partial_qddml.log.tex}
\end{stlog}

\noindent \stcmd{qddml}\xspace always includes linear regression (\stcmd{regress}) as an additional learner. The option \stcmd{noreg} disables this default behavior.

\subsubsection*{Using \stcmd{pystacked}\xspace with single learners}

The above example relying on \texttt{cvlasso} and \texttt{rforest} runs relatively slowly at 380 seconds. We can substantially improve the speed by using the lasso and random forests implementation offered by \stcmd{pystacked}\xspace, which calls the Python library \emph{scikit-learn}. This reduces the total computational time to only 29 seconds.

\begin{stlog}
\input{Applications/logs/partial_2_add.log.tex}
\end{stlog}

\noindent Note that, despite its name, \stcmd{pystacked}\xspace does not execute stacking regression in this example, since only one learner is specified in each call to \stcmd{pystacked}\xspace.

\subsubsection*{Stacking}

We next demonstrate DDML with stacking. To this end, we exploit the stacking regressor implemented in \stcmd{pystacked}\xspace. \stcmd{pystacked}\xspace allows to combine multiple base learners with learner-specific settings and covariates into a final meta learner. The learners are separated by \verb+||+. \stopt{method}{name} selects the learner, \stopt{xvars}{varlist} specifies learner-specific covariates (overwritting the default covariates \verb+$X+) and \stopt{opt}{string} 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 \citet{Ahrens2022}.

\begin{stlog}
\input{Applications/logs/partial_stack_add.log.tex}
\end{stlog}

\noindent After cross-fitting, we can obtain the MSE and stacking weights averaged over folds:
\begin{stlog}
\input{Applications/logs/partial_stack_crossfit.log.tex}
\end{stlog}

By adding the \texttt{detail} option, the user can display the stacking weights for each fold.
 Finally, we retrieve the results of DDML with stacking:
\begin{stlog}
\input{Applications/logs/partial_stack_estimate.log.tex}
\end{stlog}

\noindent The run time of this example is 83 seconds, which despite its computational complexity is faster than using \stcmd{ddml}\xspace with \texttt{cvlasso} and \texttt{rforest}.

\subsection{The market for automobiles}\label{sec:ivhd_application}
For this demonstration, we follow \citet{chernozhukov2015} who estimate a stylized demand model using instrumental variables based on the data from \citet{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, \citet{berry1995automobile} suggest exploiting other products' characteristics as instruments. Following \citet{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. \citet{chernozhukov2015} also construct an augmented set of instruments, including first-order interactions, squared and cubic terms. In the analysis below, we extend \citet{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:
\begin{stlog}
\input{Applications/logs/ivhd_start.log.tex}
\end{stlog}

\paragraph{Step 1: Initialize \texttt{ddml} model.}

\begin{stlog}
\input{Applications/logs/ivhd_init.log.tex}
\end{stlog}

Note that in the \texttt{ddml init} step, we include the option \texttt{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 \texttt{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 \stopt{learner}{string} option to match one learner for $E[D|\bm{X}]$ with a learner for $E[D|\bm{X},\bm{Z}]$, and \stopt{vname}{varname} to indicate the name of the treatment variable.

\begin{stlog}
\input{Applications/logs/ivhd_add.log.tex}
\end{stlog}

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

\begin{stlog}
\input{Applications/logs/ivhd_estimate.log.tex}
\end{stlog}

\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 ``\texttt{\_1}'' to ``\texttt{\_}$r$''.

\begin{stlog}
\input{Applications/logs/ivhd_ivreg.log.tex}
\end{stlog}

\paragraph{One-line syntax (output omitted).}

\begin{stlog}
\input{Applications/logs/ivhd_qddml.log.tex}
\end{stlog}

\section{Acknowledgments}
We thank users who tested earlier versions of the program. We also thank Jan Ditzen, Ben Jann, Eroll Kuhn and Di Liu for helpful comments, as well as participants at the German Stata Conference 2021, Italian Stata Conference 2022 and Swiss Stata Conference 2022. All remaining errors are our own.

\bibliographystyle{sj}

\section{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, \citet{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 \citet{roberts2020a} or digital trace data \citep{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 \citet{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 \citep{Belloni2012}. Lasso-based approaches for estimation of causal effects have become a popular strategy in applied econometrics \citep[e.g.][]{gilchrist2016something,dhar2022reshaping}, partially facilitated by the availability of software programs in Stata (\stcmd{pdslasso}\xspace, \citealp{Ahrens2018}; \citealp{statacorp2019}) and R (\texttt{hdm}, \citealp{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, \citet{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., \citealp{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, \citet{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 \citep{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 1.2 of \stcmd{ddml}\xspace.} \stcmd{ddml}\xspace adds to a small number of programs for causal machine learning in Stata (\citealp{Ahrens2018}, \citealp{statacorp2019}, \citealp{statacorp2021}). We briefly summarize the four main features of the program:
\begin{enumerate}

\item \stcmd{ddml}\xspace supports flexible estimators of causal parameters in five econometric models: (1) the Partially Linear Model, (2) the Interactive Model (for binary treatment), (3) the Partially Linear IV Model, (4) the Flexible Partially Linear IV Model, and (5) the Interactive IV Model (for binary treatment and instrument).

\item \stcmd{ddml}\xspace supports data-driven combinations of multiple machine learners via stacking by leveraging \texttt{pystacked} \citep{Ahrens2022}, our complementary Stata frontend relying on the Python library \emph{scikit-learn} \citep{scikit-learn,sklearn_api}.

\item Aside from \texttt{pystacked}, \stcmd{ddml}\xspace can be used in combination with many other existing supervised machine learning programs available in or via Stata. \stcmd{ddml}\xspace has been tested with \texttt{lassopack} \citep{Ahrens2019}, \texttt{rforest} \citep{Schonlau2020}, \texttt{svmachines} \citep{Guenther2018}, and \texttt{parsnip} \citep{Huntington2021}. Indeed, the requirements for compatibility with \stcmd{ddml}\xspace are minimal: Any \texttt{eclass} program with the Stata-typical ``\texttt{reg\,\,y\,\,x}'' syntax, support for \texttt{if} conditions and post-estimation \texttt{predict} is compatible with \stcmd{ddml}\xspace.

\item \stcmd{ddml}\xspace provides flexible multi-line syntax and short one-line syntax. The multi-line syntax offers a wide range of options, guides the user through the DDML algorithm step-by-step, and includes auxiliary programs for storing, loading and displaying additional information. We also provide a complementary one-line version called \stcmd{qddml}\xspace (`quick' \stcmd{ddml}\xspace), which uses a similar syntax as \stcmd{pdslasso}\xspace and \texttt{ivlasso} \citep{Ahrens2018}.
\end{enumerate}

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

\section{DDML with Conditional Unconfoundedness}\label{sec:ddml_confounding}

This section discusses DDML for the Partially Linear Model and the Interactive Model in turn. Both models are special cases of the general causal model \begin{align}\label{eq:all_causes_1}
    Y = f_0(D,\bm{X},U),
\end{align}
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 \eqref{eq:all_causes_1}, $(D, \bm{X}, U)$ jointly determine $Y$, the model is also dubbed the ``all causes model'' (see, e.g., \citealp{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.

\subsection{The Partially Linear Model (\texttt{partial})}\label{sec:plm}
The Partially Linear Model imposes the estimation model
\begin{align}\label{eq:plm}
\begin{aligned}
    Y & = \theta_0 D + g_0(\bm{X}) + U
    \end{aligned}
\end{align}
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 \citet{angrist1999empirical} imply that in the general causal model \eqref{eq:all_causes_1}, $\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{assumption:plm_orthogonal}.\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.}
\begin{assumption}[Conditional Orthogonality]\label{assumption:plm_orthogonal}
$E[Cov(U, D\vert \bm{X})]= 0$ .
\end{assumption}

To show identification of $\theta_0$, consider the score \begin{align}\label{eq:plm_score}
    \psi(\bm{W}; \theta, m,\ell) = \Big(Y - \ell(\bm{X}) - \theta\big(D - m(\bm{X})\big)\Big)\big(D - m(\bm{X})\big),
\end{align}
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
\begin{align*}
    E\left[\psi(\bm{W}; \theta_0, m_0, \ell_0)\right] = 0
\end{align*}
by Assumption \ref{assumption:plm_orthogonal}. When in addition $E[Var(D\vert \bm{X})] \neq 0$, we get
\begin{align}\label{eq:plm_alpha0}
    \theta_0 = \frac{E\left[\big(Y - \ell_0(\bm{X})\big)\big(D - m_0(\bm{X})\big)\right]}{E\left[(D - m_0(\bm{X}))^2\right]}.
\end{align}

Equation \eqref{eq:plm_alpha0} 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 \citet{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 \eqref{eq:plm_score}:
\begin{align*}
    \partial_r \left\{E\left[\psi(W; \theta_0, m_0 + r(m - m_0),\ell_0 + r(\ell - \ell_0))\right]\right\}\vert_{r=0} = 0,
\end{align*}
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 \citet{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 \textit{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$:\begin{align}\label{eq:plm_estimator}
        \hat{\theta}_n = \frac{\frac{1}{n}\sum_{i=1}^n \big(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\big)\big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)}{\frac{1}{n}\sum_{i=i}^n \big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)^2},
\end{align}
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 \citet{Chernozhukov2018}. \citet[Remark~3.1]{Chernozhukov2018} recommend `DML2' over the alternative `DML1' algorithm, which fits the final estimator by fold.}

\noindent\begin{minipage}{\linewidth}
\begin{sttech}[Algorithm 1. DDML for the Partially Linear Model.]
Split the sample $\{(Y_i, D_i, \bm{X}_i)\}_{i=1}^n$ randomly in $K$ folds of approximately equal size. Denote $I_k$ the set of observations included in fold $k$ and $I_k^c$ its complement.
\begin{enumerate}[nosep]
    \item For each $k\in\{1,\ldots, K\}$:
    \begin{enumerate}
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $Y_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{\ell}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $D_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{m}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
    \end{enumerate}
    \item Compute \eqref{eq:plm_estimator}.
\end{enumerate}
\end{sttech}
\end{minipage}

\citet{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)$. \texttt{ddml} computes the DDML estimator for the Partially Linear Model using Stata's \texttt{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 \texttt{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. \citet[][Remark~3.1]{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 \stopt{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 \stopt{fcluster}{varname}.

\subsection{The Interactive Model (\texttt{interactive})}\label{sec:interactive}

The Interactive Model is given by
\begin{align}\label{eq:interactive}
\begin{aligned}
    Y &= g_0(D,\bm{X}) + U
\end{aligned}
\end{align}
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\begin{align*}
    \theta^{\textrm{ATE}}_0 &\equiv E[g_0(1, \bm{X}) - g_0(0,\bm{X})]\\
    \theta^{\textrm{ATET}}_0 & \equiv E[g_0(1, \bm{X}) - g_0(0,\bm{X})\vert D = 1],
\end{align*}
which correspond to the average treatment effect (ATE) and average treatment effect on the treated (ATET), respectively.

Assumptions \ref{assumption:interactive_cmi} and \ref{assumption:interactive_overlap} 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.
\begin{assumption}[Conditional Mean Independence]\label{assumption:interactive_cmi}
$E[U\vert D, \bm{X}]=0$.
\end{assumption}
\begin{assumption}[Overlap]\label{assumption:interactive_overlap}
$\Pr(D=1 \vert \bm{X})\in (0,1)$ with probability 1.
\end{assumption}

Under assumptions \ref{assumption:interactive_cmi} and \ref{assumption:interactive_overlap}, we have\begin{align*}
    E[Y\vert D, \bm{X}] = E[g_0(D,\bm{X})\vert D, \bm{X}] + E[U\vert D, \bm{X}] = g_0(D,\bm{X}),
\end{align*}
so that identification of the ATE and ATET immediately follows from their definition.\footnote{In the defined Interactive Model under Assumption \ref{assumption:interactive_cmi}, 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 \eqref{eq:all_causes_1} that average over both observed and unobserved heterogeneity. See, e.g., \cite{Belloni2017}.}

In contrast to Section \ref{sec:plm}, 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 \citet{Chernozhukov2018} consider is the efficient influence function of \citet{Hahn1998}
\begin{align*}
    \psi^{\textrm{ATE}}(\bm{W};\theta, g, m) =  \frac{D(Y-g(1, \bm{X}))}{m(\bm{X})} {-} \frac{(1-D)(Y - g(0, \bm{X}))}{1 - m(\bm{X})} {+} g(1, \bm{X}) {-} g(0, \bm{X}) {-} \theta,
\end{align*}
where $\bm{W}\equiv (Y, D, \bm{X})$. Similarly for the ATET,
\begin{align*}
    \psi^{\textrm{ATET}}(\bm{W};\theta, g, m, p) =  \frac{D(Y - g(0,\bm{X}))}{p} - \frac{m(\bm{X})(1-D)(Y-g(0,\bm{X}))}{p(1-m(\bm{X}))} - \theta.
\end{align*}
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{assumption:interactive_cmi} and \ref{assumption:interactive_overlap} imply
\begin{align*}E[\psi^{\textrm{ATE}}(\bm{W};\theta_0^{\textrm{ATE}}, g_0, m_0)] &= 0 \\  E[\psi^{\textrm{ATET}}(\bm{W};\theta_0^{\textrm{ATET}}, g_0, m_0, p_0)] &= 0;\end{align*}
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
\begin{align}
\begin{split}\label{eq:ddml_ate}
   \hat{\theta}^{\textrm{ATE}}_n &=\frac{1}{n}\sum_{i=1}^n\Bigg(\frac{D_i(Y_i-\hat{g}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{m}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-D_i)(Y_i-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{m}_{I^c_{k_i}}(\bm{X}_i)} \\ &\hspace{.5\linewidth} +\hat{g}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i)\Bigg),
   \end{split}
   \\
   \begin{split}
      \hat{\theta}_n^{\textrm{ATET}} &=\frac{1}{n}\sum_{i=1}^n\left(\frac{D_i(Y_i-\hat{g}_{I^c_{k_i}}(0,\bm{X}_i))}{\hat{p}} - \frac{\hat{g}_{I^c_{k_i}}(0,\bm{X}_i)(1-D_i)(Y_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i))}{\hat{p}(1-\hat{m}_{I^c_{k_i}}(\bm{X}_i))}\right)  ,\label{eq:ddml_atet}
      \end{split}
\end{align}
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{sec:plm}. 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} \equiv \frac{1}{n}\sum_{i=1}^nD_i$ is the sample share of treated observations.

\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 Appendix \ref{app:algorithms}. Mean and median aggregation over cross-fitting repetitions are implemented as outlined in Remark~2.

\section{DDML with Instrumental Variables}\label{sec:ddml_iv}
This section outlines the Partially Linear IV Model, the Flexible Partially Linear IV Model, and the Interactive IV Model. As in the previous section, each model is a special case of the general causal model \eqref{eq:all_causes_1}. 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 \eqref{eq:plm}, while the Interactive IV Model allows for arbitrary interactions between the treatment $D$ and the controls $\bm{X}$ as in \eqref{eq:interactive}. 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 \citet{Belloni2012} and \citet{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})$.

\subsection{Partially Linear IV Model (\texttt{iv})}\label{sec:pliv}
The Partially Linear IV Model considers the same functional form restriction on the causal model as the Partially Linear Model in Section \ref{sec:plm}. Specifically, the Partially Linear IV Model maintains
\begin{align*}
\begin{aligned}
        Y&= \theta_0 D  + g_0(\bm{X}) + U, \\
\end{aligned}
\end{align*}
where $\theta_0$ is the unknown parameter of interest.\footnote{As in Section \ref{sec:plm}, the interpretation of $\theta_0$ can be generalized under stronger identifying assumptions. See \citet{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{assumption:pliv_orthogonal} and \ref{assumption:pliv_relevance} below are sufficient orthogonality and relevance conditions, respectively, for identification of $\theta_0$.
\begin{assumption}[Conditional IV Orthogonality]\label{assumption:pliv_orthogonal}
$E[Cov(U,Z\vert \bm{X})]= 0$.
\end{assumption}
\begin{assumption}[Conditional Linear IV Relevance]\label{assumption:pliv_relevance}
$E[Cov(D,Z\vert \bm{X})]\neq 0$.
\end{assumption}

To show identification, consider the score function \begin{align*}
    \psi(\bm{W}; \theta, \ell,m,r) = \Big(Y-\ell(\bm{X}) - \theta(D - m(\bm{X}))\Big)\Big(Z - r(\bm{X})\Big),
\end{align*}
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{assumption:pliv_orthogonal} 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 \begin{align}\label{eq:pliv_alpha0}
    \theta_0 = \frac{E\left[\big(Y-\ell_0(\bm{X})\big)\big(Z - r_0(\bm{X})\big)\right]}{E\left[\big(D - m_0(\bm{X})\big)\big(Z - r_0(\bm{X})\big)\right]},
\end{align}
where Assumption \ref{assumption:pliv_relevance} is used to ensure a non-zero denominator.

The DDML estimator based on Equation \eqref{eq:pliv_alpha0} is given by \begin{align}\label{eq:pliv_ddml_estimator}
    \hat{\theta}_{n} = \frac{\frac{1}{n}\sum_{i=1}^n \left(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\right)\left(Z_i-\hat{r}_{I^c_{k_i}}(\bm{X}_i)\right)}{\frac{1}{n}\sum_{i=i}^n \left(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\right)\left(Z_i-\hat{r}_{I^c_{k_i}}(\bm{X}_i)\right)},
\end{align}
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.

\subsection{Flexible Partially Linear IV Model (\texttt{fiv})}\label{sec:hdiv}
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 \citet{Belloni2012} and \citet{Chernozhukov2015a}.
In particular, the Flexible Partially Linear IV Model leverages a conditional mean independence assumption rather than an orthogonality assumption as in Section \ref{sec:pliv}. As in Section \ref{sec:pliv}, we state everything in the case of a scalar $D$.

\begin{assumption}[Conditional IV Mean Independence]\label{assumption:hdiv_iv_meanindependence}
$E[U\vert \bm{Z},\bm{X}] = 0$.
\end{assumption}
Assumption \ref{assumption:hdiv_iv_meanindependence} implies that for any function $\tilde{p}(\bm{Z},\bm{X})$, it holds that  \begin{align}\label{eq:hdiv_moment}
    E\left[\Big(Y - \ell_0(X) - \theta\big(D - m_0(\bm{X})\big)\Big)\Big(\tilde{p}(\bm{Z},\bm{X}) - E\left[\tilde{p}(\bm{Z},\bm{X})\vert \bm{X}\right]\Big)\right] = 0,
\end{align}
where $\ell_0(\bm{X})=E[Y\vert \bm{X}]$ and $m_0(\bm{X})=E[D\vert \bm{X}]$. Identification based on \eqref{eq:hdiv_moment} requires that there exists some function $\tilde{p}$ such that \begin{align}\label{eq:hdiv_ptld_relevance}
    E\left[Cov(D, \tilde{p}(\bm{Z},\bm{X})\vert \bm{X})\right] \neq 0.
\end{align}
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 \eqref{eq:hdiv_ptld_relevance}.\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{assumption:hdiv_relevance} is a consequence of this non-mean independence.
\begin{assumption}[Conditional IV Relevance]\label{assumption:hdiv_relevance}
$E[Var(E[D\vert \bm{Z},\bm{X}]\vert \bm{X})] \neq 0.$
\end{assumption}

Consider now the score function \begin{align*}
    \psi(\bm{W}; \theta, \ell, m, p) = \Big(Y-\ell(\bm{X}) - \theta(D - m(\bm{X}))\Big)\Big(p(\bm{Z}, \bm{X}) - m(\bm{X})\Big),
\end{align*}
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{assumption:hdiv_iv_meanindependence} 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 \begin{align}\label{eq:hdiv_alpha0}
    \theta_0 = \frac{E\left[\big(Y-\ell_0(\bm{X})\big)\big(p_0(\bm{Z}, \bm{X}) - m_0(\bm{X})\big)\right]}{E[\big(D-m_0(\bm{X})\big)\big(p_0(\bm{Z}, \bm{X}) - m_0(\bm{X})\big)]},
\end{align}
where Assumption \ref{assumption:hdiv_relevance} ensures a non-zero denominator.

The DDML estimator based on the moment solution \eqref{eq:hdiv_alpha0} is given by \begin{align}\label{eq:hdiv_ddml_estimator}
        \hat{\theta}_{n} = \frac{\frac{1}{n}\sum_{i=1}^n \big(Y_i-\hat{\ell}_{I^c_{k_i}}(\bm{X}_i)\big)\big(\hat{p}_{I^c_{k_i}}(\bm{Z}_i, \bm{X}_i)-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)}{\frac{1}{n}\sum_{i=i}^n \big(D_i-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)\big(\hat{p}_{I^c_{k_i}}(\bm{Z}_i, \bm{X}_i)-\hat{m}_{I^c_{k_i}}(\bm{X}_i)\big)},
\end{align}
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 \eqref{eq:hdiv_ddml_estimator} 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 \eqref{eq:hdiv_ddml_estimator}.

\noindent\begin{minipage}{\linewidth}
\begin{sttech}[Algorithm 2. LIE-compliant DDML for the Flexible Partially Linear IV Model.]
 Split the sample $\{(Y_i, D_i, \bm{X}_i, \bm{Z}_i)\}_{i=1}^n$ randomly in $K$ folds of approximately equal size. Denote $I_k$ the set of observations included in fold $k$ and $I_k^c$ its complement.
\begin{enumerate}[nosep]
    \item For each $k\in\{1,\ldots, K\}$:
    \begin{enumerate}
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $Y_i$ as the outcome and $\bm{X}_i$ as predictors. Obtain the out-of-sample predicted values $\hat{\ell}_{I^c_k}(\bm{X}_i)$ for $i\in I_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using $D_i$ as the outcome and $(\bm{Z}_i,\bm{X}_i)$ as predictors. Obtain the out-of-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ for $i\in I_k$ and in-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ for $i\in I^c_k$.
            \item Fit a CEF estimator to the sub-sample $I_k^c$ using the in-sample predicted values $\hat{p}_{I^c_k}(\bm{Z}_i,\bm{X}_i)$ as the outcome and $X_i$ as predictors. Obtain the out-of-sample predicted values $\hat{m}_{I^c_k}(\bm{X}_i)$ for $i\in I_k.$
    \end{enumerate}
    \item Compute \eqref{eq:hdiv_ddml_estimator}.
\end{enumerate}
\end{sttech}
\end{minipage}

Standard errors corresponding to $\hat{\theta}_{n}$ in  \eqref{eq:hdiv_ddml_estimator} are the same as in Section \ref{sec:pliv} 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.

\subsection{Interactive IV Model (\texttt{interactiveiv})}\label{sec:interactiveiv}

The Interactive IV Model considers the same causal model as in Section \ref{sec:interactive}; specifically
\begin{align*}
\begin{aligned}
    Y &= g_0(D, \bm{X}) + U \\
\end{aligned}
\end{align*}
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$.

The parameter of interest we target is\begin{align}\label{interactiveiv:fuzzy_late}
\begin{aligned}
   \theta_0 & =  E\left[g_0(1, \bm{X}) - g_0(0, \bm{X})\vert \: p_0(1,\bm{X}) > p_0(0,\bm{X})\right],
\end{aligned}
\end{align}
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 \citet{imbens1994identification}, ``local'' here 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. Under these stronger assumptions, the DDML LATE estimator outlined here targets the conventionally considered LATE parameter.}

Identification again leverages Assumptions \ref{assumption:hdiv_iv_meanindependence} and \ref{assumption:hdiv_relevance} 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.
\begin{assumption}[Monotonicity]\label{assumption:iv_monotonicty}
$p_0(1, \bm{X})\geq p_0(0, \bm{X})$ with probability 1.
\end{assumption}
\begin{assumption}[IV Overlap]\label{assumption:iv_overlap}
$\Pr(Z=1\vert \bm{X})\in(0,1)$ with probability 1.
\end{assumption}

Assumptions \ref{assumption:hdiv_iv_meanindependence}-\ref{assumption:iv_overlap} imply that \begin{align}\label{eq:interactiveiv_flate}
    \theta_0 = \frac{E\left[\ell_0(1, \bm{X})-\ell_0(0, \bm{X})\right]}{E\left[p_0(1, \bm{X})-p_0(0, \bm{X})\right]},
\end{align}
where $\ell_0(Z, \bm{X}) \equiv E[Y\vert Z, \bm{X}]$, verifying identification of the LATE $\theta_0$. Akin to Section \ref{eq:interactive}, however, estimators of $\theta_0$ should not directly be based on Equation \eqref{eq:interactiveiv_flate} because the estimating equations implicit in obtaining \eqref{eq:interactiveiv_flate} do not satisfy Neyman-orthogonality. Hence, a direct estimator of $\theta_0$ obtained by plugging nonparametric estimators in for nuisance functions in  \eqref{eq:interactiveiv_flate} will potentially be highly sensitive to the first step nonparametric estimation error. Rather, we base estimation on the Neyman orthogonal score function \begin{align*}
    \psi(\bm{W};\theta, \ell, p, r) &= \frac{Z(Y - \ell(1, \bm{X}))}{r(\bm{X})} - \frac{(1-Z)(Y-\ell(0, \bm{X}))}{1-r(\bm{X})} + \ell(1, \bm{X})-\ell(0, \bm{X}) \\
    &  + \left[\frac{Z(D - p(1,\bm{X}))}{r(\bm{X})} - \frac{(1-Z)(D-p(0,\bm{X}))}{1-r(\bm{X})} + p(1,\bm{X})-p(0,\bm{X})\right] \times \theta
\end{align*}
where $\bm{W} \equiv (Y, D, \bm{X}, Z)$ . Note that under Assumptions \ref{assumption:hdiv_iv_meanindependence}-\ref{assumption:iv_overlap} 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 \begin{align}
\begin{split}
    \hat{\theta}_n &= \\
    &\frac{\frac{1}{n}\sum_{i}\left(\frac{Z_i(Y_i-\hat{\ell}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{r}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-Z_i)(Y_i-\hat{\ell}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{r}_{I^c_{k_i}}(\bm{X}_i)}+\hat{\ell}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{\ell}_{I^c_{k_i}}(0,\bm{X}_i)\right)}{\frac{1}{n}\sum_{i}\left(\frac{Z_i(D_i-\hat{p}_{I^c_{k_i}}(1,\bm{X}_i))}{\hat{r}_{I^c_{k_i}}(\bm{X}_i)} - \frac{(1-Z_i)(D_i-\hat{p}_{I^c_{k_i}}(0,\bm{X}_i))}{1-\hat{r}_{I^c_{k_i}}(\bm{X}_i)}+\hat{p}_{I^c_{k_i}}(1,\bm{X}_i)-\hat{p}_{I^c_{k_i}}(0,\bm{X}_i)\right)}, \label{eq:ddml_late}
\end{split}
\end{align}
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 only using treated and untreated 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.

\section{The choice of machine learner}\label{sec:stacking}
\citet{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 \citep{bickel2009simultaneous, Belloni2012}, random forests \citep{wager2015adaptive, wager2018, Athey2019a}, neural networks \citep{schmidt2020nonparametric, farrell2021deep}, and boosting \citep{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 \emph{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 \citet{Wolpert1992} and \citet{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.

We introduce stacking for DDML in Section~\ref{eq:stacking}. Section~\ref{eq:stacking_mc} demonstrates the performance of DDML in combination with stacking using a simulation.

\subsection{Stacking}\label{eq:stacking}
Our discussion of stacking in the context of DDML focuses on the Partially Linear Model in \eqref{eq:plm}, 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:
\[  \underset{w_{k,1},\ldots,w_{k,J}}{\min} \sum_{i\in T_k} \left( Y_i - \sum_{j=1}^J w_{k,j} \hat{\ell}^{(j)}_{T_{k,v(i)}^c}(\bm{X}_i)  \right)^2,   \qquad \textrm{s.t.}\ w_{k,j}\geq 0,\ \sum_{j=1}^J |w_{k,j}|=1,\]
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)$ which observation $i$ falls into. 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 \citep{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 \texttt{pystacked} \citep{Ahrens2022}.

\paragraph{Short-stacking.} Stacking relies 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 \emph{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 is the lower computational cost. Short-stacking also produces a single set of weights for the entire sample, rather than a different set of weights in each cross-fit fold and thus facilitates interpretation. 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.}

\subsection{Monte Carlo simulation}\label{eq:stacking_mc}
To illustrate the advantages of DDML with stacking, we generate artificial data based on the Partially Linear Model
\begin{align}
Y_i &= \theta_0D_i + c_Yg(\bm{X}_i) + \sigma_Y(D_i,\bm{X}_i)\varepsilon_i \label{eq:mc1}\\
D_i &= c_Dg(\bm{X}_i) + \sigma_D(\bm{X}_i)u_i \label{eq:mc2}
\end{align}
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 \eqref{eq:mc1} and \eqref{eq:mc2}
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:
\[
\begin{array}{ll}
\textrm{DGP 1:}\quad g(\bm{X}_i) &= \displaystyle\sum\nolimits_j 0.9^jX_{ij} \\
\textrm{DGP 2:}\quad g(\bm{X}_i) &= X_{i1}X_{i2}{+}X_{i3}^2{+}X_{i4}X_{i5}{+}X_{i6}X_{i7}{+}X_{i8}X_{i9}{+}X_{i10}{+} X_{i11}^2{+}X_{i12} X_{i13} \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 3:}\quad g(\bm{X}_i) &= \mathbbm{1}\{X_{i1}>0.3\}\mathbbm{1}\{X_{i2}>0\}\mathbbm{1}\{X_{i3}>-1\}  \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 4:}\quad g(\bm{X}_i) &= X_{i1} + \sqrt{|X_{i2}|} + \textrm{sin}(X_{i3}) + 0.3X_{i4}X_{i5} + X_{i6} + 0.3X_{i7}^2  \phantom{\displaystyle\sum\nolimits_j} \\
\textrm{DGP 5:}\quad g(\bm{X}_i)&=\textnormal{same as DGP 4 with } p=7 \phantom{\displaystyle\sum\nolimits_j} \\
\end{array}\]
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 which 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 \texttt{pystacked}. We use the defaults of \texttt{pystacked} for parameter values and settings not mentioned here.}
\begin{enumerate}[noitemsep,topsep=0pt,itemsep=0pt,parsep=0pt,before=\vspace{0mm},after=\vspace{-0mm}]
    \item[1.-2.] Cross-validated lasso \& ridge with untransformed base controls
    \item[3.-4.] Cross-validated lasso \& ridge with 5th-order polynomials of base controls but no interactions (referred to as `Poly 5')
    \item[5.-6.] Cross-validated lasso \& ridge with second-order polynomials and all first order interaction terms (referred to as `Poly 2 + Inter.')
    \item[7.] Random forests (RF) with low regularization: base controls, maximum tree depth of 10, 500 trees and approximately $\sqrt{p}$ features considered at each split
    \item[8.] RF with medium regularization: same as 7., but with maximum tree depth of 6
    \item[9.] RF with high regularization: same as 7., but with maximum tree depth of 2
    \item[10.] Gradient boosted trees (GB) with low regularization: base controls, 1000 trees and a learning rate of 0.3. We enable early stopping which uses a 20\% validation sample to decide whether to stop the learning algorithm. Learning is terminated after 5 iterations with no meaningful improvement in the mean-squared loss of the validation sample.\footnote{We use a tolerance level of 0.01 to measure improvements.}
    \item[11.] GB with medium regularization: same as 10., but with learning rate of 0.1
    \item[12.] GB with high regularization: same as 10., but with learning rate of 0.01
    \item[13.] Feed-forward neural net with base controls and two layers of size 20
\end{enumerate}

\noindent We use the above set of learners as base learners for DDML with stacking approaches. Specifically, we estimate DDML using stacking regression with CLS, stacking with single-best learner, and short-stacking with CLS. We set the number of folds to $K=20$ if $n=100$, and $K=5$ if $n=1000$. For comparison, we report results for OLS using the base controls, PDS-Lasso with the 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{tab:mc1}. We delegate results for DGPs~4 and~5, including a brief discussion, to  Appendix~\ref{appendix:moremc}. DDML estimators leveraging stacking or short-stacking 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 imprecisely estimated in very 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 four stacking approaches considered. 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~\ref{tab:mc2}). 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~\ref{tab:stacking_weights_appendix} and \ref{tab:mspe_appendix} 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 compare 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.

\begin{sidewaystable}
    \centering\scriptsize\singlespacing
\begin{threeparttable}
    \caption{Bias and Coverage Rates in the Linear and Nonlinear DGPs}
    \label{tab:mc1}
        \begin{tabular}{rlcccccccccccccc}
    \toprule
    \midrule
          &       & \multicolumn{4}{c}{DGP 1} & & \multicolumn{4}{c}{DGP 2} & & \multicolumn{4}{c}{DGP 3} \\
    \cline{3-6} \cline{8-11} \cline{13-16} \\[-6pt]
     & & \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}    &&  \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}     &&  \multicolumn{2}{c}{$n=100$}    &\multicolumn{2}{c}{$n=1000$}       \\
    && MAB & CR.   &  MAB & CR   &&   MAB & CR   &  MAB & CR   &&  MAB & CR   & MAB & CR   \\ \midrule
          \multicolumn{2}{l}{Full sample:}\\
          \partialinput{1}{5}{Simul/sim_SJ/bias_1.tex}
           \multicolumn{2}{l}{DDML methods:}\\
           \multicolumn{2}{l}{\it ~~Base learners}\\
          \partialinput{5}{19}{Simul/sim_SJ/bias_1.tex}
          \multicolumn{2}{l}{\it ~~Meta learners}\\
          \partialinput{20}{23}{Simul/sim_SJ/bias_1.tex}
    \midrule
    \bottomrule
    \end{tabular}%
  \begin{tablenotes}[para,flushleft]
  \footnotesize
  \item \textit{Notes:} The table reports median absolute bias (MAB) and coverage rate of a  95\% confidence interval (CR). We employ standard errors robust to heteroskedasticity.
  For comparison, we report the following full sample estimators: infeasible Oracle, OLS, PDS-Lasso with base and two different expanded sets of covariates.
  DDML estimators use 20 folds for cross-fitting if $n=100$, and 5 folds if $n=1000$. Meta-learning approaches rely on all listed base learners.
  Results are based on 1,000 replications. Results for DGPs 4 and 5 can be found in Table~\ref{tab:mc2} in the Appendix.
  \end{tablenotes}
\end{threeparttable}
\end{sidewaystable}

\section{The program}\label{sec: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{sec:syntax_ddml}. Section~\ref{sec:options} lists the options. Section~\ref{sec:syntax_qddml} covers the simplified one-line program \stcmd{qddml}\xspace. We provide an overview of supported machine learning programs in Section~\ref{sec:supported_programs}. Finally, Section~\ref{sec:replication} adds a note on how to ensure replication with \stcmd{ddml}\xspace.

\subsection{Syntax: \stcmd{ddml}\xspace}\label{sec:syntax_ddml}
The \texttt{ddml} estimation proceeds in four steps.

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

\begin{stsyntax}
ddml init {\it model}
    \optional{,
    mname({\it name})
    vars({\it varlist})
    kfolds({\it integer})
    fcluster({\it varname})
    foldvar({\it varlist})
    reps({\it integer})
    tabfold
    vars({\it varlist})
    }
\end{stsyntax}
\noindent where {\it model} selects between the Partially Linear Model (\texttt{partial}), the Interactive Model (\texttt{interactive}), the Partially Linear IV Model (\texttt{iv}), the Flexible Partially Linear IV Model (\texttt{fiv}), and the Interactive IV Model (\texttt{interactiveiv}). This step creates a persistent Mata object with the name provided by \stopt{mname}{name} 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})}).

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

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

\begin{stsyntax}
ddml {\it cond\_exp}
    \optional{,
    mname({\it name})
    vname({\it varname})
    learner({\it name})
    vtype({\it string})
    predopt({\it string})
    }
    :
    {\it command} {\it depvar} {\it vars}
    \optional{,
    {\it cmdopt}
    }
\end{stsyntax}

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

\begin{table}[H]
\newcommand{\xmark}{}
    \centering
    \begin{tabular}{l|ccccc}\hline\hline
       {\it cond\_exp}      & \texttt{partial} &   \texttt{interactive} & \texttt{iv} & \texttt{fiv} &  \texttt{late} \\ \hline
       \verb+E[Y|X]+ &  \checkmark &  \xmark & \checkmark & \checkmark & \xmark\\
       \verb+E[Y|X,D]+ &\xmark & \checkmark &\xmark & \xmark & \xmark \\
       \verb+E[Y|X,Z]+ & \xmark & \xmark & \xmark & \xmark & \checkmark\\[.1cm]
       \verb+E[D|X]+ & \checkmark & \checkmark & \checkmark & \checkmark & \xmark \\
       \verb+E[D|Z,X]+ & \xmark & \xmark & \xmark & \checkmark & \checkmark  \\[.1cm]
       \verb+E[Z|X]+  & \xmark & \xmark & \checkmark & \xmark & \checkmark \\\hline\hline
    \end{tabular}
    \caption{The table lists the conditional expectations which need to be specified for each model.}
    \label{tab:cond_exp}
\end{table}

\subsubsection*{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.

\begin{stsyntax}
ddml crossfit
    \optional{,
    mname({\it name})
    {\ttfamily shortstack}
    }
\end{stsyntax}

\subsubsection*{Step 4: Estimate causal effects.}

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

\begin{stsyntax}
ddml estimate
    \optional{,
    mname({\it name})
    robust
    cluster({\it varname})
    vce({\it vcetype})
    att
    trim
    spec({\it string})
    rep({\it string})
    }
\end{stsyntax}

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

\begin{stsyntax}
ddml estimate
    \optional{,
    replay
    mname({\it name})
    spec({\it integer or string})
    rep({\it integer or string})
    fulltable
    notable
    allest
    }
\end{stsyntax}

\subsection{Options}\label{sec:options}

\paragraph{Step 1 options: Initialization.}
\begin{description}
    \item[\stopt{mname}{{\it name}}] name of the DDML model. Allows to run multiple DDML models simultaneously. Defaults to {\it m0}.
    \item[\stopt{kfolds}{{\it integer}}] number of cross-fitting folds. The default is 5.
    \item[\stopt{fcluster}{{\it varname}}] cluster identifiers for cluster randomization of folds.
    \item[\stopt{foldvar}{{\it varlist}}] integer variables to specify custom folds (one per cross-fitting repetition).
    \item[\stopt{reps}{{\it integer}}] number of cross-fitting repetitions, i.e., how often the cross-fitting procedure is repeated on randomly generated folds.
    \item[{{\tt tabfold}}]  prints a table with frequency of observations by fold.
\end{description}

\paragraph{Step 2 options: Adding learners.}
\begin{description}
    \item[\stopt{vname}{{\it varname}}] name of the dependent variable in the reduced form estimation.  This is usually inferred from the {\it command} line but is mandatory for the {\tt fiv} model.
    \item[\stopt{learner}{{\it varname}}] optional name of the variable to be created.
    \item[\stopt{vtype}{{\it string}}] optional variable type of the variable to be created. Defaults to {\it double}.  {\it none} can be used to leave the type field blank. (Setting \stopt{vtype}{{\it none}} is required when using \stcmd{ddml}\xspace with {\tt rforest}.)
    \item[\stopt{predopt}{{\it varname}}]  {\tt predict} option to be used to get predicted values.  Typical values could be {\tt xb} or {\tt pr}. Default is blank.
\end{description}

\paragraph{Step 3 options: Cross-fitting.}
\begin{description}
    \item[{\ttfamily shortstack}]  asks for short-stacking to be used.  Short-stacking runs constrained least squares on the cross-fitted predicted values to obtain a weighted average of several base learners.
\end{description}

\paragraph{Step 4 options: Estimation.}
\begin{description}
    \item[\stopt{spec}{{\it string}}] select specification. This can either be the specification number, {\it mse} for minimum-MSE specification (the default) or {\it ss} for short-stacking.
    \item[\stopt{rep}{{\it string}}] select cross-fitting repetitions. This can either be the cross-fit repetition number, {\it mn} for mean aggregation or {\it md} for median aggregation (the default). See Remark~2 for more information.
    \item[{\tt robust}] report SEs that are robust to the presence of arbitrary heteroskedasticity.
    \item[\stopt{cluster}{{\it varname}}] select cluster-robust variance-covariance estimator.
    \item[\stopt{vce}{{\it type}}] select variance-covariance estimator, e.g.\ \texttt{vce(hc3)} or \texttt{vce(cluster id)}. See \verb+help regress##vcetype+ for available options.
    \item[\stopt{trim}{{\it real}}] trimming of propensity scores for the Interactive and Interactive IV models. The default is 0.01 (that is, values below 0.01 and above 0.99 are set to 0.01 and 0.99, respectively).
    \item[{\tt atet}] report average treatment effect of the treated (default is ATE).
    \item[{\tt noconstant}] suppress constant term in the estimation stage (only relevant for partially linear models).
\end{description}

\subsection{Short syntax: \stcmd{qddml}\xspace}\label{sec:syntax_qddml}
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 \texttt{pdslasso} and \texttt{ivlasso} \citep{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 \texttt{regress}, which is added by default).

\subsubsection*{Syntax for Partially Linear and Interactive Model}
\begin{stsyntax}
qddml
    \depvar\
    {\it treatment\_vars}
    \textnormal{(}{\it controls}\textnormal{)}\verb+,+
    \verb+model(partial|interactive)+
    \optional{
    {\it options}
    }
\end{stsyntax}

\subsubsection*{Syntax for IV models}
\begin{stsyntax}
qddml
    \depvar\
    \textnormal{(}{\it controls}\textnormal{)}
    \textnormal{(}{\it treatment\_vars}=excluded\_instruments\textnormal{)}
    \verb+,+
    \verb+model(iv|late|fiv)+
    \optional{
    cmd({\it string})
    cmdopt({\it string})
    {\it ddml\_options}
    {\tt noreg}
    }
\end{stsyntax}

\noindent where {\it ddml\_options} options are listed in Section~\ref{sec:options} and are internally passed on to the \stcmd{ddml}\xspace sub-routines. \stopt{cmd}{string} selects the machine learning command to be used and \stopt{cmdopt}{string} allows to pass options to the estimation command. \stcmd{qddml}\xspace adds by default \stcmd{regress} to the user-specified learner; this behavior can be disabled with \stcmd{noreg}. See the \stcmd{qddml}\xspace help file for a full list of options.

\subsection{Supported machine learning programs}\label{sec:supported_programs}
\stcmd{ddml}\xspace is compatible with any supervised ML program in Stata that supports the typical ``\texttt{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:
\begin{itemize}[nosep]
    \item \texttt{lassopack} implements regularized regression, e.g.\ lasso, ridge, elastic net \citep{Ahrens2019}.
    \item \stcmd{pystacked}\xspace facilitates the stacking of a wide range of machine learners including regularized regression, random forests, support vector machines, gradient boosted trees and feed-forward neural nets using Python's \emph{scikit-learn} \citep{Ahrens2022,scikit-learn,sklearn_api}. In addition, \stcmd{pystacked}\xspace can also be used as a front-end to fit individual machine learners.

    \item \texttt{rforest} is a random forest wrapper for WEKA \citep{Schonlau2020,frank2009weka}.
    \item \texttt{svmachines} allows for the estimation of support vector machines using \texttt{libsvm} \citep{chang2011libsvm,Guenther2018}.
    \item The program \texttt{parsnip} of the package \texttt{mlrtime} provides access to R's \emph{parsnip} machine learning library through \texttt{rcall} \citep{Huntington2021,haghish2019seamless}. Using \texttt{parsnip} requires the installation of the supplementary wrapper program \texttt{parsnip2}.\footnote{Available from {\tt\url{https://github.com/aahrens1/parsnip2}.}}
\end{itemize}

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

\subsection{Inspecting results and replication}\label{sec: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_+\textit{r}, where \verb+m0+ is the default model name and \textit{r} is the cross-fitting repetition). \stcmd{ddml}\xspace models are stored on Mata \texttt{struct}s and using Mata's associative arrays. Specifically, the \stcmd{ddml}\xspace model created by \stcmd{ddml}\xspace \texttt{init} is an \texttt{mStruct}, and information relating to the estimation of conditional expectations are stored in \texttt{eStruct}s. 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:
\begin{itemize}
    \item Set the Stata seed before {\tt ddml init}. This ensure that the same random fold variable is used for a given data set.
    \item Using the same fold variable alone is usually not sufficient to ensure replication, since many machine learning algorithms involve randomization. That said, note that the fold variable is stored in memory and can be reused for subsequent estimations via the \stopt{foldvar}{varlist} option.
    \item Replication of \stcmd{ddml}\xspace results may require additional steps with some programs that rely on randomization in other software environments, e.g., R or Python. \stcmd{pystacked}\xspace uses a Python seed generated in Stata. Thus, when \stcmd{ddml}\xspace is used with \stcmd{pystacked}\xspace, setting the seed before {\tt ddml init} also guarantees that the same Python seed underlies the stacking estimation. Other programs relying on randomization outside of Stata might not behave in the same way. Thus, when using other programs, check the help files for options to set external random seeds. Try estimating each of the individual learners on the entire sample to see what settings need to be passed to them for their results to replicate.

    \item The \stcmd{ddml}\xspace \texttt{extract} utility can be used to retrieve and inspect a wide range of intermediate results and statistics from Mata structures. This can be useful for trying to track down estimation errors and replication anomalies.
    \item Beware of changing samples. Fold splits or learner idiosyncracies may mean that sample sizes vary slightly across learners, estimation samples and/or cross-fitting repetitions. Note that \stcmd{ddml}\xspace stores sample indicators for cross-fitting repetitions as Stata variables; \stcmd{ddml}\xspace \texttt{extract} with the \texttt{show(n)} option will report sample sizes by learner and fold. See the {\tt ddml extract} help file for more information.
    \item The \stcmd{ddml}\xspace \texttt{export} utility can be used to export the estimated conditional expectations, fold variables and sample indicators to a CSV format file for examination and comparison in other software environments.
\end{itemize}

\section{Applications}\label{sec:application}
We demonstrate the \texttt{ddml} workflow using two applications. In Section~\ref{sec:401k}, we apply the DDML estimator to estimate the effect of 401(k) eligibility on financial wealth following \citet{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~\ref{appendix:moreapplications}. Additional examples can also be found in the help file.
Based on \citet{berry1995automobile}, we show in Section~\ref{sec:ivhd_application} 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.

\subsection{401(k) and financial wealth}\label{sec:401k}
The data consists of $n=9915$ households from the 1991 SIPP. The application is originally due to \citet{poterba1995}, but has been revisited by \citet{Belloni2017}, \citet{Chernozhukov2018}, and \citet{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{sec:syntax_ddml}.

\begin{stlog}
\input{Applications/logs/partial_start.log.tex}
\end{stlog}

\subsubsection*{Step 1: Initialize \texttt{ddml} model.} We initialize the \stcmd{ddml}\xspace model and select the Partially Linear Model in \eqref{eq:plm}. Before initialization, we set the seed to ensure replication. This should always be done before \texttt{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.

\begin{stlog}
\input{Applications/logs/partial_init.log.tex}
\end{stlog}

\subsubsection*{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 \texttt{regress} for linear regression, \texttt{cvlasso} for cross-validated lasso and \texttt{rforest} for random forests. When using \texttt{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.

\begin{stlog}
\input{Applications/logs/partial_add.log.tex}
\end{stlog}

The flexible \stcmd{ddml}\xspace syntax allows specifying 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, for example, 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 \texttt{cvlasso}, 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 \texttt{ddml E[D|X]: reg D2 \$X} where \texttt{D2} would be the name of the second treatment variable. An example with two treatments is provided in the help file.

The auxiliary \texttt{ddml} sub-command \texttt{describe} allows to verify that the learners were correctly registered:

\begin{stlog}
\input{Applications/logs/partial_describe.log.tex}
\end{stlog}

\subsubsection*{Step 3: Perform cross-fitting.} The third step performs cross-fitting, which is the most time-intensive process. The \texttt{shortstack} option enables the short-stacking algorithm of Section~\ref{eq:stacking}.

\begin{stlog}
\input{Applications/logs/partial_crossfit.log.tex}
\end{stlog}

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_cvlasso_1+, \verb+Y3_rforest_1+, \verb+D1_reg_1+, \verb+D2_cvlasso_1+ and \verb+D3_rforest_1+. \texttt{Y} and \texttt{D} indicate outcome and treatment variable. The index \texttt{1} to \texttt{3} is a learner counter. \texttt{reg}, \texttt{cvlasso} and \texttt{rforest} correspond to the name of the commands used. The \verb+_1+ suffix  indicates the cross-fitting repetition.

After cross-fitting, we can inspect the mean-squared prediction errors by fold and learner:
\begin{stlog}
\input{Applications/logs/partial_extract.log.tex}
\end{stlog}

\subsubsection*{Step 4: Estimate causal effects.}

In this final step, we obtain the causal effect estimates. Since we requested short-stacking in Step~3, \texttt{ddml} shows 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 (denoted as \texttt{opt}).

\begin{stlog}
\input{Applications/logs/partial_estimate.log.tex}
\end{stlog}

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_cvlasso_1+ and \verb+Y3_rforest_1+ with \verb+D1_reg_1+, \verb+D2_cvlasso_1+ and \verb+D3_rforest_1+). To get all results, we add the \texttt{allcombos} option:

\begin{stlog}
\input{Applications/logs/partial_estimate_all.log.tex}
\end{stlog}

We can use the \stopt{spec}{string} option to select among the listed specifications. {\it string} is either the specification number, \texttt{ss} to get the short-stacking specification 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 \stopt{spec}{string} option can be provided either in combinations with \texttt{allcombos}, or after estimation in combination with the {\tt replay} option, for example:

\begin{stlog}
\input{Applications/logs/partial_spec.log.tex}
\end{stlog}

\subsubsection*{Manual final estimation.}
In the background, \texttt{ddml estimate} regresses \verb+Y1_reg_1+ against \verb+D1_reg_1+ with a constant. We can verify this manually:

\begin{stlog}
\input{Applications/logs/partial_reg.log.tex}
\end{stlog}

Manual estimation using {\ttfamily regress} allows the use of postestimation regression tools such as {\ttfamily avplot}.

\subsubsection*{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. For demonstration, we use \texttt{cvlasso}:

\begin{stlog}
\input{Applications/logs/partial_qddml.log.tex}
\end{stlog}

\noindent \stcmd{qddml}\xspace always includes linear regression (\stcmd{regress}) as an additional learner. The option \stcmd{noreg} disables this default behavior.

\subsubsection*{Using \stcmd{pystacked}\xspace with single learners}

The above example relying on \texttt{cvlasso} and \texttt{rforest} runs relatively slowly at 380 seconds. We can substantially improve the speed by using the lasso and random forests implementation offered by \stcmd{pystacked}\xspace, which calls the Python library \emph{scikit-learn}. This reduces the total computational time to only 29 seconds.

\begin{stlog}
\input{Applications/logs/partial_2_add.log.tex}
\end{stlog}

\noindent Note that, despite its name, \stcmd{pystacked}\xspace does not execute stacking regression in this example, since only one learner is specified in each call to \stcmd{pystacked}\xspace.

\subsubsection*{Stacking}

We next demonstrate DDML with stacking. To this end, we exploit the stacking regressor implemented in \stcmd{pystacked}\xspace. \stcmd{pystacked}\xspace allows to combine multiple base learners with learner-specific settings and covariates into a final meta learner. The learners are separated by \verb+||+. \stopt{method}{name} selects the learner, \stopt{xvars}{varlist} specifies learner-specific covariates (overwritting the default covariates \verb+$X+) and \stopt{opt}{string} 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 \citet{Ahrens2022}.

\begin{stlog}
\input{Applications/logs/partial_stack_add.log.tex}
\end{stlog}

\noindent After cross-fitting, we can obtain the MSE and stacking weights averaged over folds:
\begin{stlog}
\input{Applications/logs/partial_stack_crossfit.log.tex}
\end{stlog}

By adding the \texttt{detail} option, the user can display the stacking weights for each fold.
 Finally, we retrieve the results of DDML with stacking:
\begin{stlog}
\input{Applications/logs/partial_stack_estimate.log.tex}
\end{stlog}

\noindent The run time of this example is 83 seconds, which despite its computational complexity is faster than using \stcmd{ddml}\xspace with \texttt{cvlasso} and \texttt{rforest}.

\subsection{The market for automobiles}\label{sec:ivhd_application}
For this demonstration, we follow \citet{chernozhukov2015} who estimate a stylized demand model using instrumental variables based on the data from \citet{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, \citet{berry1995automobile} suggest exploiting other products' characteristics as instruments. Following \citet{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. \citet{chernozhukov2015} also construct an augmented set of instruments, including first-order interactions, squared and cubic terms. In the analysis below, we extend \citet{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:
\begin{stlog}
\input{Applications/logs/ivhd_start.log.tex}
\end{stlog}

\paragraph{Step 1: Initialize \texttt{ddml} model.}

\begin{stlog}
\input{Applications/logs/ivhd_init.log.tex}
\end{stlog}

Note that in the \texttt{ddml init} step, we include the option \texttt{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 \texttt{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 \stopt{learner}{string} option to match one learner for $E[D|\bm{X}]$ with a learner for $E[D|\bm{X},\bm{Z}]$, and \stopt{vname}{varname} to indicate the name of the treatment variable.

\begin{stlog}
\input{Applications/logs/ivhd_add.log.tex}
\end{stlog}

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

\begin{stlog}
\input{Applications/logs/ivhd_estimate.log.tex}
\end{stlog}

\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 ``\texttt{\_1}'' to ``\texttt{\_}$r$''.

\begin{stlog}
\input{Applications/logs/ivhd_ivreg.log.tex}
\end{stlog}

\paragraph{One-line syntax (output omitted).}

\begin{stlog}
\input{Applications/logs/ivhd_qddml.log.tex}
\end{stlog}

\section{Acknowledgments}
We thank users who tested earlier versions of the program. We also thank Jan Ditzen, Ben Jann, Eroll Kuhn and Di Liu for helpful comments, as well as participants at the German Stata Conference 2021, Italian Stata Conference 2022 and Swiss Stata Conference 2022. All remaining errors are our own.

\bibliographystyle{sj}