Dynamic covariate balancing: estimating treatment effects over time with potential local projections

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.

77,789 characters · 14 sections · 92 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractThis paper studies the estimation and inference of treatment effects in panel data settings when treatments change dynamically over time. We propose a balancing method that allows for (i) treatments to be assigned dynamically over time based on high-dimensional covariates, past outcomes, and treatments; (ii) outcomes and time-varying covariates to depend on the trajectory of all past treatments; (iii) heterogeneity of treatment effects. Our approach recursively projects potential outcomes' expectations on past histories. It then controls the bias arising from the non-experimental and sequential nature of this setting by balancing dynamically observable characteristics over time. We establish inferential guarantees of the proposed method even when the number of observable characteristics significantly exceeds the sample size. We study numerical properties of the estimator and illustrate the benefits of the procedure in an empirical application.

keywordsCausal Inference, High Dimensions, Treatment Effects, Panel Data.

Introduction

\spacingset{1.7} Researchers collect a panel of $n$ independent observations observed over a finite number of $T$ periods in an observational study. The dataset encompasses time-varying covariates, outcomes, and time-varying treatments. The primary objective is to conduct inference on the average effect of exposure to different treatment histories, such as the effect of being treated for a certain number of periods.

We consider a setting where treatments change dynamically over time, and potential outcomes, covariates and treatment may depend on past histories. Two alternative procedures can be considered in this setting. First, researchers may consider explicitly modeling how treatment effects propagate over each period through time-varying covariates and intermediate outcomes. This approach is prone to large estimation error and misspecification in high-dimensions: it requires modeling outcomes and each time-varying covariate as a function of all past covariates, outcomes, and treatment assignments. A second approach is to use inverse-probability weighting estimators for estimation and inference tchetgen2012semiparametric, vansteelandt2014structural . However, classical semi-parametric estimators are prone to instability in the estimated propensity score. There are two main reasons. First of all, the propensity score defines the joint probability of the entire treatment history and can be close to zero for moderately long treatment histories. Additionally, the propensity score can be misspecified in observational studies.

