The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
77,789 characters
Dynamic covariate balancing: estimating treatment effects over time with potential local projections
\maketitle
\begin{abstract}
This 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.
\end{abstract}
\begin{keywords}
Causal Inference, High Dimensions, Treatment Effects, Panel Data.
\end{keywords}
\section{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 \citep{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 \textit{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 \citep[][]{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 \citep[][]{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 \citep{ding2019bracketing}. Researchers remain agnostic on the propensity score.
We leverage a model on the \textit{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 \citep{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 \citep{robins2000marginal}, we identify the parameters of interest by \textit{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 \citep{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 \citep{imai2014covariate, li2018balancing, hainmueller2012entropy}.
We generalize covariate balancing in the absence of dynamics of \cite{zubizarreta2015stable, athey2018approximate, ben2018augmented, hirshberg2017augmented} to a dynamic setting. We show that balancing with potential local projections corresponds to constructing weights \textit{sequentially} in time by first balancing treated and control units' covariates in the first period and then balancing histories in the next periods \textit{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 \textit{both} the outcome model and propensity score model \citep[what known as rate-doubly robustness, see e.g.,][]{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{fig:overlap}). 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 \cite{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 \cite{acemoglu2019democracy}, we show that standard local projections and \cite{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.
\section{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,
\cite{imbens2019panel} propose balancing assuming no treatment dynamics, whereas here, treatment dynamics require different (and novel) balancing conditions. In the context of dynamics, different from \cite{imai2015robust}, who estimate a \textit{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 \cite{zhou2018residual}, who extend entropy balancing of \cite {hainmueller2012entropy} to dynamic settings, and \cite{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 \cite{kallus2018optimal} design conservative balancing conditions for the worst-case bias. Different from \cite{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 \textit{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 \cite{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 \citep[e.g.][]{bruns2023augmented}, and our analysis differs from cross-sectional studies with low dimensions in \cite{wang2020minimal}.
More broadly, this paper connects to the literature on
DiD, local projections, and dynamic treatments.
Different from the literature on DiD \citep{rambachan2023more, de2022difference, callaway2019difference, abraham2018estimating, athey2022design, caetano2022difference} or subsequent work on potential projections with DiD \citep{dube2023local}, here we allow for dynamic treatment regimes. This literature imposes the parallel trends assumption, violated with dynamic treatments \citep{marx2022parallel, ghanem2022selection}.
Different from the time-series literature \citep{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 \cite{robins2000marginal}, \cite{hernan2001marginal}, \cite{boruvka2018assessing}, \cite{blackwell2013framework}, \cite{bang2005doubly} \citep[for a review, ][]{vansteelandt2014structural}. \cite{bojinov2020panel} study IPW estimators from a design-based perspective.
Doubly robust estimators for dynamic treatments have been studied by \cite{nie2021learning, zhang2013robust, jiang2015doubly, tchetgen2012semiparametric, babino2019multiple}. Here, we focus on studying the effect of a given treatment path, as in \cite{robins2000marginal} or \cite{blackwell2013framework}, different from and complementary to studying optimal policies (e.g., \cite{murphy2003optimal}, \cite{nie2021learning}).
Specifically, studies with high-dimensional panels require correct specification of the propensity score
\citep{lewis2020double, zhu2017high, shi2018high, bodoryevaluating, belloni2016inference, chernozhukov2017orthogonal}, or impose homogeneous treatment effects \citep{high_dim_IRF, kock2015inference}.
(\cite{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 \cite{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, \cite{chernozhukov2022automatic} generalize the use of riesz representers with arbitrary non-linear outcome models, and \cite{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 \cite{chernozhukov2022automatic} and \cite{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.
\section{Dynamics and potential local projections} \label{sec:2}
\subsection{Setup}
We start with the analysis of two time periods, deferring multiple periods to Section \ref{sec:multiple}. 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.
\begin{figure}[!ht]
\centering
\begin{tikzpicture}
\coordinate (1) at (-4,3);
\coordinate (2) at (-2,3);
\coordinate (3) at (-2,5);
\coordinate (4) at (-4,5);
\coordinate (5) at ($(1)!.5!(2)$);
\coordinate (6) at ($(2)!.5!(3)$);
\coordinate (7) at ($(3)!.5!(4)$);
\coordinate (8) at ($(1)!.5!(4)$);
\coordinate (9) at ($(1)!.5!(3)$);
\coordinate (10) at (-4,0);
\coordinate (11) at (-2,0);
\coordinate (12) at (-2,2);
\coordinate (13) at (-4,2);
\coordinate (14) at ($(10)!.5!(11)$);
\coordinate (15) at ($(11)!.5!(12)$);
\coordinate (16) at ($(12)!.5!(13)$);
\coordinate (17) at ($(10)!.5!(13)$);
\coordinate (18) at ($(10)!.5!(14)$);
\coordinate (21) at (-10,4);
\coordinate (22) at (3.5,4);
\coordinate (23) at (3.5,5);
\coordinate (24) at (-10, 5);
\coordinate (31) at (-10,3.5);
\coordinate (32) at (-7,3.5);
\coordinate (33) at (-7,2.5);
\coordinate (34) at (-10, 2.5);
\coordinate (41) at (-6.5,3.5);
\coordinate (42) at (-3.5,3.5);
\coordinate (43) at (-3.5,2.5);
\coordinate (44) at (-6.5, 2.5);
\coordinate (51) at (-3,3.5);
\coordinate (52) at (0,3.5);
\coordinate (53) at (0,2.5);
\coordinate (54) at (-3, 2.5);
\coordinate (61) at (0.5,3.5);
\coordinate (62) at (3.5,3.5);
\coordinate (63) at (3.5,2.5);
\coordinate (64) at (0.5, 2.5);
\coordinate (71) at (-10,2);
\coordinate (72) at (-8.6,2);
\coordinate (73) at (-8.6,1);
\coordinate (74) at (-10, 1);
\coordinate (81) at (-8.4,2);
\coordinate (82) at (-7,2);
\coordinate (83) at (-7,1);
\coordinate (84) at (-8.4, 1);
\coordinate (91) at (-6.5,2);
\coordinate (92) at (-5.1,2);
\coordinate (93) at (-5.1,1);
\coordinate (94) at (-6.5, 1);
\coordinate (101) at (-4.9,2);
\coordinate (102) at (-3.5,2);
\coordinate (103) at (-3.5,1);
\coordinate (104) at (-4.9, 1);
\coordinate (111) at (-3,2);
\coordinate (112) at (-1.6,2);
\coordinate (113) at (-1.6,1);
\coordinate (114) at (-3, 1);
\coordinate (121) at (-1.4,2);
\coordinate (122) at (0,2);
\coordinate (123) at (0,1);
\coordinate (124) at (-1.4, 1);
\coordinate (131) at (0.5,2);
\coordinate (132) at (1.9,2);
\coordinate (133) at (1.9,1);
\coordinate (134) at (0.5, 1);
\coordinate (141) at (2.1,2);
\coordinate (142) at (3.5,2);
\coordinate (143) at (3.5,1);
\coordinate (144) at (2.1, 1);
\draw[->] (-8,4.3) -- (-1,4.3);
\node[circle] (g) at (-8,4.3) {$|$};
\node[circle] (g) at (-5,4.3) {$|$};
\node[circle] (g) at (-2,4.3) {$|$};
\node[circle] (g) at (-8,4) {$t = 0$};
\node[circle] (g) at (-5,4) {$t = 1$};
\node[circle] (g) at (-2,4) {$t = 2$};
\node[circle] (g) at (-8,5) {$X_1,$};
\node[circle] (g) at (-6.8,5) {$D_1,$};
\node[circle] (g) at (-5,5) {$(Y_1, X_2),$};
\node[circle] (g) at (-3.3,5) {$D_2,$};
\node[circle] (g) at (-2,5) {$Y_2$};
\end{tikzpicture}
\caption{Sampling process in two periods. First, baseline covariates $X_{i,1}$ realize at $t = 0$. Then, treatment $D_{i,1}$ is assigned and the outcomes and covariates $(Y_{i,1},X_{i,2})$ realize at $t = 1$. Finally, the treatment $D_{i,2}$ is assigned and, afterwards, the endline outcome $Y_{i,2}$ realizes. } \label{fig:time}
\end{figure}
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)
$$
\small
\begin{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],
\end{aligned}
$$
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{fig:seqign} 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{sec:app}).
\begin{figure}[!ht]
\centering
\scalebox{0.7}{
\begin{tikzpicture}[scale = 1.16]
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (h) at (-3,-2) {$D_1$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (e) at (-0.5,-2) {$Y_1, X_2$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (f) at (-0.5,-4) {$D_2$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=13.5mm,
align=center, circle] (d) at (2,-2) {$Y_2$};
\draw[->, -triangle 90] (e) edge (f) (h) edge (f);
\draw[->, -triangle 90] (h) edge (e) (e) edge (d) (f) edge (d) ;
\draw[->, -triangle 90] (h) edge[bend right=-30] node [left] {} (d);
\end{tikzpicture}}
\scalebox{0.7}{
\begin{tikzpicture}[scale = 1.16]
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (h) at (-3,-2) {$D_1$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (e) at (-0.5,-2) {$Y_1, X_2$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=14mm,
align=center, circle] (f) at (-0.5,-4) {$D_2$};
\node[draw, black,ultra thick, inner sep=0pt,
text width=13.5mm,
align=center, circle] (d) at (2,-2) {$Y_2$};
\draw[->, -triangle 90] (e) edge (f) (h) edge (f) ;
\draw[->, -triangle 90, red] (h) edge (e) (e) edge (d) ;
\draw[->, -triangle 90, red] (h) edge[bend right=-30] node [left] {} (d);
\draw[->, -triangle 90, red, dotted] (f) edge (d) ;
\end{tikzpicture} }
\caption{The left panel illustrates all the possible causal paths under Sequential Ignorability (Assumption \ref{ass:seqign}). Here, past treatments may affect intermediate covariates, and future treatments may depend on past treatments, covariates and outcomes. The right panel presents two estimands of interest. In particular, $\mathrm{ATE}(\mathbf{1}, \mathbf{0})$ (the effect of increasing treatments in both periods) denotes the effect mediated through all red edges, including the dotted red edge. Instead, $\mathrm{ATE}((1, 0), (0, 0))$ (the \textit{direct} effect of only increasing treatment in the first period) denotes the effect mediated through all red edges excluding the dotted red edge. } \label{fig:seqign}
\end{figure}
\subsection{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.
\begin{ass} \label{ass:noant} For $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)$.
\end{ass}
Assumption \ref{ass:noant} 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{ass:noant} allows for anticipatory effects governed by \textit{expectations} (e.g., individuals may choose treatments based on \textit{expected} future utilities), but not on the future treatment \textit{realizations} \citep[see][for a discussion]{athey2022design}.
\begin{exmp}[Observed outcomes] \label{exmp:a_0}
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{ass:noant} holds. (Assumption \ref{ass:linearity} below will impose restrictions on $\mathbb{E}[g_2(\cdot)]$.)
\qed
\end{exmp}
In the rest of our discussion, we index potential outcomes and covariates by past treatment history under Assumption \ref{ass:noant}. 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 \textit{potential history} under treatment status $d_1$ in the first period. Here, $H_{i,2}$ can include interaction terms, omitted for brevity.
\begin{ass}[Sequential Ignorability] \label{ass:seqign} 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}
$$
\end{ass}
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{exmp:a_0}, Assumption \ref{ass:seqign} 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}.
$
\subsection{Potential local projections}
Following in spirit, \cite{jorda2005estimation}, we approximate the expectation of potential outcomes as linear functions of (high-dimensional) past characteristics. Different from \cite{jorda2005estimation}, linearity is imposed on expected potential instead of realized outcomes. Modeling potential outcomes directly avoids functional form restrictions on the treatment assignment mechanism.
\begin{ass} \label{ass:linearity}
For 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}
$$
\end{ass}
Assumption \ref{ass:linearity} 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 \citep{robins2000marginal}, Assumption \ref{ass:linearity} (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.
\begin{lem}[Identification] \label{lem:identification_model1} Let Assumptions \ref{ass:noant}, \ref{ass:seqign}, \ref{ass:linearity} 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}
$$
\end{lem}
The proof is in Appendix \ref{sec:lemma1}. Lemma \ref{lem:identification_model1} builds on results in the literature on marginal structural models \citep[e.g.][]{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{lem:identification_model1} motivates a recursive estimation strategy discussed in Section \ref{sec:two_p}.
\begin{exmp}[Linear Model] \label{exmp:sem1}
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{ass:linearity}.
\qed
\end{exmp}
\begin{rem}[Linearity in high-dimensions as an approximation to the true model] \label{rem:approximation} In the same spirit of \cite{belloni2014inference},
our results also directly extend to the case where we relax Assumption \ref{ass:linearity} 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 \textit{approximate} the conditional mean function as linear. As we further discuss in Section \ref{sec:theory}, 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
\end{rem}
\begin{rem}[Comparison with standard local projections and DiD] Appendix \ref{sec:lp1} presents an extensive discussion and comparison of Lemma \ref{lem:identification_model1} 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 \textit{observed} outcomes $Y_{i,T}$ instead of potential outcomes, which therefore would also depend on the distribution of treatment assignments.
\end{rem}
\section{Estimation with dynamic balancing} \label{sec:two_p}
\subsection{Estimation with two periods}
This section studies estimation. We defer to Section \ref{sec:app} a complete guide for practice, including discussion about the model, tuning parameters, and complexity. Appendix \ref{sec:extended_discussion} 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:
\begin{itemize}
\item \textit{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{lem:identification_model1}, 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.
\vspace{2mm}
\item \textit{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.
\vspace{2mm}
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 \citep[e.g.][]{tchetgen2012semiparametric, zhang2013robust, jiang2015doubly, nie2021learning}
\begin{equation} \label{eqn:myestimator}
\small
\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 \citep{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 \eqref{eqn:myestimator} uses regression adjustments over \textit{each} period to control the high-dimensional bias.
\vspace{2mm}
A choice of the weights from previous literature are inverse probability weights (IPW). These weights for the first and second period are
\begin{equation} \label{eqn:ipw}
\small
\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{rem:aipw}). 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.
\vspace{2mm}
\item \textit{Main balancing conditions:} The main step is in the choice of the balancing weights. The core idea is to decompose
\begin{equation} \label{eq:decomposition}
\small
\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 \small $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)})$, \normalsize we have
\begin{equation} \label{eqn:eqbalance}
\small
\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 \textit{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} \label{eqn:bb}
\small
\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 \eqref{eqn:bb} coincides with the balancing term in \cite{athey2018approximate}. However,
here we also require that histories in the second period are balanced, once \textit{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.
\vspace{2mm}
\item \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.
\vspace{2mm}
\item \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$.
\end{itemize}
We formalize this discussion below (Appendix \ref{lem:33} contains a formal proof).
\begin{thm}[Balancing weights] \label{thm:residual_final} \label{lem:balancing1} Let assumptions \ref{ass:noant} - \ref{ass:linearity} hold. Then Equation \eqref{eq:decomposition} holds with $T_1$ bounded as in Equation \eqref{eqn:eqbalance}.
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.
$
\end{thm}
\begin{rem}[Estimating when-to-treat policies] \label{rem:when_to_treat} 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., \cite{robins2000marginal}, \cite{imai2015robust}, \cite{acemoglu2019democracy}.
This differs from studying the effect of an optimal treatment path as in \cite{murphy2003optimal} or \cite{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 \eqref{eqn:H_it}, 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{sec:appendix_when_to_treat} discusses an extension for this setting.
\end{rem}
\subsection{Generalization to multiple time periods} \label{sec:multiple}
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.8, page = 1,
trim={1cm 17cm 1cm 2.5cm},clip]{./alg_final.pdf}
\end{figure}
We now describe in details the procedure with finite $T$ periods. Let $d_{1:T} = (d_1, \cdots, d_T)$,
\begin{equation}
\small
\begin{aligned}
\mathrm{ATE}(d_{1:T}, d_{1:T}') = \mu_T(d_{1:T}) - \mu_T(d_{1:T}'), \quad \mu_T(d_{1:T}) = \mathbb{E}\Big[Y_T(d_{1:T}) \Big].
\end{aligned}
\end{equation}
This estimand denotes the difference in potential outcomes for two treatment histories $d_{1:T}, d_{1:T}'$. We denote
\begin{equation} \label{eqn:H_it}
\small
\begin{aligned}
H_{i,t} = \Big[D_{i,1}, \cdots , D_{i,t - 1}, X_{i,1}, \cdots, X_{i,t}, Y_{i,1}, \cdots, Y_{i,t-1} \Big] \in \mathbb{R}^{p_t}
\end{aligned}
\end{equation}
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].
$
\begin{ass} \label{ass:seqignm} For any $d_{1:T} \in \{0,1\}^T$, and $t \le T$,
\begin{itemize}
\item[(A)] (No-anticipation) The potential history $H_{i,t}(d_{1:T})$ is constant in $d_{t:T}$;
\item[(B)] (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$;
\item[(C)] (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}
\end{ass}
Assumption \ref{ass:seqignm} generalizes Assumptions \ref{ass:noant}-\ref{ass:linearity} from the two-period setting. Identification follows similarly to Lemma \ref{lem:identification_model1}. For given weights $\hat{\gamma}_{1:T}$, and coefficients $\hat{\beta}^{(1:T)}$, we estimate
\begin{equation} \label{eqn:generalest}
\small
\begin{aligned}
\hat{\mu}_T(d_{1:T}) = & \sum_{i=1}^n \left\{\hat{\gamma}_{i,T}(d_{1:T}) Y_{i,T} - \sum_{t = 2}^T \Big(\hat{\gamma}_{i,t}(d_{1:T}) - \hat{\gamma}_{i,t-1}(d_{1:T})\Big) H_{i,t} \hat{\beta}_{d_{1:T}}^{(t)} - \Big(\hat{\gamma}_{i,1}(d_{1:T}) - \frac{1}{n} \Big) X_{i,1} \hat{\beta}_{d_{1:T}}^{(1)} \right\}.
\end{aligned}
\end{equation}
Coefficients are estimated recursively as in the two periods setting (see Algorithm 2).
\begin{lem} \label{lem:lemmam1}
Let $\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}
\small
\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}
\end{lem}
The proof is in Appendix \ref{app:lem22}.
Lemma \ref{lem:lemmam1} 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{lem:lemmam2}. Note that $\bar{X}_1 \beta_{d_{1:T}}^{(1)}$ converges to $\mu_T(d_{1:T})$ under standard $\sqrt{n}$-asymptotics.
\section{Theoretical properties and inference} \label{sec:theory}
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.
\begin{ass}[Overlap and tails' conditions] \label{ass:weakoverlap}
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.
\end{ass}
Condition (i) is the overlap condition, standard in the causal inference literature. The overlap condition is sufficient (but not necessary, see the discussion of \cite{athey2018approximate} in cross-sectional settings) to show existence of a feasible solution of Algorithm 1; see Remark \ref{rem:overlap}. Condition (ii) is a tail restriction. Assumption \ref{ass:weakoverlap} 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.
\begin{thm}[Existence of feasible weights] \label{thm:overlap}
Let Assumptions \ref{ass:seqignm}, \ref{ass:weakoverlap} 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}\label{eqn:hh}
\small
\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}
\end{thm}
The proof is in Appendix \ref{app:overlap_proof}. In particular, the reader may refer to Appendix Lemma \ref{lem:firststat} 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 \textit{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{thm:thm_asym_t} 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.
\begin{cor} \label{cor:comparison} Let the conditions in Theorem \ref{thm:overlap} 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 \eqref{eqn:hh}, and $n ||\hat{\gamma}_t||^2 \le n c_t ||\hat{\gamma}_{t-1}||^2$ for a finite constant $c_t < \infty$.
\end{cor}
Corollary \ref{cor:comparison} 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.
\begin{ass} \label{ass:4m} Let the following hold: for every $t \in \{1, \cdots, T\}$, $d_{1:T} \in \{0,1\}^T$,
\begin{itemize}
\item[(i)] $\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$;
\item[(ii)] 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.
\item[(iii)] $\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}
\end{ass}
Assumption \ref{ass:4m} 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., \cite{buhlmann2011statistics}. An example and derivation for condition (i) for Lasso under sparsity is included in Example \ref{lem:suff} (Appendix \ref{aa:suff}). As we discuss in Appendix Remark \ref{rem:restricted}, the restricted eigenvalue condition is imposed for all $H_t, t \ge 1$; references that study this condition from different angles include \cite{deshpande2023online}, Section C.4.
\begin{thm}[Parametric convergence rate] \label{thm:convergence_rate_first} Let the conditions in Theorem \ref{thm:overlap} and Assumption \ref{ass:4m} 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).
$
\end{thm}
Theorem \ref{thm:convergence_rate_first} guarantees a parametric convergence rate with high-dimensional covariates.
\begin{thm}[Inference] \label{thm:thm_asym_t}
Let the conditions in Theorem \ref{thm:overlap} and Assumption \ref{ass:4m} 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}
\small
\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}.
$$
\end{thm}
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{cor:ate}), as described in Appendix \ref{sec:main_theorem}. 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.
\begin{rem}[Strict overlap assumption] \label{rem:overlap} Strict overlap (Assumption \ref{ass:weakoverlap} $(i)$) is not necessary to achieve parametric convergence rates whenever a feasible solution to Algorithm 1 exists (see for instance \cite{athey2018approximate}, Lemma 2 in cross-sectional settings). \qed
\end{rem}
\begin{rem}[Rate conditions and comparison with AIPW] \label{rem:aipw} As noted in \cite{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).
$
\end{rem}
\begin{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{aa:suff}. 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{lem:firststat}). 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.
\end{rem}
\section{Guide to practice: numerical studies and application} \label{sec:app}
\subsection{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.
\textit{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
\begin{equation} \label{eqn:average_effect}
\small
\begin{aligned}
\mathbb{E}\left[Y_{i,T}\Big(D_{1:(T - h)}, d_{T- h + 1}, \cdots, d_T\Big)\right] -
\mathbb{E}\left[Y_{i,T}\Big(D_{1:(T - h)}, d_{T- h + 1}', \cdots, d_T'\Big)\right]
\end{aligned}
\end{equation}
for given treatment histories $d_{(T-h):T}, d_{(T-h):T}'$. Equation \eqref{eqn:average_effect} 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 \cite{imai2018matching}, the focus on Equation \eqref{eqn:average_effect} 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). \\
\textit{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 \cite{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
$$
\begin{aligned}
Y_{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},
\end{aligned}
$$
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.
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.8, page = 2,
trim={1cm 15cm 1cm 2.5cm},clip]{./alg_final.pdf}
\end{figure}
\begin{rem}[Tuning parameters] \label{rem:tuning} Similarly to one-dimensional setting \citep{athey2018approximate}, Algorithm 1 requires choosing tuning parameters. A complete description is in Algorithm \ref{alg:algtuning} 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{sec:theory} 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{sec:num} we illustrate the benefits of this procedure through numerical studies. \qed
\end{rem}
\begin{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{fig:long_T2} 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
\end{rem}
\subsection{Numerical studies} \label{sec:numerics}
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
\begin{equation} \label{eqn:propmodel}
\small
\begin{aligned}
\iota_{i,t} = \eta \sum_{s=1}^t X_{i,s} \phi+ \sum_{s=1}^{t-1} \delta_s (D_{i,s} - \bar{D}_s) + \xi_{i,t}, \quad \bar{D}_s =n^{-1} \sum_{i=1}^n D_{i,s}
\end{aligned}
\end{equation}
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 \cite{athey2018approximate}.
The larger $\eta$ corresponds to weaker overlap (see Table~\ref{tab:summaries} 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: \textit{Sparse} with $\beta^{(j)} \propto 1\{j \le 10\}$, \textit{Moderate} with moderately sparse $\beta^{(j)} \propto 1/j^2$ and the \textit{Harmonic} setting with $\beta^{(j)} \propto 1/j$. We set $\| \beta \|_2 =1, \tau = 1$.
We consider the following competing methodologies:
\begin{itemize}
\item \textit{Augmented IPW}, with \textit{known} propensity score and with \textit{estimated} propensity score. The method replaces the balancing weights in Equation \eqref{eqn:myestimator} 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) \citep[e.g.][for a discussion]{nie2021learning, bodoryevaluating}. For both AIPW and IPW we consider stabilized inverse probability weights.
\item \textit{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 \cite{yiu2018covariate, yiu2020joint} that, different from ours, require information about the propensity score. We follow Section 3 in \cite{yiu2020joint} for its implementation. (We do not also compare to \citealt{imai2015robust} for MSM since it is not feasible in high-dimensions.)
\item
\textit{``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 \cite{belloni2014inference}.
\item \textit{Naive Lasso}: it runs a regression controlling for covariates and treatment assignments.
\item \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.
\item \textit{DiD switchback}: it is a DiD estimator similar to \cite{de2024difference}.
\item \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.
\end{itemize}
For Dynamic Covariate Balancing, \textit{DCB}, the choice of tuning parameters is data adaptive, and it uses a grid-search method discussed in Appendix \ref{app:algorithms} and Remark \ref{rem:tuning}. 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{tab:mse2} 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{thm:overlap}, 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{sec:num} 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}
\begin{table*}[!ht]\centering
\caption{Mean Squared Error (MSE) for estimating the average treatment effect of always vs never being under treatment of Dynamic Covariate Balancing (DCB) across $200$ repetitions with sample size $400$ and $101$ variables in time period 1. This implies that the number of variables in time period $2$ and $3$ are $203$ and $304$. Oracle Estimator is denoted with aIPW$^*$ whereas aIPWh(l) denote AIPW with high(low)-dimensional estimated propensity. CAEW (MSM) corresponds to the method in \cite{yiu2020joint}, D.Lasso is adaptation of Double Lasso \citep{belloni2014inference}.} \label{tab:mse2}
\scalebox{0.7}{\begin{tabular}{@{}lrrrcrrrcrrr@{}}\toprule
& \multicolumn{3}{c}{$\eta=0.1$} & \ & \multicolumn{3}{c}{$\eta=0.3$} & \ & \multicolumn{3}{c}{$\eta=0.5$} \\
\cmidrule{2-4} \cmidrule{6-8} \cmidrule{10-12}
&\small sparse&\small mod &\small harm && \small sparse&\small mod & \small harm && \small sparse& \small mod & \small harm \\\Xhline{.8pt}
\rowcolor{lightgray}\multicolumn{12}{c}{$T=2$} \\ \Xhline{.8pt}
\small aIPW$^*$ & $0.069$ & $0.092$ & $0.071$ &&$0.102$ & $0.104$ & $0.118$ &&0.131&$0.127$ & $0.132$\\
\cline{2-4} \cline{6-8} \cline{10-12}
\rowcolor{glaucous!60} \small DCB & $0.060$ & ${0.077}$ & $0.075$&& ${0.092}$ & ${0.076}$ & ${0.084}$&&0.099&0.077&0.085\\
\small aIPWh & ${0.064}$ & $0.091$ & ${0.070}$&&$0.180$ & $0.204$ & $0.218$ &&$0.265$ & $0.312$ & $0.368$\\
\small aIPWl & $0.260$ & $0.229$ & $0.212$ && $0.157$ & $0.201$ & $0.165$&& $0.214$ & $0.234$ & $0.213$ \\
\small IPWh & $2.37$ & $1.78$ & $2.80$&&$10.19$ & $6.49$ & $11.72$ &&$15.25$ & $8.09$ & $16.67$\\
\small Seq.Est. & $0.932$ & $1.333$ & $0.692$ && $1.388$ & $1.787$ & $1.152$&& $1.759$ & $1.795$ & $1.664$ \\
\small Lasso & $0.247$ & $0.410$ & $0.132$ && $0.509$ & $0.710$ & $0.298$ && $0.762$ & $0.948$ & $0.560$ \\
\small CAEW & $0.432$ & $0.444$ & $0.517$ &&$1.934$ & $1.274$ & $1.974$&& $3.376$ & $2.168$ & $4.423$ \\
\small Dyn.D.Lasso & $0.124$ & $0.118$ & $0.256$ && $0.208$ & $0.147$ & $0.430$&& $0.218$ & $0.153$ & $0.554$ \\
\small DiD Switchback & $2.06$ & $1.71$ & $1.60$ && $14.52$ & $6.98$ & $20.72$ && $38.79$ & $6.98$ & $52.48$ \\
\small Simple LP & $1.28$ & $1.40$ & $1.37$ && $1.613$ & $1.682$ & $1.573$ && $1.777$ & $1.689$ & $1.808$ \\
\Xhline{.8pt} \rowcolor{lightgray}\multicolumn{12}{c}{$T=3$} \\ \Xhline{.8pt}
\small aIPW$^*$ & $0.226$ & $0.296$ & $0.261$&& $0.403$ & $0.251$ & $0.339$&& $0.472$ & $0.496$ & $0.562$\\
\cline{2-4} \cline{6-8} \cline{10-12}
\rowcolor{glaucous!60} \small DCB & ${0.155}$ & ${0.208}$ & ${0.199}$ && ${0.257}$ & ${0.217}$ & ${0.329}$ && $0.294$ & $0.267$ & $ 0.455$\\
\small aIPWh & $0.201$ & $0.273$ & $0.280$&& $0.595$ & $0.747$ & $0.835$&& $0.999$ & $1.328$ & $1.607$ \\
\small aIPWl & $0.823$ & $0.625$ & $0.829$ &&$0.623$ & $0.704$ & $0.638$&& $1.078$ & $1.396$ & $1.234$ \\
\small IPWh & $11.03$ & $8.09$ & $12.84$&&$34.65$ & $20.34$ & $39.37$ &&$47.65$ & $23.30$ & $45.47$\\
\small Seq.Est. & $2.608$ & $4.016$ & $2.316$ && $3.722$ & $5.269$ & $3.818$&& $5.279$ & $6.829$ & $5.467$ \\
\small Lasso & $0.409$ & $0.492$ & $0.514$ && $0.559$ & $0.732$ & $0.507$ && $1.290$ & $1.315$ & $1.174$ \\
\small CAEW & $3.580$ & $2.446$ & $4.279$ && $18.50$ & $12.07$ & $22.85$ && $30.07$ & $18.71$ & $33.01$ \\
\small Dyn.D.Lasso & $0.471$ & $0.344$ & $0.679$ && $0.694$ & $0.378$ & $1.182$ && $0.964$ & $0.383$ & $1.594$ \\
\small DiD Switchback & $24.9$ & $27.98$ & $20.52$ && $21.07$ & $7.503$ & $38.33$ && $59.23$ & $22.15$ & $89.07$ \\
\small Simple LP & $7.38$ & $7.54$ & $7.38$ && $8.180$ & $8.154$ & $8.294$ && $9.131$ & $9.087$ & $9.217$ \\
\bottomrule
\end{tabular} }
\end{table*}
\subsection{Empirical illustration}
In this section, we present an empirical application for studying the effect of democracy on GDP growth using data from \cite{acemoglu2019democracy}. \cite{acemoglu2019democracy} studied dynamic treatment effects of democracy under GDP growth under sequential ignorability \citep[Assumption 1 in][]{acemoglu2019democracy}. Figure \ref{fig:treatment_path2} 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 \cite{acemoglu2019democracy}.
We use the same treatment specification as in \cite{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 \cite{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 \cite{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}
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.7]{./figures/treatment_status-eps-converted-to.pdf}
\caption{The figure illustrates the dynamics of treatments. }
\label{fig:treatment_path2}
\end{figure}
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.35]{./figures/overlap_acemoglu-eps-converted-to.pdf}
\caption{Estimated probability of treatment for one year (left-panel) and two consecutive years (right-panel). Estimation is performed via logistic regression with a pooled regression with year, region fixed effects and four lagged outcomes. The right panel also controls for the past treatment assignment. The figure illustrates the sensitivity of inverse probability weights to longer time horizon, motivating more stable balancing weights.} \label{fig:overlap}
\end{figure}
\begin{table}[!htbp]
\centering
\caption{Diagnostics across horizons with time horizon $h= 3$: inverse of weights dispersion for DCB weights and $1/||w_t||^2$ for IPW where $w_t$ is the stabilized inverse probability weight. Estimation of IPW is performed via logistic regression with a pooled regression with year, region fixed effects and four lagged outcomes. Larger number indicates better precision.}
\label{tab:diagnostics_h}
\begin{tabular}{lccc|ccc}
\toprule
$h = 3$ & \multicolumn{3}{c|}{$\mathbb{E}[Y_T(\mathbf{1}_3, D_{T-3:1})]$} & \multicolumn{3}{c}{$\mathbb{E}[Y_T(\mathbf{0}_3, D_{(T-3):1}))]$} \\
\cmidrule(lr){2-4}\cmidrule(lr){5-7}
& $1/||\hat{\gamma}_1||^2$ & $1/||\hat{\gamma}_2||^2$ & $1/||\hat{\gamma}_3||^2$ & $1/||w_1||^2$ & $1/||w_2||^2$ & $1/||w_3||^2$ \\
\midrule
DCB & 67 & 64 & 62 & 34 & 34 & 34 \\
IPW & 73 & 17 & 10 & 32 & 23 & 23 \\
\bottomrule
\end{tabular}
\end{table}
Figure \ref{fig:acemoglu} 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 \cite{acemoglu2019democracy}, and results are robust across the two specifications for DCB.
We compare our method to (i) the \textit{linear estimator} reported by \cite{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 \textit{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 \textit{future} assignments. Therefore, the causal effects estimated by the local projection differ from the target long-run effect, which instead \textit{fixes} future treatment assignments.
The effect estimated as in \cite{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 \cite{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 \cite{acemoglu2019democracy} with and without unit fixed effects report almost identical results.
To investigate differences with (A)IPW methods, the right panel in Figure \ref{fig:acemoglu} presents comparisons in terms of the imbalance over the lagged outcome at time $t - 1$ when using balancing or inverse probability weights. As \cite{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{fig:overlap} complements Figure \ref{fig:acemoglu} showing instability of inverse probability weights, and Figure \ref{fig:dispersion} in the Appendix show that DCB weights present less dispersion than IPW weights. As a helpful diagnostic, in Table \ref{tab:diagnostics_h} 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.
\begin{figure}[!ht]
\centering
\includegraphics[scale=0.35]{./figures/main_acemoglu-eps-converted-to.pdf}
\caption{Left-hand side: pooled regression from $t \in \{1989, \cdots, 2010\}$. Gray region denotes the $90\%$ confidence band for the least parsimonious model. DCB and DCB2 refer to two separate specification, with DCB corresponding to the more parsimonious one. LP denotes a local projection on $t$ periods before. (no fe) indicates a specification without country fixed effects. Right-hand side reports $\log(|I(dcb)| + 1)/\log(|I(aipw)| + 1) - 1 \approx |I(dcb)|/|I(aipw)| - 1$, where $I(\cdot)$ denotes the imbalance (as in Lemma \ref{lem:balancing1}) in the lagged outcome usin DCB or IPW weights. } \label{fig:acemoglu}
\end{figure}