Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
104,366 characters · 14 sections · 37 citation commands
Double/Debiased Machine Learning for Dynamic Treatment Effects via $g$-Estimation
\etocdepthtag.toc{mtchapter} \etocsettagdepth{mtchapter}{subsection} \etocsettagdepth{mtappendix}{none}
Improving outcomes often requires multiple treatments: patients may need a course of drugs to manage or cure a disease, soil may need multiple additives to improve fertility, companies may need multiple marketing efforts to close the sale. To make data-driven decisions, policy-makers need precise estimates of what will happen when a new policy is pursued. Because of its importance, this topic has been studied by many communities and under multiple regimes and formulations; examples include the field of reinforcement learning Sutton1998, longitudinal analysis in biostatistics hernan2010causal, the dynamic treatment regime in causal inference and adaptive clinical trials lei2012smart.
This paper offers a new method for estimating and making inferences about the counterfactual effect of a new treatment policy. Our method is designed to work with observational data, in an environment with multiple treatments (either discrete or continuous), and a high-dimensional state. Valid causal inference is necessary to correctly attribute changes in outcomes to the different treatments applied. But it is more challenging than in a static context, since there are multiple causal pathways from treatments to subsequent outcomes (e.g. directly, or by changing future states, or by affecting intermediate outcomes, or by influencing future treatments).
Our work bridges many distinct literatures. The first is the econometrics literature on semi-parametric inference Neyman:1979,robinson:88,newey1990semiparametric,Ai2003,Chernozhukov2016locally,chernozhukov2018double. We extend this literature, which has typically focused on static treatment regimes, to dynamic treatment regimes and to estimation of counterfactual dynamic policies. We propose an estimation algorithm that estimates the dynamic effects of interest, from observational data, at parametric root-$n$ rates. We prove asymptotic normality of the estimates, even when the state is high-dimensional and machine learning algorithms are used to control for indirect effects of the state, as opposed to effects stemming from the treatments. Our formulation can be viewed as a dynamic version of Robinson's classic partial linear model in semi-parametric inference robinson:88, where the controls evolve over time and are affected by prior treatments. Our estimation algorithm builds upon and generalizes the recent work of chernozhukov2018double,semenova:17 to estimate not only contemporaneous effects, but also dynamic effects over time. In particular, we propose a sequential residualization approach, where the effects at every period are estimated in a Neyman orthogonal manner and then peeled-off from the outcome, so as to define a new “calibrated” outcome, which will be used to estimate the effect of the treatment in the previous period.
In doing this, we build on results in the semi-parametric inference literature in bio-statistics bickel1993efficient,van2003unified,tsiatis2007semiparametric on estimating causal effects in structural nested models (see Vansteelandt2016,Robins2004,hernan2010causal,Chakraborty2013 for recent overviews). In particular, our identification and estimation strategy for dynamic treatment effects is a variant of the well-studied $g$-estimation framework for structural nested mean models (SNMMs) robins1986new,robins1992g,robins1994correcting,robins1997toward,robins2000marginal,Robins2004,lok2012impact,vansteelandt2014structural. The works cited above have provided doubly robust estimators for this setting. One practical challenge faced by these estimation approaches is that the nuisance functions required are hard to estimate when there are many treatments or treatments are continuous and they do not yield and obvious way to perform cross-fitting, which is crucial for enabling the use of machine learning approaches for nuisance estimation. Moreover, typically $g$-estimation has been used in practice with either logistic regresion or ordinary least squares as nuisance function estimators.
Our approach addresses these challenges by providing a Neyman orthogonal (aka locally robust) $g$-estimation algorithm for linear structural nested mean models robins1994correcting,Robins2004, that allows for both continuous and discrete treatments in each time period. We propose a two stage algorithm, where in the first stage a sequence of regression and classification models are fitted and evaluated in a cross-fitting manner and in the second stage a simple linear system of equations or a simple square loss minimization problem is solved. This approach also allows for an easy sample-splitting/cross-fitting approach, which allows the use of arbitrary machine learning approaches in the first stage.
While Neyman orthogonality is a weaker condition than double robustness, it is sufficient to achieve robustness to bias introduced by using machine learning to estimate the nuisance parameters. This approach thus permits machine learning for the nuisance model estimation, which is practically important in high-dimensional state space/control variable settings. Our work extends and enriches the $g$-estimation framework for dynamic treatment regimes with arbitrary target dynamic policies in multiple ways:
Many of our aforementioned advancements stem from the fact that the Neyman orthogonal moment that we propose at each stage of $g$-estimation is the gradient of a strongly convex square loss. This allows for stronger finite sample learning rates, as well as the use of techniques from the recent framework of orthogonal statistical learning foster2019orthogonal,chernozhukov2018plugin. We extend these techniques to also account for the recursive nature of the final stage estimation and the fact that each of these steps in the final stage cannot be viewed in isolation, since Neyman orthgoonality does not hold within structural parameters fitted at different steps of the final stage and we also have sample re-use.
Our work is also closely related to the work on doubly robust estimation in offline policy evaluation in reinforcement learning nie2019learning,Thomas2016,petersen2014targeted,kallus2019efficiently,kallus2019double from the causal machine learning community. However, one crucial point of departure from all of these works is that we formulate minimal parametric assumptions that allow us to avoid non-parametric rates or high-variance estimation. Typical fully non-parametric approaches in off-policy evaluation in dynamic settings, requires the repeated estimation of inverse propensity weights at each step, which leads to a dependence on quantities that can be very ill-posed in practice, such as the product of the ratios of states under the observational and the target policy. Moreover, these approaches are typically restricted to discrete states and actions. Our goal is to capture settings where the treatment can also be continuous (investment level, price discount level, drug dosage), and the state contains many continuous variables and is potentially high dimensional. Inverse propensity weighting approaches, though more robust in terms of the assumptions that they make on the model, can be quite prohibitive in these settings even in moderately large data sets.
Our work is also somewhat related to the online debiasing literature deshpande2018accurate,deshpande2019online,Zhang2020,zhan2021policy,Hadad2021,bibaut2021post,zhan2021off, since we perform inference from adaptively collected data. However, our inferential target is different (dynamic effects vs same-period effects). Moreover, in our main setting, we assume multiple independent small chains of samples (each from a separate treated unit), rather than one long time series, while in our single time series setting, we make stationarity assumptions that avoid the main problems from adaptive data collection that are being handled by that literature. It is an interesting avenue for future work to combine techniques from this literature and lift some of our stationarity conditions in the estimation of dynamic treatment effects from a single time series.
We begin by presenting the main algorithm of this work in the context of a high-dimensional partially linear Markovian data generating process. In Section (ref) we show that our results generalize to more complex dynamic treatment models, known in the bio-statistics and causal inference literature as Structural Nested Mean Models (SNMMs), but for simplicity of exposition we focus on this simplified version, as it captures all the complexities needed to highlight our main contributions.
Consider a partially linear state space Markov decision process $\{X_t, T_t, Y_t\}_{t=1}^{m}$, where $X_t\in \mathbb{R}^p$ is the state at time $t$, $T_t\in \mathbb{R}^d$ is the action or treatment at time $t$ and $Y_t \in \mathbb{R}$ is an observed outcome of interest at time $t$. We assume that these variables are related via a linear Markovian process:
where $\eta_t, \zeta_t$ and $\epsilon_t$ are exogenous mean-zero random shocks, independent of all contemporaneous and lagged treatments and states, that for simplicity we assume are each drawn i.i.d. across time. Moreover, for simplicity we assume $T_0=X_0=0$. In Section (ref), we will substantially drop the heavy assumptions on the exogenous shocks and merely assume that the observational process satisfies the quite permissive notion of conditional sequential exogeneity, which essentially only requires that conditional on the past history, the treatment is randomized in an exogenous manner and there is no unobserved confounder at each stage of the observational decision process.
The structural parameter $A$ is a $p\times d$ matrix that governs how past treatments affect next period's states. The structural parameter $B$ is a $p\times p$ matrix that governs how past states affect next period's states, i.e. how the system evolves in the absence of treatments. The function $p(T_{t-1}, X_t, \zeta_t)$ is the observational dynamic policy that determines the distribution of next period's treatments as a function of past periods treatments and states. Finally, the parameter vector $\theta_0\in \mathbb{R}^d$ is the contemporaneous effect of the treatments and $\mu\in R^p$ is the contemporaneous effect of the states on the outcome.
Our goal is to estimate the effect of a change in the treatment policy on the final outcome $Y_m$. Since throughout the analysis we will primarily care about the final outcome, we denote it for simplicity as:
A similar analysis can be derived for any other linear combination of the outcomes at all periods. More concretely, suppose that were to make an intervention and set each of the treatments to some sequence of values: $\{\tau_1, \ldots, \tau_m\}$, then what would be the expected difference in the final outcome $Y_m$ as compared to some baseline policy? For simplicity and without loss of generality, we will consider the baseline policy to be setting all treatments to zero. We will denote this expected difference as: $V(\tau_1, \ldots, \tau_m)$. Equivalently, we can express the quantity we are interested in do-calculus: if we denote with:
In Section (ref), we will also analyze the estimation of the effect of adaptive counterfactual treatment policies, where the treatment at each step can be a function of the state.
Our first observation is that we can decompose the quantity $V(\tau)$ into the estimation of the dynamic treatment effects: if we were to make an intervention and increase the treatment at period $t$ by $1$ unit, then what is the change $\psi_{t}$ in the outcome $Y_{m}$, for $t \in \{1,\ldots, m\}$; assuming that we set all subsequent treatments to zero (or equivalently to some constant value). This quantity is the effect in the final outcome, that does not go through the changes in the subsequent treatments, due to our observational Markovian treatment policy, but only the part of the effect that goes through changes in the state space $X_t$, that is not part of our decision process. This effect can also be expressed in terms of the constants in our Markov process as:
Thus to estimate the function $V$, it suffices to estimate the dynamic treatment effects: $\psi_1, \ldots, \psi_{m}$. We first start by showing that the parameters $\psi_1, \ldots, \psi_m$ are identifiable from the observational data. Identification is not immediately obvious. If we write the final outcome as a linear function of all the treatments and the initial state by repeatedly expanding the structural equations, we will arrive at an equation of the form:
However, the shocks $\{\eta_j\}_{j=t+1}^m$ are heavily correlated with the treatments $T_2,\ldots, T_m$, since the shocks at period $t$ affect the state at period $t+1$, which in turn affects the observed treatment at period $t+1$. As a result the natural moment condition that the sum of shocks is conditionally mean zero $\mathbb{E}[Y - \sum_{t=0}^{m} \psi_t' T_{t} - \mu'B^{m} X_1 \mid X_1, T_1, \ldots, T_m]=0$ is not valid. In other words, if we view the problem as a simultaneous treatment problem, where the final outcome is the outcome, then we essentially have a problem of unmeasured confounding (implicitly because we ignored the confounding through intermediate states, sometimes referred as the “treatment-confounder feedback”).
However, note that if we apply this recursive expansion process up until any period $t=\{1,\ldots, m\}$, then we can write:
Since the random shocks $\{\eta_j\}_{j=t+1}^m$ and $\epsilon_m$ are independent of $T_t, X_t$ and mean zero, we thus have that the following conditional moment restriction is satisfied:
This leads to an identification strategy of the dynamic effects via a recursive peeling process, which as we show in the next section, leads to an estimation strategy that achieves parametric rates.
We now address the estimation problem. We assume that we are given access to $n$ i.i.d. samples from the Markovian process, i.e. we are given $n$ independent time-series, and we denote sample $i$, with $\{X_{t}^{i}, T_t^i, Y_t^i\}$. Our goal is to develop an estimator of the function $V$ or equivalently of the parameter vector $\psi=(\psi_{1},\ldots, \psi_m)$. We will consider the case of a high-dimensional state space, i.e. $p\gg n$, but low dimensional treatment space and a low dimensional number of periods $m$, i.e. $d, m\ll n$ is a constant independent of $n$. We want to estimate the parameters $\psi$ at $\sqrt{n}$-rates and in a way that our estimator is asymptotically normal, so that we can construct asymptotically valid confidence intervals around our dynamic treatment effects and our estimate of the function $V$. The latter is a non-trivial task due to the high-dimensionality of the state space. For instance, the latter would be statistically impossible if we were to take the direct route of estimating the whole Markov process (i.e. the high-dimensional quantities $A, B, \mu$): if these quantities have a number of non-zero coefficients that grows with $n$ at any polynomial rate, then known results on sparse linear regression, preclude their estimation at root-n rates (see e.g. Wainwright15). However, we are not really interested in these low-level parameters of the dynamic process, but solely on the low dimensional parameter vector $\theta$. We will treat this problem as a semi-parametric inference problem and develop a Neyman orthogonal estimator for the parameter vector Neyman:1979,robinson:88,Ai2003,Chernozhukov2016locally,chernozhukov2018double.
In particular, we consider a sequential version of the double machine learning algorithm proposed in chernozhukov2018double. In the case of a single time-period, i.e. $m=0$, then chernozhukov2018double, recommends the following estimator for $\psi_m := \theta_0$: using half of your data, fit a model $\hat{q}_0(X_0)$ of $\mathbb{E}[Y_0 \mid X_0]$, i.e. that predicts the outcome $Y_0$ from the controls $X_0$ and a model $\hat{p}_0(X_0)$ for $\mathbb{E}[T_0 \mid X_0]$. Then estimate $\theta_0$ on the other half of the data, based on the estimating equation:
where $\tilde{Y}_0=Y_0-\hat{q}_0(X_0)$ and $\tilde{T}_0 = T_0 - \hat{p}_0(X_0)$ are the residual outcome and treatment.
We propose a sequential version of the double machine learning process that we call Dynamic DML for dynamic double/debiased machine learning. Intuitively our algorithm proceeds as follows:
We provide a formal statement of the Dynamic DML process in Algorithm (ref), which also describes more formally the sample splitting and cross-fitting approach that we follow in order to estimate the nuisance models $p$ and $q$ required for calculating the estimated residuals.
We show that subject to the first stage models of the conditional expectations achieving a small (but relatively slow) estimation error, then the recovered parameters are root-n-consistent and asymptotically normal.
In the static case, the estimate that is based on the estimating Equation (ref) is a special case of a broader class of moment based $Z$-estimators, where the true parameter $\theta$ is known to satisfy a vector of moment restrictions: $\mathbb{E}[m(W;\theta, \nu)]=0$, where $W$ is the vector of all random variables and $\nu\in {\mathcal V}$ is an unknown (potentially infinite dimensional) nuisance parameter, which we do not care about, but on which the moment conditions depends.A plug-in $Z$-estimate $\hat{\theta}$ is a solution to an empirical analogue of a vector of moment equations $\mathbb{E}_n[m(W;\theta, \hat{\nu})] = \frac{1}{n}\sum_{i=1}^n m(W^i;\theta, \hat{\nu})$, where $\hat{\nu}$ is some estimate of the nuisance parameter obtained via some separate statistical learning process and potentially on a separate sample.
A vector of moments satisfies Neyman orthogonality (aka local robustness) if:
Neyman orthogonality implies that small perturbations to the nuisance functions around their true values only has a second order effect on the moment function and hence cannot impact a lot the target parameter estimate. Neyman orthogonality (accompanied with sample splitting techniques) allows one to estimate the target parameter at $\sqrt{n}$-rates and with an asymptotic normal distribution, subject only to much slower mean squared error rates for the nuisance components. Thus allowing machine learning approaches to be used for nuisance estimation.
Our asymptotic normality proof relies on showing that one can re-interpret our Dynamic DML algorithm as a $Z$-estimator based on a set of moments that satisfy the property of Neyman orthogonality. Our finite sample $\ell_2$-error result uses the fact that each step of the recursive process in the final stage is a minimization of a strongly convex loss and crucially invokes Neyman orthogonality in arguing that the true parameter is an approximate first order optimal of the population analogue of the strongly convex loss. The further difficulty in the finite sample result is that the loss at each stage is biased due to the errors propagating from previous stage estimates, with respect to which Neyman orthogonality is not satisfied and which are not constructed in a cross-fitting manner, which introduces sample re-use considerations. We provide a recursive upper bound on the $\ell_2$-errors and complete the theorem by induction.
To present the theorems, we introduce some notation. Let ${h}=\{{p}_{j,t}, {q}_t\}_{1\leq t\leq j\leq m}$ denote the vector of all nuisance functions and $h^*=\{{p}_{j,t}^*, {q}_t^*\}_{1\leq t\leq j\leq m}$ their corresponding true values. Moreover, let $\hat{\psi}=(\hat{\psi}_1,\ldots,\hat{\psi}_m)$ denote the vector of dynamic effect parameter estimates, and $\psi^*=(\psi_1^*, \ldots, \psi_m^*)$ their corresponding true values. We provide both finite sample $\ell_2$-error rates and asymptotic normality of our estimates (proofs in Appendix (ref) and (ref)), subject to mean squared error guarantees for the nuisance functions. At the end of the section we discuss how the required guarantees for the nuisance functions can be easily satisfied using estimation algorithms such as the Lasso and under sparsity conditions.
\paragraph{Concrete Rates for Lasso Nuisance Estimates.} Suppose that the observational policy $p$ is also linear, i.e.
for some $d\times p$ matrix $\Gamma$. Then all the models $q_{t}$ and $p_{j, t}$ are high-dimensional linear functions of their input arguments, i.e. $q_{t}(x)= \phi_t'x$ and $p_{j, t}(x) = \Pi_{j, t} x$. If these linear functions satisfy a sparsity constraint then under standard regularity assumptions we can guarantee if we use the Lasso regression to estimate each of these functions that w.p. $1-\delta$, the estimation error of all nuisance models is $O\left(s\sqrt{\frac{\log(p/\delta)}{n}}\right)$, where $s$ is an upper bound on the number of non-zero coefficients. One sufficient regularity condition is that the expected co-variance matrix of every period's state has full rank, i.e. $\mathbb{E}[X_{t} X_{t}']\succeq \lambda I$ (we note that for an MSE rate we do not require the minimum eigenvalue condition, albeit then a computationally inefficient, support enumeration based estimation algorithm needs to be used and the Lasso results would not apply). Thus the requirements of the main theorems of this section would be satisfied as long as the sparsity grows as $s=o(n^{1/4})$, so that the error from the nuisance estimates is of second order importance. These sparsity conditions are for instance satisfied if only $s$ coordinates of the high-dimensional state, which coevolve separately from the remainder states, have any effect on the final outcome (i.e. are outcome-relevant), and similarly if only $s$ coordinates of the high-dimensional state, which coevolve separately from the remainder states, enter the observational policy.
\paragraph{The Constants $\lambda, \kappa$.} Given the potentially cryptic nature of some of the constants in our main theorems, we connect them here to some more low level quantities in the data generating process for the case of linear policies, i.e. under Equation (ref). There are two main constants $\kappa$ and $\lambda$ that govern our finite sample error rates. Especially parameter $\kappa$ greatly impacts the estimation rate, as a function of the number of iterations $m$, since for $\kappa < 1$ the estimation error grows polynomially with the number of rounds, while for $\kappa > 1$ it grows exponentially. Thus understanding which regime occurs in a setting of interest is of great importance in understanding whether long time sequences can be tolerated. For simplicity of the calculations and exposition, we will further assume a scalar treatment (or binary treatment), i.e. $d=1$, in which case we will write: $p(x, \zeta)=\gamma'x + \zeta$. Moreover, we note that in this case the matrix $A$ in Equation (ref) is a column vector and we will denote it with $\alpha$.
First we note that the parameter $\lambda$ can be easily characterized under the linear policy as:
Thus $\lambda=\mathbb{E}[\zeta^2]$ is the variance of the exploration/randomization of the observational policy deployed at each round $t$. Now let us analyze the constant $\kappa$. First we need to understand the conditional covariances $\ensuremath{\mathtt{Cov}}(T_j, T_t\mid X_t)$. Note that, if we denote with $\Delta = \alpha\, \gamma' + B$, then we can write by recursively expanding the linear Markovian expressions and using the fact that noise shocks are exogenous and jointly independent (see Appendix (ref) for details):
and we conclude that:
Note that the matrix $\Delta:=\alpha\gamma' + B$ captures the relationship between treatment $X_t$ and $X_{t-1}$ under the observational policy. Observe that under standard linear system theory $\|\Delta\|_{op}\leq 1$ is required for the system of states to be stable and not rapidly growing, which we would expect in many settings. Thus in practicy, under some stationarity of the states we would expect $\|\Delta\|_{op}$ to be small and therefore $\kappa$ to also be small. On the other hand if the data generating process is far from stationary and has long-range dependencies (i.e. a small change in the initial state at period $1$ can have a tremendous impact on the state at period $M$), then the estimation algorithm will incur an exponential dependence on the time horizon $m$.
Thus far we have assumed that we are working with $n$ independent time series, each of duration $m$. Though this is applicable to many settings where we have panel data with many units over time, in some other settings it is unreasonable to assume that we have many units over time, but rather that we have the same unit over a long period. In this case, we would want to do asymptotics as the number of periods grows. Our goal is still to estimate the dynamic treatment effects, i.e. the effect $\theta_{\kappa}$ of a treatment at period $t$ on an outcome in period $t+\kappa$, for $\kappa\in \{0, \ldots,m\}$) for some fixed look-ahead horizon $m$.
These quantities can allow us to evaluate the effect of counterfactual treatment policies on the discounted sum of the outcomes, i.e. $\sum_{t=0}^{\infty} \gamma^t Y_t$ for $\gamma<1$. We can write the counterfactual value function for any non-adaptive policy as: $V(\tau) = \sum_{t=0}^{\infty} \gamma^t \sum_{q\leq t} \theta_{t-q} \tau_{q}$. Assuming outcomes are bounded, the effect $\sum_{q\leq t} \theta_{t-q} \tau_{q}$ on any period $t$ can be at most some constant. Thus taking $m$ to be roughly $\log_{\gamma}(n)$, suffices to achieve a good approximation of the effect function $V(\tau)$, since the rewards vanish after that many periods, i.e. if we let: $V_{m}(\tau) = \sum_{t=0}^{m} \gamma^t \sum_{q\leq t} \theta_{t-q} \tau_{q}$, then observe that: $\|V_{m}(\tau) - V(\tau)\| \leq O(\gamma^{m})$. Thus after $m=\log_{1/\gamma}(n)$, we have that the approximation error is smaller than $1/\sqrt{n}$. Thus it suffices to learn the dynamic treatment effect parameters for a small number of steps. To account for this logarithmic growth, we will make the dependence on $m$ explicit in our theorems below.
For any $m$, we will estimate these parameters by splitting the time-series into sequential $B=n/m$ blocks of size $m$. Then we will treat each of these blocks roughly as independent observations and apply our dynamic DML algorithm to estimate parameters $\hat{\psi}_t$ and observe that under the markovian stationary nature of the DGP, we have that $\hat{\psi}_t = \hat{\theta}_{m-t}$. We denote the resulting estimate as $\hat{\theta}$. The main challenge in our proofs is dealing with the fact that these blocks are not independent but serially correlated. However, we can still apply techniques, such as martingale Bernstein concentration inequalities and martingale Central Limit Theorems to achieve the desired estimation rates.
The other important change that we need to make is in the way that we fit our nuisance estimates. To avoid using future samples to train models that will be used in prior samples (which would ruin the martingale structure), we instead propose a progressive nuisance estimation fitting approach, where at every period, all prior blocks are used to train the nuisance models and then they are evaluated on the next block. We present a formal description of this progressive splitting process in Algorithm (ref) and prove asymptotic normality of the resulting estimate.
\paragraph{Finite sample $\ell_2$-error.} We note that the proof of the latter theorem also provides a finite sample $\ell_2$-error bound. However we omit a separate such theorem for succinctness.
\paragraph{Conditions on $m$ and dependence on $\kappa$.} We provide some exposition on the conditions on $m$ as a function of the eigenvalues of the linear dynamical system in the case of a single treatment and a linear policy, e.g. $f(x,t)=\gamma'x$, as described in the corresponding remark at the end of Section (ref). Note that as long as the observational state transition matrix $\Delta=\alpha\gamma'+B$ has a small maximum eigenvalue $\ll 1$ and that $\|\alpha\|_2\|\gamma\|_2\ll 1$, then $\kappa \ll 1$, irrespective of the value of $m$. In other words, in this setting the correlations among the randomizations in the treatments are vanishing in an exponential manner and hence even if we estimate parameters in a long-chain, the estimation errors do not propagate in a manner that explodes exponentially with the length of the path. In that case the conditions on $m$ in Theorem (ref) simplify to:
Note that the number of nuisance functions also grows quadratically with $m$. Thus we should expect that $\epsilon_B(\hat{h})$ to also grow as $m^2$ times the convergence rate of each of the nuisance components. If each nuisance component estimate (denoted here $\hat{f}$) satisfies that $\frac{2}{B} \sum_{b=B/2}^{B} \mathbb{E}[(\hat{f}(Z_b) - f^*(Z_b))^2\mid {\mathcal F}_b]$ convergences at a rate of $\frac{1}{B^{1/2+\epsilon}}$, then a sufficient condition for all the latter properties is:
This allows for $m$ to grow polynomially with $n$, i.e. it suffices that $m = O(n^{\epsilon/(4+\epsilon) - \delta})$, for any $\delta>0$. Thus we can achieve very small approximation error if we are interested in a discounted reward, with discount $\gamma$, as in the beginning of a section, where with simply $m=\log_{1/\gamma}(n)$ we could achieve an error of $\gamma$. Thus we have that as long as the observational policy is such that the linear system is stable, estimation of long-term discounted rewards is feasible via Algorithm (ref).
We present a more formal treatment of the extension of our main algorithm to $g$-estimation of structural nested models in biostatistics robins1986new. Consider an arbitrary time-series process $\{X_t, T_t\}_{t=1}^{m}$, with $X_t\in {\mathcal X}_t$ and $T_t\in {\mathcal T}_t$. Let $Y$ denote some final outcome of interest. For any time $t$, let $\bar{X}_t=\{X_1,\ldots, X_t\}$ and $\bar{T}_t=\{T_1,\ldots, T_t\}$, denote the sequence of the variables up until time $t$ and similarly, let $\underline{X}_t = \{X_t, \ldots, X_m\}$ and $\underline{T}_t=\{T_t,\ldots, T_m\}$. We will also denote with $\bar{x}_t, \bar{\tau}_t, \underline{x}_t, \underline{\tau}_t$, corresponding realizations of the latter random sequences. Let $\pi=(\pi_1, \ldots, \pi_m)$ denote any dynamic policy, such that for each $t$, $\pi_t$ maps a history $\bar{x}_t, \bar{\tau}_{t-1}$ into a next period action $\tau_{t}$. For any such dynamic policy, let $Y^{(\pi)}$ denote the counterfactual outcome under policy $\pi$. For any static policy $\tau\in \times_{t=1}^m {\mathcal T}_t$, we will overload notation and let $Y^{(\tau)}$ denote the counterfactual outcome under this static treatment policy. Moreover, for any two policies (static or dynamic) we will be denoting with $(\bar{\pi}'_t, \underline{\pi}_{t+1})$, the policy that follows $\pi'$ up until time $t$ and then continues with policy $\pi$. We let $0\in {\mathcal T}_t$ denote a baseline policy value, which could be appropriately instantiated based on the context.
We assume that the data generating process satisfies the following sequential conditional randomization condition:
Identification of mean counterfactual outcomes $\mathbb{E}[Y^{(\pi)}]$ for a target policy of interest $\pi$ can be expressed in terms of the following conditional expectation functions:
which corresponds to the mean change in outcome if we go to all units which received treatment $\bar{\tau}_t$ up until time $t$ and had observed state history $\bar{x}_t$ and we remove their last treatment, while we subsequently always continue with the target policy $\pi$. These functions are known as the blip functions Chakraborty2013,Robins2004 and can be shown to be non-parametrically identifiable, assuming sequential conditional exogeneity and a sequential analogue of the positivity (aka overlap) assumption Robins2004.
Theorem 3.1 of Robins2004 combines a telescoping sum argument and the sequential randomization condition to express counterfactual outcomes in terms of blip functions. We restate this result here, adapting it to our notation and providing a proof for completeness:
Importantly, the conditioning set in Equation (ref) contains the observed $t$ periods treatment. Intuitively, each term $\rho_j$, removes from the outcome the blip effect of the observed action $T_j$ and adds the blip effect of the target action $\pi(\bar{X}_j, \bar{T}_{j-1})$. Lemma (ref), together with conditional sequential exogeneity also implies that the following set of moment restrictions must be satisfied (the following is an adaptation of Theorem 3.2 of Robins2004 to our notation and we include its proof for completeness).
To achieve parametric estimation rates for the quantities of interest, we will need to further make a semi-parametric assumption, i.e. that the blip functions take a low-dimensional parametric form:
Then we can identify $\psi^*$ by finding a parameter vector $\psi$ that satisfies the subset of the moment restrictions of the form:
Moreover, we can also subtract from $H_t(\psi)$, the conditional expectation $\mathbb{E}[H_t(\psi^*)\mid \bar{X}_t, \bar{T}_{t-1}]$ while maintaining the moment condition:
This is the doubly robust moment condition proposed by robins1994correcting,Robins2004, where it is shown that an estimator of $\psi$ based on this moment is correct if either the estimate of $\mathbb{E}[H_t(\psi^*)\mid \bar{X}_t, \bar{T}_{t-1}]$ or the estimate of $\mathbb{E}[\phi_t(\bar{X}_t, \bar{T}_t)\mid \bar{X}_t, \bar{T}_{t-1}]$ is correct. However, estimating $\mathbb{E}[H_t(\psi^*)\mid \bar{X}_t, \bar{T}_{t-1}]$, requires knowing the true parameter vector $\psi^*$, and since this is unknown, any feasible estimation approach must first construct preliminary estimates of $\psi^*$, which is computationally cumbersome and introduces another source of error. This issue has been discussed as one of the main points not to use the doubly robust correction in practice in $g$-estimation of structural nested models (see e.g. the discussion at the end of Section 6.1 of vansteelandt2014structural). For a binary treatment and when the target policy is the all-zero policy hernan2010causal (see Technical Point 21.5) note that $\mathbb{E}[H_t(\psi^*)\mid \bar{X}_t, \bar{T}_t]$ can be estimated by regressing the outcome of the population that received zero subsequent treatment on the history. However, such a population can be quite small in practice and can have severe co-variate imbalances compared to the overall population. Moreover, this approach only applies to the case of a binary treatment and a static target policy.
In this work, we show that we can achieve a Neyman orthogonal moment for identifying $\psi^*$, which is sufficient for robustness to biases stemming from machine learning models used to train the nuisance components, while avoiding the cumbersome part of estimating the nuisance $\mathbb{E}[H(\psi^*)\mid \bar{X}_t, \bar{T}_{t-1}]$. Moreover, our approach leads to a strongly convex loss for the parameters at each step of the recursive process, which is beneficial for finite sample guarantees and subsequently for generalizing it to linear models with parameter heterogeneity with respect to exogenous co-variates.
In particular, instead of subtracting $\mathbb{E}[H_t(\psi^*)\mid \bar{X}_t, \bar{T}_{t-1}]$, we subtract $\mathbb{E}[H_t(\psi)\mid \bar{X}_t, \bar{T}_{t-1}]$, i.e. the function that we subtract is dynamically dependent on the estimate of $\psi$. Even though this moment is not doubly robust, we will show that it remains locally robust LRSP, aka Neyman orthogonal chernozhukov2018double. To define our moment restrictions and connect them to the main results of the paper we first define several convenient random variables: for any $j\in [m]$, let
and for any $1\leq t\leq j$, let
Then, note that for any $\psi$:
Moreover, note that since the second term in $Q_t$ only depends on the conditioning set $\bar{X}_t, \bar{T}_{t-1}$, then:
Thus we conclude that the true parameter $\psi^*$ must be satisfying the moment restrictions:
These are exactly of the same form as the moment restrictions considered in the definition of the Dynamic DML algorithm and its analysis. Thus the results we presented so far directly extend to the estimation of structural nested mean models, with a linear parameterization of the blip functions, simply by using the different definition of the residual variables $\tilde{Y}_t$ and $\tilde{T}_{t,t}$ and letting $\psi_{t}=\theta_{m-t}$ in the definition of Algorithm (ref) and in Theorems (ref) and (ref).
Moreover, note that the nuisance models that are trained in the first stage have a larger conditioning set which includes all past history of states and treatments $\bar{X}_{t}, \bar{T}_{t-1}$. If we made further restrictions such that there was a "funnel state" $S_t$ at each period that summarizes the history and such that any dependence of the future to the past is going through that funnel state, then conditioning only on that state would have been sufficient. This is what we essentially did in the linear Markovian model. Moreover, the linear Markovian model with a static policy $\tau$, is a special case where the blip functions take the simple form: $\gamma_t(\bar{x}_t, \bar{\tau}_t)=\theta_{m-t}'\tau_t$ and are target policy independent. When the blip functions are target policy independent, then any target policy can be used to estimate the structural parameters $\psi^*$. In our main development, we essentially used the baseline zero policy as a target policy to estimate the structural parameters.
Thus our Dynamic DML algorithm extends to the estimation of the structural parameters in a structural nested mean model for any target dynamic policy $\pi$ and any user defined baseline policy $\bar{0}$ (referred to as a $(\pi,\bar{0})$-double regime structural nested mean model). It allows for the estimation of the nuisance functions with arbitrary machine learning algorithms, subject to a relatively slow mean-squared-error condition and reduces estimation to simple regression and classification oracles in the first stage, with only a simple linear system of equations in the second phase, which can also be solved in linear time in a recursive manner.
For completeness, we present the generalization of Algorithm (ref) to SNMMs in Algorithm (ref) and we re-state the main Theorems in this broader context.
We present here the generalized version of the asymptotic normality theorem for SNMMs and we defer the finite sample guarantee to the next section where we will analyze a more general setting of heterogeneous dynamic effects and present a finite sample bound. The proof is identical to that of Theorem (ref) and so we omit it.
Using Equation (ref) of Lemma (ref), we have that asymptotic normality and asymptotic linearity of the structural parameters, derived in Theorem (ref), also implies asymptotic normality of the plug-in off-policy value estimate for the target dynamic policy $\pi$.
We note that our moment condition that identifies each parameter $\psi_t$ is the derivative of a square loss and can be written as the solution to the square loss minimization problem defined in Equation (ref). This allows us to generalize the Dynamic DML to the case where we allow non-parametric heterogeneity in the parameters $\psi_t$, with respect to an exogenous fixed covariate vector of each sample, denoted as $X_0$, i.e. $\gamma_t(\bar{x}_t, \bar{\tau}_t)=\psi_t(x_0)'\phi(\bar{x}_t, \bar{\tau}_t)$, for a known feature map $\phi$ and unknown heterogeneous parameters $\psi$. Thus we can essentially generalize the $g$-estimation approach to SNMMs to allow for infinite or high dimensional parameters of the blip functions, as long as the input to these infinite dimensional parameters is fixed and not changing endogenously by the treatments (e.g. fixed characteristics of a unit). This can be achieved by simply minimizing recursively the square loss:
over arbitrary function spaces $\Psi_t$ or by using any other machine learning techniques that achieve small excess risk with respect to the latter square loss problem (e.g. regularized least squares, early stopping, etc). In the latter, we denoted with $\hat{h}$ an estimate of the vector of all nuisance functions (e.g. estimated in the first stage of the Dynamic DML Algorithm) and $\hat{\underline{\psi}}_{t+1}=(\hat{\psi}_{t+1}, \ldots, \hat{\psi}_m)$ the estimates of the target structural parameters constructed in previous iterations of the recursion.
Using techniques from the recently introduced orthogonal statistical learning framework foster2019orthogonal,chernozhukov2018plugin, we show in the appendix that this estimation method provides mean-squared-error guarantees on the recovered heterogeneous parameters $\hat{\psi}_t$, that are robust to errors in the nuisance functions. This heterogeneous extension can also be viewed as an analogue of the RLearner meta-learner algorithm nie2017quasi, generalized to the dynamic treatment regime setting. The formal description of the algorithm appears in Algorithm (ref), where we also allow for the function space over which we are optimizing the loss function $\ensuremath{{\cal L}}_{D,t}$, denoted as $\Psi_{t}^n$ to not necessarily be equal to $\Psi_t$ (which we know containts $\psi_t^*$) and to be changing with the sample size.
\paragraph{Norm notation.} To state our main results we will introduce some norm notation. For any vector valued function $\psi$, taking as input a random variable $X$ and having output in $\mathbb{R}^r$, we will denote with:
for any $u,v>1$. If $\psi$ is a parameter vector, then we will overload notation and let $\|\psi\|_{u, v}=\|\psi\|_{u}=\left(\sum_{j=1}^{r}\psi_j^u\right)^{1/u}$. If $u$ or $v$ equals $\infty$, then this would designate the sup norm, e.g. $\|\psi\|_{\infty, v} := \mathbb{E}[\max_{j\in r} \psi_j(X)^v]^{1/v}$ and $\|\psi\|_{u,\infty}=\sup_{x\in {\mathcal X}} \|\psi(x)\|_u$. For any $u,v$, we will denote with $\bar{u},\bar{v}$ the parameters that correspond to the dual norm, i.e. $1/u + 1/\bar{u}=1$ and similarly for $\bar{v}$. For any two functions $f, g$, taking as input random variables $X,Y$ we will use the shorthand notation: $\|f\circ g\|_{u,v} = \mathbb{E}[\|f(X)\|_u^v \cdot \|g(Y)\|_{u}^v]^{1/v}$. For an $n\times m$ matrix $A$, we will use the matrix norms:
and for $u=v=2$, we denote with $\|A\|_{op}=\|A\|_{2,2}$, the spectral or operator norm of $A$.
\paragraph{Algorithm agnostic robustness to nuisance.} We first prove a general bound on the estimation that is independent of the estimation process that is run in the second stage of Algorithm (ref). This lemma will be useful in multiple subsequent theorems in this and subsequent sections.
\paragraph{MSE for plug-in empirical risk minimization.} The main result of this section applies the latter lemma with $u=v=2$, $\sigma=\lambda$ and with $\psi_t^{*,n}=\operatorname*{arg\,inf}_{\psi_t\in \Psi_t^n} \|\psi_t - \psi_t^*\|_{2,2}$, so as to get a guarantee on the mean-squared-error of the heterogeneous structural parameters, as a function of the statistical complexity of the function spaces $\{\Psi_t^n\}_{t=1}^m$ and their bias with respect to the true heterogeneous effect parameters. Note that in this case, $c_{t,j} := \sup_{x_0\in {\mathcal X}_0} \|\mathbb{E}[C_{t,j} \mid X_0=x_0]\|_{op}$ and $M:=\sup_{t\in [m], x_0\in {\mathcal X}_0, \psi_t\in \Psi_t} \|\psi_t(x_0)\|_2$ and the convergence rates on pairs of nuisances will be with respect to the quantity $\|\|\Delta(\hat{f},\hat{g})\|_{op}\|_2$.
We will combine the conclusion of Lemma (ref) with a plug-in excess risk bound for the case when the estimate $\hat{\psi}$ is produced via running empirical risk minimization in the second stage of the Heterogeneous Dynamic DML algorithm, as described in Algorithm (ref). To state the theorem we will use the notion of the critical radius, which is a measure of statistical complexity of a function space. For any function space ${\mathcal F}$, with functions having range in $[-1,1]$, we consider the localized Rademacher complexity as:
We denote as the critical radius $\delta_n>0$ of ${\mathcal F}$, any solution to the inequality:
For each function space $\Psi_t$, we denote with $\Psi_{t,i}$ the marginal function space corresponding to the $i$-th coordinate output of the functions in $\Psi_t$. Moreover, we denote with $\Psi_{t,i}-\psi_{t,i}^*=\{\psi_{t,i} - \psi_{t,i}^*: \psi_{t,i}\in \Psi_{t,i}\}$. Finally, we define the star hull of a function space as: $\text{star}({\mathcal F}) := \{\tau f: f\in {\mathcal F}, \tau\in [0,1]\}$. The critical radius is a well-established concept in modern statistical learning theory and has been characterized for many function spaces. Moreover, for many function spaces it yields minimax optimal statistical learning rates. See wainwright_2019 for an overview.
Analogous results hold with high probability and exponential tail, if we make such exponential tail assumptions also on the guarantees provided by the nuisance functions. We omit them for succinctness.
\paragraph{Partial double robustness.} Note that most conditions on the nuisance functions have a doubly robust flavor, i.e. we need the product of two different nuisance function errors to be small. In particular, observe that by Jensen's inequality, for any $f, g$, we have that:\footnote{Since: $\mathbb{E}\left[\mathbb{E}[\|\hat{f}(Z_f) - f^*(Z_f)\|_2 \|\hat{g}(Z_g) - g^*(Z_g)\|_2\mid X_0]^2\right] \leq \mathbb{E}\left[\|\hat{f}(Z_f) - f^*(Z_f)\|_2^2 \|\hat{g}(Z_g) - g^*(Z_g)\|_2^2\right]$}
while when $X_0$ is the empty set (i.e. no heterogeneity), then we have:
in which case the conditions on the nuisance estimates boil down to the same as those in Theorem (ref) in the expository section. Thus we need that either one or the other nuisance to be modeled and estimated accurately. The only exception is the functions $p_{t,t}$, which also need to satisfy that: $\|\hat{p}_{t,t} - p_{t,t}^*\|_{2,4}^4 = o_p(\delta_{n/2}^2)$. Thus one step ahead treatment propensities, need to be more accurately estimated than the remainder of the nuisance functions.
\paragraph{Alternative norm bounds.} We note that when the $\ell_{2,2}$ norm of $\hat{\psi}_t-\psi_t^*\in \Psi_t$ is lower bounded by some fraction of its $\ell_{2,\infty}$ norm, i.e. the sup norm over $X_0$ (e.g. if $\psi_t$ is a linear class and $\mathbb{E}[X_0\, X_0']\succeq \mu I$, a special case of which is when there is no heterogeneity), then invoking Lemma (ref) in the proof of Theorem (ref) with $v=\infty$ instead of $v=2$, we can get a result of the form:
subject to a weaker $\ell_{1}$ norm convergence for the products of the nuisance components:
By applying a Holder inequality,\footnote{Since: $\left\|\|\Delta(\hat{p}_{t,t}, \hat{p}_{j,t})\|_{op}\right\|_{1}\leq \mathbb{E}\left[\|\hat{f}(Z_f) - f^*(Z_f)\|_2 \|\hat{g}(Z_g) - g^*(Z_g)\|_2\right]\leq \|\hat{f}-f^*\|_{2,2}\, \|\hat{g}-g^*\|_{2,2}$} the latter is satisfied whenever:
Recovering again qualitatively the same norm convergence conditions as in Theorem (ref). When $\Psi_t$ is a parametric class with a bounded domain of parameters, then $\delta_n=O\left(n^{-1/2}\right)$, in which case, the latter requirement is satisfied if each nuisance function $f\in \{\hat{p}_{j,t}, \hat{q}_t\}_{1\leq t\leq j\leq m}$, satisfies that $\|\hat{f}-f^*\|_{2,2}=o_p(n^{-1/4})$, recovering the typical conditions in the Neyman orthogonality literature with parametric target estimands chernozhukov2018double.
For more general function classes $\Psi_t$, where the $\ell_{2,2}$ norm of $\hat{\psi}_t$ is not related to its $\ell_{2,\infty}$ norm and without any further restrictions on the correlations of errors among the nuisance components in each of the product term conditions in Theorem (ref), then with a simple Holder inequality applied to the nuisance constraints, as in Equation (ref) we can still derive a sufficient condition of the same form as in Equation (ref), but with the $\ell_{2,2}$ norms replaced by the slightly stronger $\ell_{2,4}$ norms. Hence, it suffices that each nuisance function satisfies $\|\hat{f}-f^*\|_{2,4}=o_p(\sqrt{\delta_n})$, which for parametric classes would be $o_p(n^{-1/4})$. If errors in the nuisance components are un-correlated then an $\ell_{2,2}$ norm convergence would suffice for most nuisances, since each product of nuisance error terms can be upper bounded as:
with the only $\ell_{2,4}$ norm required for the nuisances $\{p_{t,t}\}_{t=1}^m$, i.e. the one step ahead observational propensity models.
\paragraph{Uniform consistency.} Looking at Lemma (ref), we note that if the $\ell_{2,2}$ norm of $\hat{\psi}_t-\psi_t^{*,n}$ is lower bounded by some fraction of its $\ell_{2,\infty}$, then a uniform consistency result can be derived. Since both of these functions lie in $\Psi_t^n$, if we choose the classes $\Psi_t^n$, such that their elements satisfy this property and such that as $n$ grows, the function spaces $\Psi_t^n$ uniformly approximate $\Psi_t$, then we can achieve a uniform consistency theorem. Moreover, in this case we can invoke the lemma with $v=\infty$, which would only require the weaker norm guarantees on the nuisances.
The results we have discussed so far assume that the linear feature map that parameterizes the blip functions is low dimensional, i.e. $r\ll n$. Observe for instance that Theorem (ref) is vacuous when $r=\Omega(n)$. In this section we examine the case where $r\gg n$ and provide guarantees under sparsity conditions on the true structural parameters. For simplicity, we will not consider non-parametric heterogeneity of the sparse coefficients with respect to some initial state $X_0$, i.e. we consider blip functions of the form: $\gamma_t(\bar{X}_t, \bar{T}_t;\psi_t)=\psi_t'\phi(\bar{X}_t,\bar{T}_t)$, with $\psi_t\in \mathbb{R}^r$ and $r\gg n$. However, we note that here the high-dimensionality of the feature map already offers a lot of modelling flexibility and one could encode heterogeneity of the blip effect through the feature map.
Apart from the explicit dependence on $r$, when the feature map is high dimensional then the $\ell_{2,2}$ norm of the errors of the nuisance functions in Theorem (ref) can accumulate across their $r$-dimensional components. Instead, we would ideally only require a bound on the maximum error across the $r$ dimensions of each nuisance function. Then an exponential tail bound on the MSE of each coordinate would imply a bound on the maximum that scales only logarithmically with $r$. To achieve this we can invoke Lemma (ref) with $u=1, v=\infty$ and $\psi^{*,n}=\psi^*$ and $X_0$ an empty set. Then we get the following recursive bound on the MSE of the heterogeneous structural parameters:
where $c_{t,j} := \|\mathbb{E}[C_{t,j}]\|_{\infty}$ and $M:=\sup_{t\in [m], \psi_t\in \Psi_t} \|\psi_t\|_1$ and we used the short-hand norm notation $\|A\|_{\infty}=\max_{i,j} |A_{i,j}| = \|A\|_{\infty,1}$.
However, we see that we incur a dependency on the $\ell_{1}$ norm of the error $\hat{\psi}_t - \psi_t^*$. Thus we need to be able to relate the $\ell_{1}$ with the $\ell_{2}$ norm of the error of our estimate. This is a restricted cone property on our estimate. We will thus invoke sparsity assumptions on the true parameter $\psi_t^*$ and augment the empirical risk minimization step with an $\ell_1$ penalty on $\psi_t$. This enforces the estimate to be primarily supported on the $s$ relevant dimensions. Within such a restricted cone the $\ell_1$ and the $\ell_2$ norm are equivalent with a constant that depends only on the sparsity level and not the ambient dimension.
Furthermore, the explicit dependence in Theorem (ref) also stems from the fact that we invoked a multi-dimensional contraction inequality, across the $r$ output dimensions of $\psi_t$. However, observe that the loss function depends on each $\psi_t$ only through single indices of the form $\psi_t'\tilde{T}_{j,t}$. Thus we could instead control the critical radius of these $m$ index spaces at each iteration of the recursion, which would avoid the explicit dependence on $r$. Together these insights yield the following theorem.
We consider data drawn from the DGP presented in Equation (ref), with a linear observational policy:
with $X_0, T_0 = 0$ and $\epsilon_t, \zeta_t, \eta_t$ standard normal r.v.'s. (recall that $d$ is the number of treatments and $p$ the number of state variables). We consider the instance where: $A_{ij} = .5$, for all $i\in [p]$, $j\in [d]$, $B = .5\, I_p$, $C = .2\, I_d$, $D[:, 1:2] = .4$, $D[:, 3:p]= 0$, $\mu[1:2] = .8$. We consider settings where the effect is constant, i.e. $\theta_0\in \mathbb{R}^d$ or heterogeneous, where:
for some known low dimensional subset $S$ of the states.
We compare the dynamic DML to several benchmarks. The results are presented in Figures (ref), comparing the estimates of the dynamic DML algorithm to a number of other alternatives on a single instance. They fall into two categories. In the “static” set of approaches, each of the contemporaneous and lag effects is estimated one at a time, either by direct regression or (static) DML. So for example, to estimate the one period lag effect $\theta_1$, we would regress $Y_t$ on $T_{t-1}$, with controls $X$. We consider direct regression with no controls (“no-ctrls”), direct and DML with controls from the inital period (i.e $X_{t-1}$, “init-ctrls” and “init-ctrls-dml”) and direct and DML with controls from the same period as the outcome (i.e $X_t$, “fin-ctrls” and “fin-ctrls-dml”). As an alternative to all of these, we try a “direct” dynamic approach, where initially we estimate $\theta_0$ using a lasso regression of $Y_t$ with all the controls and past treatments, and return the coefficient on $T_t$, and then “peel” off the estimated effect as in the main text before running another Lasso regression of $Y_t - \theta_0 T_t$ on $T_{t-1}$ to get the first lag effect etc. So this approach incorporates the peeling effect but doesn't do any orthogonalization.
The point estimate for $\theta_0$, $\theta_1$ and $\theta_2$ are depicted in the three panels of Figure (ref) and the error bars correspond to the constructed confidence interval. For all three, the dynamic DML is relatively close to the truth and the confidence interval contains the truth. The remaining approaches are not, although for the contemporaneous effect the approaches with final period controls have similar performance - it is really in the lagged effects that the differences become most apparent. Subsequently we run multiple experiments to evaluate the performance of DynamicDML. In each setting, we run $1000$ Monte Carlo experiments, where each experiment draws $N = 500$ samples from the above DGP and then estimated the effects and lag effects based on our Dynamic DML algorithm. Figure (ref) considers the case of two treatments, and shows that the algorithm performs well in terms of giving reasonable coverage guarantees - for a nominal 95% coverage, actual coverage varies from 91% to 94.5%. The right panel shows that the average estimates are close to the truth. In Figure (ref) we repeat the experiments with $N=2000$, and actual coverage is now tightly in the range 94% to 95%, and the average estimates remain relatively unbiased. We also find qualitatively similar performance for the case of estimating heterogeneous treatment effects and policy values (see Figures(ref) and (ref)).