This is a common problem both in social sciences and bio-statistics. For example, in a survey of all articles in 2021 top-5 economics journals, more than $20\%$ of studies with time-varying treatments exhibit treatment dynamics.\footnote{This is based on the authors' calculation. Top-5 economics journals are American Economic Review, Econometrica, Journal of Political Economy, Quarterly Journal of Economics, Review of Economic Studies.} On the other hand, typical approaches in economics and related disciplines often employ, Difference-in-Differences designs. If units can dynamically choose treatments in response to their previous outcomes (or treatments), this will lead to violations of the parallel trends assumption required by such designs ghanem2022selection, marx2022parallel. We, therefore, introduce an approach that is valid when treatment decisions at period $t$ can depend on the history of outcomes and treatments prior to $t$. The second challenge is that treatment dynamics are difficult to estimate. Individuals may select into treatment arbitrarily based on high-dimensional covariates, outcomes, and treatments, e.g., when maximizing future expected utilities heckman2007dynamic. This motivates a method that does not impose modeling assumptions on selection into treatment mechanisms (i.e., propensity score).

This paper studies the estimation and inference of the effects of treatment histories when potential outcomes (and covariates) depend on present and past treatments. Individuals dynamically select into treatment based on past (time-varying) covariates, outcomes, and treatments. There are no unobserved confounders after controlling for high-dimensional past characteristics ding2019bracketing. Researchers remain agnostic on the propensity score.

We leverage a model on the potential outcomes' conditional expectations as an (approximately) linear function of previous potential outcomes and (high-dimensional) covariates in each period. Our model is motivated by local projection frameworks jorda2005estimation, montiel2021local. Local projections impose a (linear) model on observed outcomes conditional on each period observables and do not require estimating how each time-varying covariate changes in response to treatments -- which would be prone to large estimation error in high dimensions. However, different from standard local projections, our model is imposed on expected potential instead of observed outcomes. This difference is important here because of treatments' serial correlation and selection into treatment based on past outcomes and covariates: a model on realized outcomes imposes restrictions on the distribution of the treatment assignments, whereas a potential outcome model does not. Building on the literature on marginal structural models robins2000marginal, we identify the parameters of interest by recursively projecting outcomes' conditional expectations over past histories, allowing for dynamic selection into treatment.

Our estimation method, Dynamic Covariate Balancing (DCB), estimates the parameters of the model by using recursive penalized projections through lasso hastie2015statistical. It then reweights observations to guarantee balance between treated and control units. Balancing covariates is intuitive and common in practice: in cross-sectional studies, treatment, and control units are comparable when the two groups have similar characteristics imai2014covariate, li2018balancing, hainmueller2012entropy. We generalize covariate balancing in the absence of dynamics of zubizarreta2015stable, athey2018approximate, ben2018augmented, hirshberg2017augmented to a dynamic setting. We show that balancing with potential local projections corresponds to constructing weights sequentially in time by first balancing treated and control units' covariates in the first period and then balancing histories in the next periods reweighted by the weights obtained in the previous period. The estimated balancing weights solve a sequence of quadratic programs to minimize the weights' variance.

Our estimation procedure guarantees a vanishing bias of order faster than $n^{-1/2}$ and a parametric rate of convergence of the estimated treatment effect in high-dimensional settings. In addition, the optimization problem over the set of balancing weights admits a feasible solution, with the true propensity score being one such solution (and without requiring knowledge of it). This result highlights the benefits of balancing over propensity score reweighting here: the proposed balancing weights have a smaller variance than inverse probability weights and -- by leveraging an (approximate) high-dimensional linear outcome model -- do not require the correct specification of the propensity score.\footnote{Typical methods in high dimensions require conditions on the product of the rates of estimators for the propensity score and coefficients of the linear model to be faster than $n^{-1/4}$, and also require consistent estimation of both the outcome model and propensity score model athey2017efficient. Compared to estimating the propensity score with a semi-parametric model, our guarantees do not depend on the estimation error of the propensity score (only require that the estimation error of the coefficients is $o(n^{-1/4})$), by leveraging the high dimensional linear outcome model. } This is an advantage especially in dynamic settings: the propensity score defines the joint probability that units are assigned to a given treatment history, and therefore inverse probability weights can exhibit large variance in finite sample (see e.g., Figure (ref)). Finally, we provide guarantees for inference. Relative to cross-sectional studies, our dynamic structure necessitates novel considerations for identification, balancing, and derivations, that require analyzing joint distributions of correlated residuals from sequential projections.

We illustrate our method in an empirical application using data from acemoglu2019democracy on studying the effects of democracy on economic growth. Here, the authors assume a dynamic selection model. Whereas effects are in magnitude and sign consistent with acemoglu2019democracy, we show that standard local projections and acemoglu2019democracy's linear regression lead to significantly smaller point estimates compared to our approach. We also show that (A)IPW methods lead to a more substantial imbalance (and bias) compared to DCB due to the instability of the propensity score in both high and low-dimensions.

Related Literature

The goal of this paper is to conduct inference on dynamic treatment effects, while being robust to the misspecification of the propensity score. To achieve this goal, we leverage a high dimensional linear model and derive the first dynamic balancing equations within the local projection model proposed in this paper.

In the econometrics and statistics literature, imbens2019panel propose balancing assuming no treatment dynamics, whereas here, treatment dynamics require different (and novel) balancing conditions. In the context of dynamics, different from imai2015robust, who estimate a single set of balancing weights over all possible combinations of time periods and covariates, here the number of moment conditions grows linearly with $T$ and not exponentially. Unlike zhou2018residual, who extend entropy balancing of hainmueller2012entropy to dynamic settings, and li2024toward who propose a single set of balancing by regressing each covariate on past information, we do not estimate one model for each covariate in the past (which can be prone to large estimation error in high dimensions). DCB explicitly characterizes the high-dimensional model's bias in a dynamic setting to avoid overly conservative moment conditions, while kallus2018optimal design conservative balancing conditions for the worst-case bias. Different from yiu2018covariate, we do not require estimating the propensity score. Our insight with respect to all these references (in low and high dimensions) is that with a linear model and sequential weights, balancing reduces to few and novel dynamic restrictions. This insight is even more relevant with high-dimensional covariates, which none of these references study with dynamics.

Compared to cross-sectional studies, we generalize balancing in athey2018approximate, ben2018augmented, and consider an arbitrary class of weights. Therefore, our residual balancing procedure does not reduce to linear estimators as in settings with linear balancing weights bruns2023augmented, and our analysis differs from cross-sectional studies with low dimensions in wang2020minimal.

More broadly, this paper connects to the literature on DiD, local projections, and dynamic treatments. Different from the literature on DiD rambachan2023more, de2022difference, callaway2019difference, abraham2018estimating, athey2022design, caetano2022difference or subsequent work on potential projections with DiD dube2023local, here we allow for dynamic treatment regimes. This literature imposes the parallel trends assumption, violated with dynamic treatments marx2022parallel, ghanem2022selection. Different from the time-series literature montiel2021local, stock2018identification, rambachan2019nonparametric, this paper uses information from panel data and allows for arbitrary dependence of outcomes, covariates, and treatment assignments over time.

References in bio-statistics include robins2000marginal, hernan2001marginal, boruvka2018assessing, blackwell2013framework, bang2005doubly vansteelandt2014structural. bojinov2020panel study IPW estimators from a design-based perspective. Doubly robust estimators for dynamic treatments have been studied by nie2021learning, zhang2013robust, jiang2015doubly, tchetgen2012semiparametric, babino2019multiple. Here, we focus on studying the effect of a given treatment path, as in robins2000marginal or blackwell2013framework, different from and complementary to studying optimal policies (e.g., murphy2003optimal, nie2021learning).

Specifically, studies with high-dimensional panels require correct specification of the propensity score lewis2020double, zhu2017high, shi2018high, bodoryevaluating, belloni2016inference, chernozhukov2017orthogonal, or impose homogeneous treatment effects high_dim_IRF, kock2015inference. (lewis2020double also illustrate bounds on misspecification). More generally, prior works that formally study properties of dynamic doubly-robust methods in high dimensions require product of rates conditions for the estimated propensity score and conditional mean function, and consistent estimation of both; see for example follow up work by bradic2021high who provide tight rates of convergence of dynamic AIPW. Different from above, our framework does not require consistent estimation of the propensity score.

Finally, in both works subsequent to the first version of this paper, chernozhukov2022automatic generalize the use of riesz representers with arbitrary non-linear outcome models, and zhang2021dynamic study doubly robustness to model misspecification through moment restrictions. Different from these references, here we do not require conditions on the balancing weights motivated by our goal of allowing for inference with a possibly completely misspecified propensity score function, whereas chernozhukov2022automatic and zhang2021dynamic require functional form restrictions on the balancing weights to obtain a product of rates conditions. Our focus on the high-dimensional linear model (which we view as a linear approximation to conditional expectations in high dimensions) is motivated by its large use in applications.

Dynamics and potential local projections

Setup

We start with the analysis of two time periods, deferring multiple periods to Section (ref). We observe a panel with $n$ $i.i.d.$ copies of $ \Big( X_{i,1}, D_{i, 1}, Y_{i, 1}, X_{i, 2}, D_{i, 2}, Y_{i, 2} \Big)$, each distributed according to $\mathcal{P}$. Here $ D_{i,1}, D_{i,2} \in \{0,1\}$ denote binary treatments at time $t = 1,t = 2$, respectively, $X_{i,t}, Y_{i,t}$ denote covariates and the outcome at time $t$. We allow for any nonstationarity and dependencies that may occur over time within each unit. When indices are not specified, such as in $D_t$, this refers to the collective observations for all $n$ units.

figure[figure omitted — 3,123 chars of source]

We consider potential outcomes that are functions of the entire treatment history with $Y_{i,2}(d_1, d_2)$ denoting the potential outcome at time $t = 2$, under treatment $d_1$ in the first and $d_2$ in the second period. Our goal is to conduct inference on the estimand(s) $$

aligned\mathrm{ATE}(d_{1:2}, d_{1:2}') = \mu_2(d_1, d_2) - \mu_2(d_1', d_2'), \quad \mu_2(d_1, d_2) = \mathbb{E}\Big[Y_{i,2}(d_1, d_2)\Big],

$$ for given treatment histories $(d_1, d_2), (d_1', d_2')$. For example, researchers may be interested in estimating $ \mathrm{ATE}((1,1), (0,0)), $ which denotes the \textit{total} effect of treating an individual for two consecutive periods \citep{athey2022design}; or the \textit{direct} effect $\mathrm{ATE}((1, 0), (0,0))$. Figure (ref) shows that the overall treatment effects capture the direct effect of the treatment on the outcomes and the indirect effect. For longer histories, one could also consider weighted combinations of relevant treatment effects, omitted for brevity (see Section (ref)).

figure[figure omitted — 2,368 chars of source]

Dynamic treatment assignments

Treatment histories can impact both outcomes and covariates at intermediate stages. Let $Y_{i,1}(d_1, d_2)$ represent the intermediate potential outcome and $X_{i,2}(d_1, d_2)$ represent the potential covariates following a sequence of $d_1$ then $d_2$. Here, $X_{i,1}$ refers to the baseline covariates.

assFor $d_1 \in \{0,1\}$, let $Y_{i,1}(d_1, 1) = Y_{i,1}(d_1, 0)$, $X_{i,2}(d_1,1) = X_{i,2}(d_1,0)$.

Assumption (ref) is a no-anticipation restriction: (i) intermediate potential outcomes only depend on past but not future treatments; (ii) the treatment status at $t = 2$ has no contemporaneous effect on covariates.

Assumption (ref) allows for anticipatory effects governed by expectations (e.g., individuals may choose treatments based on expected future utilities), but not on the future treatment realizations athey2022design.

exmp[Observed outcomes] Consider a dynamic model of the form $$ \small \begin{aligned} Y_{i,2} = g_2\Big(Y_{i,1}, X_{i,1}, X_{i,2}, D_{i,1}, D_{i,2}, \varepsilon_{i,2}\Big), \quad Y_{i,1} = g_1\Big(X_{i,1}, D_{i,1}, \varepsilon_{i,1}\Big), \quad X_{i,2} = g_0\Big(X_{i,1}, D_{i,1}, \varepsilon_{i, X}\Big) \end{aligned} $$ for some arbitrary functions $g_2(\cdot), g_1(\cdot), g_0(\cdot)$ and unobservables $(\varepsilon_{i,2}, \varepsilon_{i,1}, \varepsilon_{i,X})$, with $ \varepsilon_{i,2} \perp D_{i,2} | Y_{i,1}, X_{i,1}, X_{i,2}, D_{i,1}, $ and $ (\varepsilon_{i,X}, \varepsilon_{i,1}) \perp D_{i,1} | X_{i,1}. $ We can write $$ \small \begin{aligned} Y_{i,2}(d_1, d_2) = g_2\Big(Y_{i,1}(d_1), X_{i,1}, X_{i,2}(d_1), d_1, d_2, \varepsilon_{i,2}\Big), \end{aligned} $$ where $Y_{i,1}(d_1) = g_1\Big(X_{i,1}, d_1, \varepsilon_{i,1}\Big), X_{i,2} = g_0\Big(X_{i,1}, d_1, \varepsilon_{i, X}\Big)$. Since $g_1(\cdot), g_0(\cdot)$ are not functions of $d_2$, Assumption (ref) holds. (Assumption (ref) below will impose restrictions on $\mathbb{E}[g_2(\cdot)]$.) \qed

In the rest of our discussion, we index potential outcomes and covariates by past treatment history under Assumption (ref). We define $ H_{i,2} = \Big[D_{i,1}, X_{i,1}, X_{i,2}, Y_{i,1}\Big], $ the vector of past treatment assignments, covariates, and outcomes in the previous period. We refer to $ H_{i,2}(d_1) = \Big[d_1, X_{i,1}, X_{i,2}(d_1), Y_{i,1}(d_1)\Big] $ as the potential history under treatment status $d_1$ in the first period. Here, $H_{i,2}$ can include interaction terms, omitted for brevity.

ass[Sequential Ignorability] Assume that for all $(d_1, d_2) \in \{0,1\}^2$ , $$ \small \begin{aligned} (A) \quad &Y_{i,2}(d_1, d_2) \perp D_{i,2} \Big | D_{i,1}, X_{i,1}, X_{i,2}, Y_{i,1}, \quad (B) &\Big(Y_{i,2}(d_1, d_2), H_{i,2}(d_1)\Big) \perp D_{i,1} \Big | X_{i,1}, \end{aligned} $$

Sequential ignorability states that treatment in the first period is unconfounded conditional on baseline covariates, and the treatment in the second period is unconfounded conditional on all observable characteristics at $t= 2$. It assumes no unobserved factors after controlling for high dimensional observable characteristics and arbitrary past information. Note that we could also state (A), conditioning on $D_{i,1} = d_1$ and potential history $H_{i,1}(d_1)$.

In Example (ref), Assumption (ref) holds if $ D_{i,2} \perp \varepsilon_{i,2} \Big| D_{1, i}, X_{i,1}, X_{i,2}, Y_{i,1}, \quad D_{i, 1} \perp (\varepsilon_{i,1}, \varepsilon_{i,2}) \Big| X_{i,1}. $

Potential local projections

Following in spirit, jorda2005estimation, we approximate the expectation of potential outcomes as linear functions of (high-dimensional) past characteristics. Different from jorda2005estimation, linearity is imposed on expected potential instead of realized outcomes. Modeling potential outcomes directly avoids functional form restrictions on the treatment assignment mechanism.

assFor some $\beta_{d_1, d_2}^{(1)} \in \mathbb{R}^{p_1}, \beta_{d_1, d_2}^{(2)} \in \mathbb{R}^{p_2}$ $$ \small \begin{aligned} & \mathbb{E}\Big[Y_{i,2}(d_1, d_2)\Big| X_{i,1} = x_1\Big] = x_1\beta_{d_1, d_2}^{(1)}, \quad \\ & \mathbb{E}\Big[Y_{i,2}(d_1, d_2) \Big| X_{i,1} = x_1, X_{i,2} = x_2, Y_{i,1} = y_1, D_{i,1} = d_1\Big] = \Big[d_1, x_1, x_2, y_1\Big] \beta_{d_1, d_2}^{(2)}. \end{aligned} $$

Assumption (ref) allows for heterogeneity in $(d_1, d_2)$, and the dimensions $p_1, p_2$ can grow with $n$ (because of additional covariates and/or covariates transformations). As for MSMs robins2000marginal, Assumption (ref) (i) does not require estimating a structural model for each time-varying-covariate, that would be prone to large estimation error in high dimensions; and (ii) it is agnostic on the treatment assignment mechanism because the model is imposed on potential outcomes. Coefficients can vary with time in the model.

lem[Identification] Let Assumptions (ref), (ref), (ref) hold. Then $$ \small \begin{aligned} & \mathbb{E}\Big[Y_{i,2} \Big| H_{i,2}, D_{i,2} = d_2, D_{i,1} = d_1\Big] = \mathbb{E}\Big[Y_{i,2}(d_1, d_2) \Big| H_{i,2}, D_{i,1} = d_1\Big] = H_{i,2}(d_1) \beta_{d_1, d_2}^{(2)} \\ & \mathbb{E}\Big[\mathbb{E}\Big[Y_{i,2} \Big| H_{i,2}, D_{i,2} = d_2, D_{i,1} = d_1\Big] \Big| X_{i,1}, D_{i,1} = d_1\Big] = \mathbb{E}\Big[Y_{i,2}(d_1, d_2) \Big| X_{i,1} \Big] = X_{i,1} \beta_{d_1, d_2}^{(1)}. \end{aligned} $$

The proof is in Appendix (ref). Lemma (ref) builds on results in the literature on marginal structural models robins2000marginal, bang2005doubly, tran2019double, kallus2020double, where, here, we make a connection between marginal structural models and local projections in economics as a contribution of independent interest. Lemma (ref) motivates a recursive estimation strategy discussed in Section (ref).

exmp[Linear Model] Let $X_{i,1}, X_{i,2}$ also contain an intercept. Let $\mathbb{E}\Big[Y_{i,1}(d_1) \Big| X_{i,1}\Big] = X_{i,1} \alpha_{d_1}$, $\mathbb{E}\Big[ X_{i,2}(d_1) \Big|X_{i,1}\Big] = W_{d_1} X_{i,1}$, $\mathbb{E}\Big[Y_{i,2}(d_1, d_1) \Big| X_{i,1}, X_{i,2}, Y_{i,1}, D_{i,1} = d_1\Big] = \\ \Big(X_{i,1}, X_{i,2}(d_1), Y_{i,1}(d_1)\Big) \beta_{d_1, d_2}^{(2)}, $ for some arbitrary parameters $\alpha_{d_1} \in \mathbb{R}^{p_1}$ and $\beta_{d_1, d_2}^{(2)} \in \mathbb{R}^{p_2}$; $W_{d_1}, V_{d_1} $ denote unknown matrices in $ \mathbb{R}^{p_2 \times p_1}$ . The model satisfies Assumption (ref). \qed

rem[Linearity in high-dimensions as an approximation to the true model] In the same spirit of belloni2014inference, our results also directly extend to the case where we relax Assumption (ref) and assume only approximate linearity up to an order $\mathcal{O}_p(r_p), $ where $r_p$ is an arbitrary sequence which depends on $p$ with $r_p = o(n^{-1/2})$. This setting embeds empirical applications where many covariates (and their transformation) can approximate the conditional mean function as linear. As we further discuss in Section (ref), we consider a high-dimensional setting where we will only require weak conditions on the estimated coefficients $||\hat{\beta}_{d_1, d_2}^{(t)} - \beta_{d_1, d_2}^{(t)}||_1 = O_p(n^{-1/4})$ with unbounded covariates and $||\hat{\beta}_{d_1, d_2}^{(t)} - \beta_{d_1, d_2}^{(t)}||_1 = o_p(1/\log(n))$ with bounded covariates, common for standard high-dimensional estimators. \qed

rem[Comparison with standard local projections and DiD] Appendix (ref) presents an extensive discussion and comparison of Lemma (ref) with standard local projections and DiD common in economics that, as we show, would return biased estimates in a dynamic context. The reason is because standard local projections impose a linear model on the observed outcomes $Y_{i,T}$ instead of potential outcomes, which therefore would also depend on the distribution of treatment assignments.

Estimation with dynamic balancing

Estimation with two periods

This section studies estimation. We defer to Section (ref) a complete guide for practice, including discussion about the model, tuning parameters, and complexity. Appendix (ref) presents a more detailed description of each step to construct the estimator.

Consider a two periods setting first, with $T = 2$. The estimator for $\mu(d_1,d_2)$ (and symmetrically for $\mu(d_1', d_2')$) proceeds in the following steps:

itemize• Estimation of the coefficients: We estimate the coefficients $\hat{\beta}_{d_1, d_2}^{(2)}$ by regressing $Y_2$ onto $H_2$ (controlling for $D_2 = d_2$). Following Lemma (ref), we then regress $H_2 \hat{\beta}_{d_1,d_2}^{(2)}$ onto $X_1$ (controlling for $D_1 = d_1$) to estimate $\hat{\beta}_{d_1,d_2}^{(1)}$. These regression may allow for high-dimensional coefficients as for lasso. The full algorithm is in Algorithm 2. • Sequential estimation via regression adjustments: Because linearity may only hold as we control for high-dimensional covariates, we cannot directly use our predictions $H_2 \hat{\beta}^{(2)}$ or $X_1 \hat{\beta}^{(1)}$ for valid causal inference. We instead must guarantee a vanishing high-dimensional bias through reweighting. For weights in each period $\hat{\gamma}_2(d_1,d_2)$ and $\hat{\gamma}_1(d_1,d_2) \in \mathbb{R}^n$, denoting $\bar{X}_1$ the sample mean of $X_1$, an equivalent of a AIPW-type estimator takes the form tchetgen2012semiparametric, zhang2013robust, jiang2015doubly, nie2021learning \begin{equation} \begin{aligned} \hat{\mu}_2(d_1, d_2; \hat{\gamma}_1, \hat{\gamma}_2) &= \hat{\gamma}_2(d_{1:2})^\top \Big(Y_2 - H_2 \hat{\beta}_{d_{1:2}}^{(2)}\Big) + \hat{\gamma}_1(d_{1:2})^\top \Big(H_2 \hat{\beta}_{d_{1:2}}^{(2)} - X_1 \hat{\beta}_{d_{1:2}}^{(1)} \Big) + \bar{X}_1 \hat{\beta}_{d_{1:2}}^{(1)}. \end{aligned} \end{equation} We will omit the arguments $(\hat{\gamma}_1, \hat{\gamma}_2)$ in $\hat{\mu}_2$ whenever clear. Its construction directly follows from properties of influence functions tchetgen2012semiparametric. To gain insights, note that a simple estimator for $\mathbb{E}[Y_2(d_1, d_2)]$ is $\bar{X}_1 \hat{\beta}_{d_1, d_2}^{(1)}$. This estimator is consistent and asymptotically normal for low dimensional $\hat{\beta}_{d_1, d_2}^{(1)}$, but not in high-dimensional settings. Instead, the estimator in (ref) uses regression adjustments over each period to control the high-dimensional bias. A choice of the weights from previous literature are inverse probability weights (IPW). These weights for the first and second period are \begin{equation} \begin{aligned} \frac{1\{D_{i,1} = d_1\}}{n P(D_{i,1} = d_1 | X_{i,1})}, \quad \frac{1\{D_{i,1} = d_1\}}{n P(D_{i,1} = d_1 | X_{i,1})} \times \frac{1\{D_{i,2} = d_2\}}{P(D_{i,2} = d_2 | Y_{i,1}, X_{i,1}, X_{i,2}, D_{i,1})}. \end{aligned} \end{equation} However, in high dimensions, AIPW weights require the correct specification of both the propensity score, which in practice may be unknown, and also of the conditional mean function (see Remark (ref)). Also, in small sample, IPW are sensitive to poor overlap (high variance). Motivated by these considerations, we leverage the linear structure to replace IPW with more stable balancing weights that we introduce here. • Main balancing conditions: The main step is in the choice of the balancing weights. The core idea is to decompose \begin{equation} \begin{aligned} \hat{\mu}_2(d_1, d_2) = \bar{X}_1\beta_{d_1, d_2}^{(1)}+ T_1 + T_2 +T_3, \end{aligned} \end{equation} where $\bar{X}_1\beta_{d_1, d_2}^{(1)}$ converges to $\mu(d_1,d_2)$ under standard $\sqrt{n}$-asymptotics, and $$ \small \begin{aligned} T_2 = \hat{\gamma}_2(d_1, d_2)^\top \Big[Y_2 - H_2 \beta_{d_1,d_2}^{(2)}\Big] , \quad T_3= \hat{\gamma}_1(d_1, d_2)^\top \Big[H_2 \beta_{d_1,d_2}^{(2)} - X_1 \beta_{d_1, d_2}^{(1)} \Big] \end{aligned} $$ do not depend on the estimation error for $\beta$. The remaining term $T_1$ is the key component that determines the high-dimensional bias due to the estimation error of $\hat{\beta}$. In particular, writing $T_1 = \Big( \hat{\gamma}_1(d_1, d_2)^\top X_1 - \bar{X}_{1}\Big) (\beta_{d_1, d_2}^{(1)} - \hat{\beta}_{d_1, d_2}^{(1)}) + \Big(\hat{\gamma}_2(d_1, d_2)^\top H_2 - \hat{\gamma}_1(d_1, d_2)^\top H_2\Big) (\beta_{d_1, d_2}^{(2)} - \hat{\beta}_{d_1, d_2}^{(2)})$, we have \begin{equation} \begin{aligned} T_1 \le& \underbrace{\| \hat{\beta}_{d_1,d_2}^{(1)} - \beta_{d_1, d_2}^{(1)}\| _1 \Big| \Big| \bar{X}_1 - \hat{\gamma}_1(d_1, d_2)^\top X_1 \Big| \Big|_{\infty}}_{(i)} + \underbrace{\| \hat{\beta}_{d_1,d_2}^{(2)} - \beta_{d_1, d_2}^{(2)}\| _1 \Big| \Big| \hat{\gamma}_2(d_1, d_2)^\top H_2 - \hat{\gamma}_1(d_1, d_2)^\top H_2 \Big| \Big|_{\infty}}_{(ii)}. \end{aligned} \end{equation} The estimation error depends on the product between the imbalance of covariates characterized by the expressions in $(i), (ii)$ and the estimation error of the coefficients, in the spirit of strong doubly-robustness properties. Therefore a key insight is that to make the estimation error of $\hat{\beta}$ asymptotically negligible over each period, we want to guarantee that \begin{equation} \begin{aligned} \Big| \Big| \bar{X}_1 - \hat{\gamma}_1(d_1, d_2)^\top X_1 \Big| \Big|_{\infty}, \quad \Big| \Big| \hat{\gamma}_2(d_1, d_2)^\top H_2 - \hat{\gamma}_1(d_1, d_2)^\top H_2 \Big| \Big|_{\infty} \end{aligned} \end{equation} are sufficiently small. The first term in Equation (ref) coincides with the balancing term in athey2018approximate. However, here we also require that histories in the second period are balanced, once reweighted by the weights in the previous period. This motivates the dynamic (sequential) balancing weights. To our knowledge, this paper is the first to derive such dynamic balancing conditions. • \textit{Additional conditions on the weights:} The remaining two terms $T_2, T_3$ are mean zero under two conditions: (i) the weights $\hat{\gamma}_2$ are only functions of $H_2$ (and therefore $D_1,D_2$) but not functions of $Y_2$, and weights $\hat{\gamma}_1$ are only functions of $(D_1,X_1)$; (ii) $\hat{\gamma}_{i,2}$ differs from zero only for units with $\{D_{i,2} = d_2, D_{i,1} = d_1\}$ and $\hat{\gamma}_{i,1}$ differs from zero only for units with $D_{i,1} = d_1$. A special case of weights satisfying such conditions are IPW. • \textit{Complete algorithm:} Algorithm 1 presents the algorithmic details in generic $T$ periods. We choose weights sequentially: each period $t$ we minimize the $l_2$ norm of the weights to find stable weights. Such weights are non-zero only for units with treatment history equal to the target path $d_{1:t}$ up-to time $t$; they sum to one and are bounded not to assign too large weight to a few units. The main balancing condition requires that the history $H_{i,t}$ are balanced once reweighting $H_{i,t}$ by the balancing weights estimated in the previous period $t-1$.

We formalize this discussion below (Appendix (ref) contains a formal proof).

thm[Balancing weights] Let assumptions (ref) - (ref) hold. Then Equation (ref) holds with $T_1$ bounded as in Equation (ref). In addition, suppose that $\hat{\gamma}_1$ is measurable with respect to the sigma algebra $\sigma(X_1, D_1)$ and $\hat{\gamma}_2$ is measurable with respect to the sigma algebra $\sigma(X_1, X_2, Y_1, D_1, D_2)$. Suppose in addition that $\hat{\gamma}_{i,1}(d_1, d_2) = 0$ if $D_{i,1} \neq d_1$ and $\hat{\gamma}_{i,2}(d_1, d_2) = 0$ if $(D_{i,1}, D_{i,2}) \neq (d_1, d_2)$. Then $ \mathbb{E}\Big[T_2 \Big| X_1, D_1, Y_1, X_2, D_2\Big] = \mathbb{E}\Big[T_3 \Big| X_1, D_1\Big] = 0. $

rem[Estimating when-to-treat policies] The goal and focus of this paper is to study the effect of a specific counterfactual treatment-assignment sequence, often the object of interest in applications, e.g., robins2000marginal, imai2015robust, acemoglu2019democracy. This differs from studying the effect of an optimal treatment path as in murphy2003optimal or nie2021learning. To gain further insights on the latter problem, define $\pi_t: \mathbb{R}^{p_t} \mapsto \{0,1\}$ a binary policy decision with $\pi(H_{i,t}) \in \{0,1\}$ indicating whether to treat individual $i$ given an observed history $H_{i,t}$ as defined in Equation (ref), and $$ \small \begin{aligned} V(\pi) = \mathbb{E}\left[Y_{i,2}(\pi)\right], \quad Y_{i,2}(\pi):= Y_{i,2}\Big(\pi_1(X_{i,1}), \pi_2(H_{i,2}(\pi_1(X_{i,1}))\Big) \end{aligned} $$ as the expected value of potential outcome $Y_{i,2}$ evaluated under a treatment trajectory $\pi_1(X_{i,1}), \pi_2(H_{i,2}(\pi(X_{i,1})))$. Different from our framework, where the choice of the balancing weight depends on the model for $\mathbb{E}[Y_{i,2}(d_1, d_2) | X_{i,1}]$, in this case the balancing weights must depend on the model for $\mathbb{E}[Y_{i,2}(\pi)| X_{i,1}]$. In the latter case, we must integrate over $\pi_2(H_{i,2}(\pi_1))$ as a function of $H_{i,2}$. Appendix (ref) discusses an extension for this setting.

Generalization to multiple time periods

figure[figure omitted — 116 chars of source]

We now describe in details the procedure with finite $T$ periods. Let $d_{1:T} = (d_1, \cdots, d_T)$,

equation[equation omitted — 181 chars of source]

This estimand denotes the difference in potential outcomes for two treatment histories $d_{1:T}, d_{1:T}'$. We denote

equation[equation omitted — 196 chars of source]

the vector containing information from time one to time $t$, after excluding the treatment assigned in the present period $D_t$. Interaction components may also be considered, omitted here for brevity. We let the potential history be $ H_{i,t}(d_{1:(t-1)}) = \Big[d_{1:(t-1)}, X_{i,1:t}(d_{1:(t-1)}), Y_{i,1:(t-1)}(d_{1:(t-1)}) \Big]. $

assFor any $d_{1:T} \in \{0,1\}^T$, and $t \le T$, \begin{itemize} • (No-anticipation) The potential history $H_{i,t}(d_{1:T})$ is constant in $d_{t:T}$; • (Sequential ignorability) $\Big(Y_{i,T}(d_{1:T}), H_{i,t+1}(d_{1:(t+1)}), \cdots, H_{i,T-1}(d_{1:(T-1)})\Big) \perp D_{i,t} | H_t$; • (Potential projections) For some $\beta_{d_{1:T}}^{(t)} \in \mathbb{R}^{p_t}$, $$ \small \begin{aligned} \mathbb{E}\Big[Y_{i,T}(d_{1:T}) | D_{i, 1:(t-1)} = d_{1:(t-1)}, X_{i,1:t}, Y_{i,1:(t-1)}\Big] = H_{i,t}(d_{1:(t-1)}) \beta_{d_{1:T}}^{(t)}. \end{aligned} $$ \end{itemize}

Assumption (ref) generalizes Assumptions (ref)-(ref) from the two-period setting. Identification follows similarly to Lemma (ref). For given weights $\hat{\gamma}_{1:T}$, and coefficients $\hat{\beta}^{(1:T)}$, we estimate

equation[equation omitted — 383 chars of source]

Coefficients are estimated recursively as in the two periods setting (see Algorithm 2).

lemLet $\hat{\gamma}_{i,T}(d_{1:T}) = 0$ if $D_{i,1:T} \neq d_{1:T}$. For $\varepsilon_{i,t}(d_{1:T}) = Y_{i,T}(d_{1:T}) - H_{i,t}(d_{1:(t-1)}) \beta_{d_{1:T}}^{(t)}$, \begin{equation} \begin{aligned} \hat{\mu}_T(d_{1:T}) - \bar{X}_1 \beta_{d_{1:T}}^{(1)} = &\underbrace{\sum_{t = 1}^T \Big(\hat{\gamma}_t(d_{1:T}) H_t - \hat{\gamma}_{t-1} (d_{1:T})H_t\Big) (\beta_{d_{1:T}}^{(t)} - \hat{\beta}_{d_{1:T}}^{(t)})}_{(I_1)} + \underbrace{\hat{\gamma}_T^\top(d_{1:T}) \varepsilon_T}_{(I_2)} \\ & + \underbrace{\sum_{t=2}^T \hat{\gamma}_{t-1}(d_{1:T})\Big(H_t \beta_{d_{1:T}}^{(t)} - H_{t-1} \beta_{d_{1:T}}^{(t-1)}\Big)}_{(I_3)} \end{aligned} \end{equation}

The proof is in Appendix (ref). Lemma (ref) decomposes the estimation error into three components. First, ($I_1$), depends on the estimation error of the coefficient and on balancing properties of the weights. ($I_1$) suggests imposing balancing conditions on \\ $ \Big| \Big| \hat{\gamma}_t(d_{1:T}) H_t - \hat{\gamma}_{t-1}(d_{1:T}) H_t \Big| \Big|_{\infty} $ each period. The components characterizing the estimation error are $ (I_2)=\hat{\gamma}_T(d_{1:T})^\top \varepsilon_T$, and ($I_3$), which are mean zero as described in Appendix Lemma (ref). Note that $\bar{X}_1 \beta_{d_{1:T}}^{(1)}$ converges to $\mu_T(d_{1:T})$ under standard $\sqrt{n}$-asymptotics.

Theoretical properties and inference

Next, we study the theoretical properties of the estimator in finite $T$ periods. We consider a high dimensional regime where the dimension covariates in each period $p_1, \cdots, p_T$ can grow to infinity, as long as $\log (\max_t p_t n)/n^{1/4} \to 0$. We impose the following conditions.

ass[Overlap and tails' conditions] Assume that (i) ${P(D_{i,t} = d_t | H_{i,t})} \in (\delta,1 - \delta), \delta \in (0,1)$ for each $t \in \{1, \cdots, T\}$; and (ii) $H_{i,t}^{(j)}, \forall j$ is Sub-Gaussian given $H_{i,t-1}$ and $X_{i,1}^{(j)}, j \in \{1, \cdots, p_1\}$ is Sub-Gaussian.

Condition (i) is the overlap condition, standard in the causal inference literature. The overlap condition is sufficient (but not necessary, see the discussion of athey2018approximate in cross-sectional settings) to show existence of a feasible solution of Algorithm 1; see Remark (ref). Condition (ii) is a tail restriction. Assumption (ref) can be relaxed by assuming that the product of the inverse probability weights times the covariates is sub-exponential at the expense of more tedious derivations.

thm[Existence of feasible weights] Let Assumptions (ref), (ref) hold. Consider $\delta_t(n,p_t) \ge c_{0,t} {n^{-1/2}}{\log^{3/2}(p_tn)}$ for a finite constant $c_{0,t} < \infty$, and $C_{n,t} \ge \frac{\bar{c}}{n \delta^t}$ for some sufficiently large constant $\bar{c} \in (0,\infty)$. Then with probability $\eta_n \rightarrow 1$, for each $t \in \{1, \cdots, T\}, T < \infty$, for some $N >0$, $n > N$, there exists a $\hat{\gamma}_t^*(\hat{\gamma}_{t-1})$ satisfying the constraint in Algorithm 1 as a function of the solution $\hat{\gamma}_{t-1}$ of Algorithm 1, where \begin{equation} \begin{aligned} \hat{\gamma}_{i,0}^* = 1/n, \quad \hat{\gamma}_{i,t}^*(\hat{\gamma}_{t-1}) = \hat{\gamma}_{i,t-1} w_{i,t}^* \Big/ \sum_{i=1}^n \hat{\gamma}_{i,t-1} w_{i,t}^*, \quad w_{i,t}^*:= \frac{1\{D_{i,t} = d_t\}}{P(D_{i,t} = d_t | H_{i,t})} \end{aligned} \end{equation}

The proof is in Appendix (ref). In particular, the reader may refer to Appendix Lemma (ref) for additional details on the construction of the feasible weights. The algorithm thus finds weights that minimize the small sample variance, with the IPW weights $w_{i,t}^*$ (reweighted by the solution in the previous iteration $\hat{\gamma}_{i,t}$) being one possible solution. Minimizing the $l_2$ norms of the weights is a natural objective when the goal is to minimize the variance of the ATE estimator: following Theorem (ref) below, under homoskedasticity of the residuals from each projection, the variance is proportional to a weighted sum of $||\hat{\gamma}_t||^2$. However, homoskedasticity is not necessary for our results.

corLet the conditions in Theorem (ref) hold. Then for some $N > 0$, $n > N$, with probability $\eta_n \rightarrow 1$, there exist a sequence of feasible solutions $\Big\{\hat{\gamma}_t\Big\}_{t=1}^T$ such that $ n ||\hat{\gamma}_{t}||^2 \le n||\hat{\gamma}_{t}^*(\hat{\gamma}_{t-1})||^2$ for each $t$, where $\hat{\gamma}_{i,t}^*(\hat{\gamma}_{t-1})$ is as defined in Equation (ref), and $n ||\hat{\gamma}_t||^2 \le n c_t ||\hat{\gamma}_{t-1}||^2$ for a finite constant $c_t < \infty$.

Corollary (ref) provides a desiderable stability property. It shows that the $l_2$ norm of the weights is upper bounded by the (stabilized) IPW weights reweighted by the solution in the previous step. From basic concentration inequalities, for $n$ sufficiently large, $\mathbb{E}[\hat{\gamma}_{i,t}^{*2}(\hat{\gamma}_{t-1})] \approx \mathbb{E}\Big[\frac{ \hat{\gamma}_{i,t-1}^2 }{P(D_{i,t} = d_t | H_{t})}\Big]$. In addition, the last part of the corollary shows that the weights' norm is controlled over each period by the norm in the previous period up to a finite multiplicative constant.

assLet the following hold: for every $t \in \{1, \cdots, T\}$, $d_{1:T} \in \{0,1\}^T$, \begin{itemize} • $\max_t \| \hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)} \| _1 \delta_t(n,p_t) = o_p(1/\sqrt{n})$, $\delta_t(n,p_t) \ge c_{0, t} {n^{-1/2}}{\log^{3/2}(p_tn)}$ for a finite constant $c_{0, t}$; In addition let either of the two following conditions hold: (a) $\max_t \| \hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)} \|_1 = O_p(n^{-1/4})$; or (b) $\max_t \| \hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)} \|_1 = o_p(1/\log(n))$ and $||H_t||_{\infty} \le \bar{H}$ almost surely for a finite constant $\bar{H} < \infty$ for all $t \ge 1$; • Let $\nu_{i,t} = (H_{i,t +1} \beta_{d_{1:T}}^{(t+1)} - H_{i,t} \beta_{d_{1:T}}^{(t)}), \varepsilon_{i,T} = Y_{i,T} - H_{i,T} \beta_{d_{1:T}}^{(T)}$. For a finite constant $C$, $\mathbb{E}[\varepsilon_{i, T}^4 | H_{i, T}, D_{i, T}] < C$, $\mathbb{E}[\nu_{i,t}^4 | H_{i,t-1}, D_{i,t-1}] < C$, almost surely. In addition, $Y_{i,T}$ is a sub-gaussian random variable. • $\mbox{Var}(\varepsilon_{i,T} | H_{i,T}, D_{i, T}), \mbox{Var}(H_{i,t} \beta_{d_{1:T}}^{(t)} - H_{i,t-1} \beta_{d_{1:T}}^{(t-1)} | H_{i,t-1}, D_{i,t-1}) > u_{min}$, almost surely, for some constant $u_{min} > 0$. \end{itemize}

Assumption (ref) imposes the consistency in estimating the outcome models. Condition (i) is attained for many high-dimensional estimators, such as the lasso method, under sparsity and restricted eigenvalues restrictions; see, e.g., buhlmann2011statistics. An example and derivation for condition (i) for Lasso under sparsity is included in Example (ref) (Appendix (ref)). As we discuss in Appendix Remark (ref), the restricted eigenvalue condition is imposed for all $H_t, t \ge 1$; references that study this condition from different angles include deshpande2023online, Section C.4.

thm[Parametric convergence rate] Let the conditions in Theorem (ref) and Assumption (ref) hold. Then, whenever $\log(n(\sum_t p_t))/n^{1/4} \to 0$ with $n,p_1, \cdots,p_T\to \infty$, it follows that $ \hat{\mu}_T(d_{1:T}) - \mu_T(d_{1:T}') = \mathcal{O}_P\Big(n^{-1/2}\Big). $

Theorem (ref) guarantees a parametric convergence rate with high-dimensional covariates.

thm[Inference] Let the conditions in Theorem (ref) and Assumption (ref) hold. Then, whenever $\log(n \sum_t p_t)/n^{1/4} \to 0$, as $n, p_1, \cdots, p_T \rightarrow \infty$, \begin{equation} \begin{aligned} { \frac{\sqrt{n} \Big(\hat{\mu}(d_{1:T}) - \mu(d_{1:T})\Big)}{\hat{V}_T(d_{1:T})^{1/2}}} \rightarrow_d \mathcal{N}(0,1) \end{aligned} \end{equation} where $$ \small \begin{aligned} & \hat{V}_T(d_{1:T}) = \\ & \sum_{i = 1}^n \left\{ n\hat{\gamma}_{i,T}^2(d_{1:T}) (Y_{i, T} - H_{i,T} \hat{\beta}_{d_{1:T}}^{(T)})^2 + \sum_{t=1}^{T-1} n \hat{\gamma}_{i,t}^2(d_{1:t}) (H_{i,t+1} \hat{\beta}_{d_{1:T}}^{t+1} - H_{i,t} \hat{\beta}_{d_{1:T}}^{t})^2 + \frac{1}{n} (\bar{X}_1 \hat{\beta}_{d_{1:T}}^{(1)} - X_{i,1} \hat{\beta}_{d_{1:T}}^{(1)})^2\right\} \end{aligned}. $$

Inference on ATE follows as a direct corollary for two histories $d_{1:T}, d_{1:T}'$ with $d_1 \neq d_1'$ (see Theorem (ref)), as described in Appendix (ref). Also, for inference conditional on baseline covariates in the first period $X_1$, the relevant variance is $\hat{V}_T(d_{1:T}) - \frac{1}{n} \sum_i (\bar{X}_1 \hat{\beta}_{d_{1:T}}^{(1)} - X_{i,1} \hat{\beta}_{d_{1:T}}^{(1)})^2$ (since we condition on $X_1$) and for the corresponding ATE is the sum of these two variances.

rem[Strict overlap assumption] Strict overlap (Assumption (ref) $(i)$) is not necessary to achieve parametric convergence rates whenever a feasible solution to Algorithm 1 exists (see for instance athey2018approximate, Lemma 2 in cross-sectional settings). \qed

rem[Rate conditions and comparison with AIPW] As noted in hirshberg2021augmented in static settings, the advantages of balancing weights compared to AIPW is to be able to estimate causal effects of interest under essentially the same conditions for AIPW on the conditional mean function, but weaker conditions on the balancing weights. Specifically, in high-dimensional settings, AIPW requires conditions of the form $||\hat{e} - e|| = o_P(n^{-1/4}), ||\hat{\beta} - \beta|| = o_P(n^{-1/4})$, where $e$ denotes the propensity score. Here, we only require $||\hat{\beta} - \beta||_1 = O_p(n^{-1/4})$ with sub-gaussian histories $H_t$ and $||\hat{\beta} - \beta||_1 = o_p(1/\log(n))$ with uniformly bounded covariates and no condition on the propensity score. We could think of balancing weights as inheriting a “product-of-rate" condition, where here the estimation error takes the form: $ ||\hat{\beta} - \beta||_1 \delta(n,p). $

rem[Propagation of error over multiple periods] A natural question is how estimation error varies as the number of periods increase. To shed light on this question, note that the estimation bias for $\hat{\mu}(d_{1:T})$ is bounded above by $$ \sum_{t=1}^T ||\hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)} ||_1 \Big| \Big| \sum_{i=1}^n \hat{\gamma}_{i,t-1} H_{i,t} - \sum_{i=1}^n \hat{\gamma}_{i,t} H_{i,t} \Big| \Big|_{\infty}. $$ Therefore, multiple time periods may affect the error through the estimation error of the coefficient, and through the weights $\hat{\gamma}_t$ in the balancing component. The properties of estimated coefficients for Lasso are presented in Appendix (ref). For the balancing component instead, existence of a feasible solution requires that the balancing component grows at rate $\sqrt{\log(t)}$, so that effectively the constants are of order $K_{1,t} = \log^{1/2}(t)$ (see Appendix Lemma (ref)). Intuitively, as we move along the path over multiple time periods it becomes harder to guarantee approximate balance, reflecting into weaker constraints on the balancing set.

Guide to practice: numerical studies and application

Implementation guide

The complete Algorithm 1 is implemented off-the-shelf in the R-package {\tt DynBalancing}.

It requires researchers to specify four main parameters: the length $h$ of the treatment history considered (i.e., carry-over effects), two treatment histories of length $h$, $d_{(T-h):T}, d_{(T-h):T}'$ to compare, the model used to estimate the coefficients ({\tt linear} or {\tt fully interacted}) as described in Algorithm 2, and whether to consider a {\tt pooled} regression.

Choosing the length of the treatment history with long panel With short panels, selecting the length of the treatment history $h = T$ is natural. With long panels, this may reduce the effective sample size or be infeasible (as the effective sample becomes “thinner"). This is because, as for IPW, the weights at time $t$ can be non zero only for those units observed over a given treatment path up to time $t$. Therefore, we recommend selecting a treatment history $h$ shorter than the number of periods $T$ (i.e., $h < T$), and estimate causal effects of the form

equation[equation omitted — 245 chars of source]

for given treatment histories $d_{(T-h):T}, d_{(T-h):T}'$. Equation (ref) estimates the effect of exposing an individual to two different histories over the last $h$ periods and average over previous assignments. Our analysis and estimation follow similarly to Algorithm 1, with the difference that we construct balancing weights starting from period $T - h$ and proceed sequentially until time $T$ (observable characteristics before time $T - h$ can be used as additional controls). As in imai2018matching, the focus on Equation (ref) makes our procedure robust to long panels.

As a rule of thumb, as we illustrate in our application, we recommend report results for different choices of $h$ (say $h \in \{1, \cdots, 10\}$ in a long panel); our package reports and plots the estimated effects along-side standard errors which can help disentangle the trade-offs between identification of long-run effects against precision. In addition, it is useful to report $1/(n ||\hat{\gamma}_t||^2)$ as a measure of effective sample size at time $t$ for different values of $t$, which can help guide the choice of $h$ (larger choices of $1/(n ||\hat{\gamma}_t||^2)$ indicates more accurate treatment effects estimates). \\ Choice of the model specification ({\tt linear} or {\tt fully interacted}) The estimation error $||\hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)}||_1$ depends on modeling assumptions. For the {\tt fully interacted} model, $||\hat{\beta}_{d_{1:T}}^{(t)} - \beta_{d_{1:T}}^{(t)}||_1$ scales exponentially with $T$ as it considers all possible interactions with the treatment assignments $d_1, \cdots, d_T$. The {\tt linear} model avoids that the effective sample size shrinks exponentially in $T$ but imposes homogeneity restrictions of treatment effects as for example in acemoglu2019democracy, by modeling treatment effects as additive and linear. See Algorithm 2 for more details. In addition, when {\tt pooled} is true, we consider a regression $$

alignedY_{i,t}(d_{1:t}) = \beta_0 + \beta_1 d_t + \beta_2 Y_{i,t-1}(d_{1:(t-1)}) + X_{i,t}(d_{1:(t-1)}) \gamma + \tau_t + \varepsilon_{i,t},

$$ where $\tau_t$ denotes fixed effects, pooling together effects estimated in different periods. We then cluster standard errors at the individual level to allow for correlation over time.

figure[figure omitted — 116 chars of source]

rem[Tuning parameters] Similarly to one-dimensional setting athey2018approximate, Algorithm 1 requires choosing tuning parameters. A complete description is in Algorithm (ref) and uses a data-adaptive procedure (i.e., researchers do not need to specify the tuning parameters). In a nutshell, we choose $\delta_t(n,p) = \log^{3/2}(p_t n)/n^{1/2}$ (here $p_t$ is the dimension of covariates at time $t$) as prescribed by the theoretical analysis in Section (ref) and for simplicity we set $C_{n,t} = \log(n) n^{-2/3}$. To guarantee balance with many covariates, we iterate over a grid of values for $K_1$, and select the smallest constant within this grid such that the optimization program admits a feasible solution (in Algorithm 1 we also refine the algorithm to weight more the constraints for which the coefficients are non-zero). This approach minimizes the estimator's bias and, within the set of weighting estimators with the smallest bias, selects the one with the smallest variance. Below and in Appendix (ref) we illustrate the benefits of this procedure through numerical studies. \qed

rem[Computational complexity] Algorithm 1 is a sequence of $T$ quadratic programs with linear constraints. Its complexity scales only polynomially with $n, p$. Appendix Figure (ref) shows that the computational time is between a few seconds and a few minutes for $T \in \{1, \cdots, 10\}$ on a personal laptop (including choosing the tuning parameters). \qed

Numerical studies

Next, we collect results from numerical experiments. We estimate $ \mathbb{E}\Big[Y_{i,T}(1, \cdots, 1) - Y_{i,T}(0, \cdots,0)\Big], T \in \{2,3\}. $ We let the baseline covariates $X_{i,1}$ be drawn from as i.i.d. $\mathcal{N}(0, \Sigma)$ with $\Sigma^{(i,j)} = 0.5^{|i-j|}$. Covariates in the subsequent period are generated according to an auto-regressive model $ \{X_{i,t}\}_j =0.5 \{X_{i,t-1}\}_j + \mathcal{N}(0, 1), j=1,\cdots,p_t. $ Treatments are drawn from a logistic model that depends on all previous treatments and past covariates: $ D_{i,t} \sim \mbox{Bern}\Big((1 + e^{\iota_{i,t}})^{-1}\Big) $ with

equation[equation omitted — 231 chars of source]

and $\xi_{i,t} \sim \mathcal{N}(0,1)$, for $t \in \{1, 2,3\}$. Here, $\eta, \delta$ controls the association between covariates and treatment assignments. We consider values of $\eta \in \{0.1, 0.3, 0.5\}$, $\delta_1 = 0.5, \delta_2 = 0.25$. We let $\phi \propto 1/j$, with $\|\phi \|_2^2 = 1$, similarly to balancing conditions presented in athey2018approximate. The larger $\eta$ corresponds to weaker overlap (see Table (ref) in the Appendix).

We generate the outcome as $ Y_{i,t}(d_{1:t}) = \sum_{s = 1}^t \Big(X_{i,s} \beta + \lambda_{s, t} Y_{i,s-1} + \tau d_s\Big) + \varepsilon_{i,t}(d_{1:t}), \quad t=1,2,3, $ where elements of $\varepsilon_{i,t}(d_{1:t})$ are i.i.d. $ \mathcal{N}(0,1)$ and $\lambda_{1,2} = 1, \lambda_{1,3}, \lambda_{2,3} = 0.5$. We consider three different settings: Sparse with $\beta^{(j)} \propto 1\{j \le 10\}$, Moderate with moderately sparse $\beta^{(j)} \propto 1/j^2$ and the Harmonic setting with $\beta^{(j)} \propto 1/j$. We set $\| \beta \|_2 =1, \tau = 1$.

We consider the following competing methodologies:

itemize• Augmented IPW, with known propensity score and with estimated propensity score. The method replaces the balancing weights in Equation (ref) with the (estimated or known) propensity score. Estimation of the propensity score is performed using a logistic regression (denoted as aIPWl) and a penalized logistic regression (denoted as aIPWh) nie2021learning, bodoryevaluating. For both AIPW and IPW we consider stabilized inverse probability weights. • CAEW (MSM): Although our balancing weights in Algorithm 1 are novel (with or without regression adjustment), we can compare to other balancing procedures. In particular, we consider Marginal Structural Model (MSM) with balancing weights computed using the method in yiu2018covariate, yiu2020joint that, different from ours, require information about the propensity score. We follow Section 3 in yiu2020joint for its implementation. (We do not also compare to imai2015robust for MSM since it is not feasible in high-dimensions.) • “Dynamic" Double Lasso: it estimates the effect of each treatment assignment separately, after conditioning on the present covariate and past history for each period using the double lasso discussed in one period from belloni2014inference. • Naive Lasso: it runs a regression controlling for covariates and treatment assignments. • \textit{Sequential Estimation}: it estimates the conditional mean in each time period sequentially using the lasso method, and it predicts end-line potential outcomes as a function of the estimated potential outcomes in previous periods. • \textit{DiD switchback}: it is a DiD estimator similar to de2024difference. • \textit{Simple LP} (Local Projection): it projects $Y_T$ onto baseline covariates $X_{i,1}$ and treatment $D_{i,1}$ and take the coefficient multiplying $D_{i,1}$ as the estimated effect, while penalizing the coefficients for $X_{i,1}$ via Lasso.

For Dynamic Covariate Balancing, DCB, the choice of tuning parameters is data adaptive, and it uses a grid-search method discussed in Appendix (ref) and Remark (ref). We estimate coefficients as in Algorithm 1 for DCB and (a)IPW, with a linear model in treatment assignments. Estimation of the penalty for the lasso methods is performed via cross-validation.

We consider $\mathrm{dim}(\beta) = \mathrm{dim}(\phi) = 100$ and set the sample size to be $n = 400$. We set $p_1=101, p_2=203, p_3=305$ as number of covariates in each period.

In Table (ref) we collect results for the average mean squared error for estimating the average treatment effect in two and three periods. Throughout all simulations, the proposed method significantly outperforms any other competitor, with one single exception for $T = 2$, good overlap and harmonic design. It also outperforms using known propensity score, consistently with our findings in Theorem (ref), where we show that the propensity score is a feasible solution of DCB weights (and in the absence of knowledge of the propensity score).

Finally, in Appendix (ref) we consider more extensive simulation studies with a longer time horizon, a misspecified (non-linear) model, low and high dimensional settings among additional simulation designs.

\definecolor{glaucous}{rgb}{0.38,0.51,0.71}

table*[table* omitted — 3,704 chars of source]

Empirical illustration

In this section, we present an empirical application for studying the effect of democracy on GDP growth using data from acemoglu2019democracy. acemoglu2019democracy studied dynamic treatment effects of democracy under GDP growth under sequential ignorability acemoglu2019democracy. Figure (ref) illustrates the dynamics of treatments. Many units switch treatment over time, violating standard event studies designs.

The data (available at \url{https://www.journals.uchicago.edu/doi/suppl/10.1086/700936}) consist of a collection of countries observed between $1960$ and $2010$. We consider observations starting from $1989$. After removing missing values, we run regressions with 141 countries. The outcome is the log-GDP in the country $i$ in period $t$ as in acemoglu2019democracy. We use the same treatment specification as in acemoglu2019democracy, which is binary. We study the effect of exposing countries at time $t$ to democracy for in $s$ years before (and including) $t$ versus not exposing them to democracy for the previous $s$ years. Namely, the estimand is the $s$-long run effect of democracy, after averaging over past assignments. We let $s \in \{1, \cdots, 20\}$ to study the impact from one to twenty years of democracy.

For each country, we condition on lag outcomes in the past four years as in the preferred specification of acemoglu2019democracy, and past four treatments. We consider a pooled regression and two alternative specifications. The first is parsimonious and includes dummies for different regions (continents) and different intercepts for different periods. The second one controls for the past four outcomes, past four treatments, for the geographical region, and colonial history as in acemoglu2019democracy. Coefficients are estimated as in Algorithm 2 with {\tt model} $=$ {\tt linear}.

\definecolor{aero}{rgb}{0.49,0.73,0.91} \definecolor{airsuperiorityblue}{rgb}{0.45,0.63,0.76} \definecolor{babyblueeyes}{rgb}{0.63,0.79,0.95} \definecolor{beaublue}{rgb}{0.74,0.83,0.9} \definecolor{glaucous}{rgb}{0.38,0.51,0.71}

figure[figure omitted — 197 chars of source]

figure[figure omitted — 571 chars of source]

table[table omitted — 889 chars of source]

Figure (ref) collects our results. Democracy has a statistically insignificant effect over the first few years and a statistically significant positive impact on long-run GDP growth after three years. Point estimates are in sign and magnitude consistent with what found by acemoglu2019democracy, and results are robust across the two specifications for DCB.

We compare our method to (i) the linear estimator reported by acemoglu2019democracy (Table 2, Column 3), where dynamic effects are estimated by propagating the effect over past outcomes at each period (we consider two specifications, with and without unit fixed effects -- both report similar results); (ii) the simple local projection, that projects the outcome on the treatment and the past outcome $s$ periods before, with and without country fixed effects, time fixed effects and controlling for lagged outcome at time $s$.

The simple local projection approach reports small point estimates compared to other methods. This result is consistent with our theoretical discussion: local projections average over the distribution of future assignments. Therefore, the causal effects estimated by the local projection differ from the target long-run effect, which instead fixes future treatment assignments. The effect estimated as in acemoglu2019democracy is larger than the local projection when including country fixed effects, but significantly smaller than the effect estimated through DCB. Therefore, the specification in acemoglu2019democracy may capture some but not all the long-run effects. After controlling for imbalance with DCB, average treatment effects are twice as large. The results from acemoglu2019democracy with and without unit fixed effects report almost identical results.

To investigate differences with (A)IPW methods, the right panel in Figure (ref) presents comparisons in terms of the imbalance over the lagged outcome at time $t - 1$ when using balancing or inverse probability weights. As acemoglu2019democracy note, the lags outcome may capture most variation in treatment. Therefore, an imbalance in lagged GDP may suggest the presence of bias. We report the relative improvement in absolute imbalance (average across the potential outcomes under treatment and control) and observe substantial gain over using inverse probability weights. Such gains illustrate the advantage of balancing in small sample.

Figure (ref) complements Figure (ref) showing instability of inverse probability weights, and Figure (ref) in the Appendix show that DCB weights present less dispersion than IPW weights. As a helpful diagnostic, in Table (ref) we report $1/(||\hat{\gamma}_t||^2)$ for the DCB method as well as for IPW weights, where $\hat{\gamma}_t$ is replaced by the AIPW weight. This measure is indicative of the level of precision and effective sample size (a larger number indicates better precision). We find substantial improvements of DCB over IPW especially for weights estimated for longer time horizons.

figure[figure omitted — 736 chars of source]