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.
168,020 characters
Causal Inference under Dynamic Selection: Time-Varying Covariates and Latent Heterogeneity
\title{Causal Inference under Dynamic Selection: Time-Varying Covariates and Latent Heterogeneity \thanks{
I thank Ben Deaner, Aureo de Paula, Chen-Wei Hsiang, Ge Sun, Zihao Wang, Andrei Zeleneev, and the participants in the UCL Econometrics Brownbag Seminar for their valuable comments.
}}
\author{ Weisheng \textsc{Zhang}
\thanks{
University College London: \textsf{[email removed]}.}
}
\maketitle
\begin{abstract}
I study dynamic treatment effects in panel data under staggered adoption when treatment timing depends jointly on unobserved time-invariant heterogeneity and time-varying pretreatment covariates, including lagged outcomes. Untreated potential outcomes follow a nonparametric dynamic panel model that allows flexible interactions between time-varying covariates and latent heterogeneity. I use pretreatment outcome histories to find individuals with similar time-invariant latent factors, and the key requirement is that these histories are sufficiently informative about those latent factors. I develop an identification strategy for the dynamic average treatment effect on the treated (ATT) and propose kernel-based doubly robust estimators for the dynamic ATT. I further combine double cross-fitting with undersmoothing and show that, under suitable regularity conditions, the proposed estimators are $\sqrt{n}$-consistent, asymptotically normal, and asymptotically unbiased. The simulation study demonstrates that the proposed method provides accurate inference across a wide range of data-generating processes. I illustrate the method with an application to the U.S. family planning program studied by \citet{bailey2012reexamining} and reestimate its effect on fertility rates.
\end{abstract}
\thispagestyle{empty}
\clearpage
\setcounter{page}{1}
\clearpage
\section{Introduction}\label{sec:introduction}
Self-selection is a fundamental challenge in causal inference, where treated and untreated units may be systematically different. Panel data can help address this problem because they contain repeated observations of the same units over time and allow researchers to control for unobserved time-invariant heterogeneity (fixed effects). Accordingly, causal panel methods, including difference-in-differences (DiD), synthetic control, and matrix completion, allow for selection on unobserved time-invariant factors and have become increasingly popular in empirical research.
In many empirical applications, treated and untreated units differ not only in time-invariant characteristics but also in \emph{time-varying} covariates, such as lagged outcomes. For example, participants in job training programs may differ from nonparticipants in ability and may also experience a decline in earnings before training \citep{ashenfelter1985using}. Failing to control for pretreatment earnings may therefore produce a spurious positive estimate of the program's effect on wages. As another example, counties with higher recent fertility rates may have greater demand for subsidized contraceptive services and therefore be more likely to apply for federal funding. In this case, controlling for pretreatment time-varying covariates is necessary to disentangle the treatment effect from the dynamic effects arising from lagged outcomes.
However, the validity of most existing causal panel methods builds on the assumption that treatment selection depends \emph{only} on time-invariant characteristics. For example, \citet{ghanem2022selection} and \citet{marx2024parallel} show that, except in a few special cases, the parallel trends assumption in DiD fails when treatment depends on both time-invariant characteristics and lagged outcomes. Other causal panel methods face similar limitations when treatment assignment also depends on time-varying pretreatment covariates \citep{arkhangelsky2024causal}.
In this paper, I propose a new method for causal inference with panel data when treatment selection depends jointly on time-invariant latent factors and time-varying pretreatment covariates. I establish nonparametric identification of the dynamic average treatment effect on the treated (ATT), develop a doubly robust estimator, and provide a corresponding inference procedure. The method is designed to accommodate staggered adoption, which is common in empirical applications.
The proposed method allows the outcome process to depend on time-varying covariates, including lagged outcomes. Such dynamic dependence is ubiquitous in economics (e.g., due to intertemporal decision-making, adjustment costs, and persistent shocks) but is often abstracted away in causal panel methods~\citep{arkhangelsky2024causal}. It is important to accommodate dynamic dependence when lagged outcomes affect both treatment assignment and future outcomes. Otherwise, researchers may fail to distinguish treatment effects from dynamic effects driven by lagged outcomes.
In addition, I do not impose functional-form restrictions on untreated potential outcomes. Many popular causal panel methods rely on linearity assumptions that may be too restrictive in applications. For example, DiD relies on the parallel trends assumption, and most synthetic control and matrix completion methods require that untreated potential outcomes follow a linear factor model. In contrast, my method builds on a nonparametric model in which outcomes depend flexibly on their histories and time-invariant latent factors, eliminating concerns about the functional-form specification of untreated potential outcomes in empirical applications.
The key technical challenge is that time-invariant latent factors are unobserved. Since untreated outcomes follow a possibly nonlinear model, these factors can neither be differenced out as in DiD nor estimated using standard linear factor methods. \citet{feng2023optimal}, \citet{deaner2025inferring}, and \citet{athey2025identification} address this challenge in a pure nonlinear factor model without dynamic dependence. They construct pairwise pseudo-distances from long pretreatment outcome histories. Each pseudo-distance measures the similarity between two individuals' pretreatment histories and serves as a proxy for the similarity between pairwise latent factors. However, with dynamic dependence, outcome histories reflect both the effect of time-invariant heterogeneity and the dynamic effects of time-varying covariates, so their methods no longer apply directly.
To address dynamic dependence, I establish that the pairwise pseudo-distance remains informative about time-invariant latent factors in dynamic settings under weak dependence. This extends the applicability of the pseudo-distance approach to dynamic panel models. The intuition is that the dynamic effects decay over time, whereas individual latent factors continue to affect the entire outcome path. As a result, individuals with similar latent factors exhibit similar long-run pretreatment histories. Since the pseudo-distance measures similarity between pretreatment histories, it can serve as a proxy for latent similarity.
The identification of the dynamic ATT uses matching based on time-varying covariates and the pseudo-distance. The key assumption for identification is the informativeness condition, which requires that if the pseudo-distance between two individuals is small, then the distance between their latent factors is also small. Under this assumption, although time-invariant latent factors are unobserved, individuals with similar latent factors are identified by pseudo-distance. I then match treated units to untreated units with similar covariates and a small pseudo-distance, impute untreated counterfactual outcomes using the matched units, and identify the dynamic ATT. Under staggered adoption, this matching strategy is combined with a backward recursive procedure to accommodate dynamic selection. The identification result is new under selection on both time-invariant latent heterogeneity and time-varying covariates and requires only pairwise latent similarity rather than identification of the latent factors or their distribution.
I propose a doubly robust estimator for the dynamic ATT. I use Nadaraya-Watson estimators based on both covariates and the sample pseudo-distance to nonparametrically estimate counterfactual outcomes and propensity scores, thereby accounting for selection on observed time-varying covariates and time-invariant latent factors. When pretreatment histories are sufficiently long, these estimators achieve the same convergence rates as if time-invariant latent factors were observed.
To obtain valid inference in a broader range of settings, I apply double cross-fitting to the proposed estimator. Unlike standard cross-fitting in \citet{chernozhukov2018double}, double cross-fitting constructs the counterfactual outcome and propensity score estimators using separate training samples. This reduces the dependence between their estimation errors and, when combined with undersmoothing, yields a faster convergence rate for the proposed estimator. I establish the root-$n$ consistency, asymptotic normality, and asymptotic unbiasedness of the proposed estimator, and show that these properties hold under more general conditions than those required for standard cross-fitting.
I provide extensive simulation evidence to evaluate the finite-sample performance of the proposed estimator. The estimator performs well in terms of bias, standard deviation, and coverage across different data-generating processes and sample sizes, even when the panel length is moderate. Based on these simulation results, practical guidance is provided for bandwidth selection for applied researchers.
I then apply the proposed method to the U.S. family planning programs studied by \citet{bailey2012reexamining} and reestimate their effects on fertility rates. First, I find strong evidence that counties with higher lagged fertility rates are more likely to be treated, which motivates the use of the proposed method. Second, my method estimates statistically significant and persistent reductions in fertility. This finding is consistent with the results in \citet{bailey2012reexamining} and those obtained using DiD. Moreover, the estimated reductions are larger than the DiD estimates, and the gap widens over longer post-treatment horizons. This difference reflects the importance of accounting for dynamic selection based on lagged fertility rates.
\paragraph{Literature review} This paper contributes to the large and rapidly growing literature on causal inference with panel data by allowing treatment selection to depend on time-varying covariates. I refer readers to \citet{de2023two} and \citet{roth2023s} for recent surveys of the DiD literature and to \citet{arkhangelsky2024causal} for a broader survey of causal inference with panel data. The implications of treatment selection for the validity of causal panel methods have received increasing attention. Earlier discussions of this issue in DiD include \citet{ashenfelter1985using} and \citet{abadie2005semiparametric}, and more recently, \citet{ghanem2022selection} and \citet{marx2024parallel} show that, except in a few special cases, the parallel trends assumption fails when treatment decisions respond to previous outcomes. \citet{arkhangelsky2024causal} further emphasize that many causal panel results are established under strict exogeneity, which rules out treatment selection based on time-varying outcome shocks. Although \citet{arkhangelsky2023large} show the restriction can be relaxed in synthetic control in simultaneous-treatment settings, their arguments do not directly extend to staggered adoption.
A related strand of the causal panel literature addresses treatment selection on time-varying covariates. For DiD, these methods rely on conditional parallel trends \citep{abadie2005semiparametric,wooldridge2025two}, with extensions to doubly robust estimation \citep{sant2020doubly}, high-dimensional covariates \citep{chang2020double}, and treatment-affected covariates \citep{caetano2022difference}. However, conditional parallel trends may fail when treatment selection depends jointly on observed covariates and time-invariant latent factors. In addition, for dynamic treatment effects, \citet{lewis2020double} and \citet{viviano2026dynamic} allow selection on high-dimensional covariates but assume that treatment assignment is independent of time-invariant latent factors conditional on these covariates. \citet{marx2024heterogeneous} and \citet{botosaru2025time} focus on dynamic treatment effects under linear panel models. In contrast, the proposed approach uses a long pretreatment panel to accommodate selection on both time-varying covariates and time-invariant latent factors without imposing functional-form restrictions on either the propensity score or the untreated outcome process.
Lastly, \citet{gulek2025synthetic} combine synthetic control and instrumental variables to address unmeasured time-varying confounding, while my analysis focuses on treatment selection based on observed time-varying covariates and does not require an external instrument.
Beyond the literature on causal panel methods, this paper connects to the broader econometric literature on panel-data models. This literature has paid close attention to individual heterogeneity and outcome dynamics in both linear settings \citep[e.g.,][]{arellano1991some,blundell1998initial,alvarez2003time} and nonlinear settings \citep[e.g.,][]{blundell2002individual,honore2000panel,bonhomme2012functional, bonhomme2015grouped, bonhomme2022discretizing}, but has primarily focused on estimating parametric models. By contrast, the causal panel literature allows for rich heterogeneity in treatment effects but often abstracts away from dynamic dependence in outcomes \citep{arkhangelsky2024causal}. The proposed method bridges these two literatures by allowing for both flexible interactions between time-varying covariates and latent heterogeneity and dynamic dependence in untreated potential outcomes, without imposing a parametric outcome model.
The nonparametric model for untreated potential outcomes connects this paper to the nonparametric panel literature. For example, \citet{chernozhukov2013average}, \citet{hoderlein2012nonparametric}, and \citet{chernozhukov2026linear} study average structural and causal effects under time-homogeneity conditions. The proposed method allows for rich time heterogeneity in untreated potential outcomes without imposing time-homogeneity conditions. Another strand of literature uses observations from different periods as measurements or proxies for (possibly time-varying) latent confounders \citep[e.g.,][]{hu2012nonparametric,sasaki2015heterogeneity, deaner2018proxy}.
These methods only need short panels but rely on completeness conditions. In comparison, although the proposed method requires a long panel, it imposes no completeness conditions and avoids solving an inverse problem.
\bigskip
Methodologically, this paper builds on the recent literature on nonlinear factor models for causal inference with panel data, particularly \citet{deaner2025inferring}.
The pseudo-distance approach to measuring similarity between latent factors was first developed by~\citet{zhang2017estimating} for graphon estimation. Since then, different forms of the pseudo-distance have been used in a range of applications, including nonparametric graphon estimation \citep{zeleneev2020identification,nowakowicz2024nonparametric}, network imputation \citep{sun2026flexible}, controlling for unobserved heterogeneity using network data \citep{auerbach2022identification}, nonlinear factor model estimation \citep{feng2023optimal}, and causal inference \citep{feng2020causal, deaner2025inferring,athey2025identification,hoshino2024estimating}. However, the validity of most pseudo-distance methods is established only under pure factor models. In panel settings, restricting attention to such models rules out dynamic dependence. I extend the scope of the pseudo-distance approach to dynamic models in which outcomes depend on both latent factors and lagged outcomes. This extension is new and is of independent interest for the study of dynamic nonseparable panel models beyond causal inference.
The proposed estimator builds on the doubly robust estimation of the ATT. Doubly robust estimation combines Neyman-orthogonal moments with cross-fitting to obtain valid inference on target parameters in the presence of nonparametric nuisance parameters \citep{chernozhukov2018double}.
In the causal panel literature, doubly robust methods have been developed for DiD and TWFE models \citep{arkhangelsky2022doubly,arkhangelsky2024design,sant2020doubly}, dynamic treatment effects \citep{chernozhukov2022automatic,lewis2020double}, and linear factor models \citep{abadie2024doubly}. Recently, \citet{feng2020causal} and \citet{deaner2025inferring} extend doubly robust methods to nonlinear factor models. The proposed estimator builds directly on the doubly robust framework of \citet{deaner2025inferring} and likewise requires only weak smoothness conditions.
The double cross-fitting procedure was introduced by \citet{newey2018cross} and further developed by \citet{mcclean2026double}. I extend this method to dynamic ATT estimation with nonparametric regressions on \emph{latent factors}. Combining double cross-fitting and undersmoothing extends the inference theory in~\citet{deaner2025inferring} and yields valid inference under weaker dimensionality restrictions. This provides a new inference result for doubly robust estimation in the causal panel literature.
\bigskip
The rest of the paper is organized as follows. Section~\ref{sec:selection} presents the basic setup and the selection mechanism. Section~\ref{sec:overview_method} provides an overview of the proposed method. Section~\ref{sec:identification} establishes identification. Section~\ref{sec:estimation} derives the asymptotic distribution of the proposed estimators. Section~\ref{sec:extension} extends the analysis to staggered adoption. Section~\ref{sec:simulation} provides simulation evidence, and Section~\ref{sec:empirical} applies the method to the setting studied by \citet{bailey2012reexamining} to reexamine the effects of U.S. family planning programs on fertility rates.
\section{Selection and Model}\label{sec:selection}
This section introduces the basic setup of this paper. I discuss the treatment-selection problem in panel-data causal inference, focusing on settings where treatment assignment depends on both pretreatment information and latent heterogeneity. I then introduce a general outcome model that accommodates individual and time fixed effects and dynamic dependence. To fix ideas, I first focus on simultaneous-treatment settings (also called block assignment), and extend the analysis to staggered adoption in Section~\ref{sec:extension}.
\subsection{Data}
The panel data ${(Y_{it},X_{it},D_i,\alpha_i)}$ are indexed by $i=1,\ldots,n$ and $t=0,1,\ldots,T$. Here, $Y_{it}\in\mathbb{R}$ denotes the outcome, and $X_{it}\in\mathbb{R}^{d_X}$ denotes a vector of time-varying covariates that may affect treatment selection. An important case is when $X_{it}$ includes lagged outcomes, e.g., $X_{it}:=Y_{it-1}$. Treatment status is denoted by $D_i\in \{0,1\}$, and multidimensional time-invariant latent heterogeneity is captured by $\alpha_i\in\mathcal{A}\subset\mathbb{R}^{d_\alpha}$. I consider a simultaneous-treatment setting in which treatment is assigned only once, at time $T_0$. Treatment is absorbing, so once a unit is treated at $T_0$, it remains treated thereafter. I assume that $X_{iT_0}$ is realized before treatment assignment at $T_0$ and therefore constitutes pretreatment information.
Let $Y_{it}(0)$ and $Y_{it}(1)$ denote the untreated and treated potential outcomes, respectively. Because treatment occurs only at $T_0$, all units are untreated before that date. Then, under the non-anticipation condition,
\begin{align*}
Y_{it}=
\begin{cases}
Y_{it}(0), & t = 0, \ldots, T_0 - 1,\\
D_iY_{it}(1)+(1-D_i)Y_{it}(0), & t = T_0, \ldots, T.
\end{cases}
\end{align*}
I also allow for an unobserved time factor $\gamma_t$, which captures aggregate shocks, such as macroeconomic conditions, that affect all units in the same period. For notational simplicity, let $\Gamma_t := (\ldots,\gamma_{t-1},\gamma_t)$ denote the realized time-factor path up to time $t$.
In addition, define $\mathbb{E}_T(\cdot):=\mathbb{E}(\cdot\mid\Gamma_T)$ and $\mathbb{P}_T(\cdot):=\mathbb{P}(\cdot\mid\Gamma_T)$ as the expectation and probability conditional on the realized time-factor path $\Gamma_T$.
For each $t = T_0, \ldots, T$, the estimand of interest is the dynamic ATT
\footnote{
The definition differs from the standard ATT, $\mathbb{E}\left(Y_{it}(1)-Y_{it}(0)\mid D_i=1\right)$, only because I condition explicitly on the realized time fixed effects. Since these effects are fixed over the sample period, this is purely notational.
}:
\begin{align*}
\mathrm{ATT}(t) := \mathbb{E}_T\left(Y_{it}(1) - Y_{it}(0) \mid D_i = 1\right).
\end{align*}
\subsection{Selection Mechanism}
Self-selection is common in economic analyses using observational data, where treated and untreated units often differ systematically, and failing to account for such selection can lead to biased estimates of causal effects. One common source of self-selection is unobserved time-invariant characteristics (fixed effects). A classic example is that individuals with higher unobserved ability may obtain more education and also earn higher wages, leading to an overestimate of the return to education. Panel data provide a way to account for selection on unobserved time-invariant characteristics because they contain repeated observations for the same units. The recent and rapidly growing literature on causal panel methods, including DiD, synthetic control, and matrix completion methods, builds on this feature to control for selection on fixed effects.
In many empirical applications, treated and untreated units differ not only in time-invariant characteristics but also in observed \emph{time-varying} characteristics, such as lagged outcomes. A well-known example is Ashenfelter's dip, in which participants in job training programs typically experience a decline in earnings prior to training \citep{ashenfelter1978estimating, ashenfelter1985using}. Failing to control for this pre-treatment earnings decline may lead to a positive estimated treatment effect on wages even when the job training program has no effect. Similar selection arises in local policy adoption. For example, counties with higher recent fertility rates may have greater demand for federally subsidized contraceptive services and may therefore be more likely to apply for the program. In this case, controlling for selection on pretreatment covariates is necessary to disentangle the treatment effect from the dynamic effects arising from lagged outcomes.
However, the validity of many existing causal panel methods usually builds on the assumption that treatment selection depends only on time-invariant characteristics. For DiD, it has been established that, except in a few special cases, the parallel trends assumption fails when treatment decisions respond to both time-invariant characteristics and previous outcomes \citep{ghanem2022selection, marx2024parallel}. Similar issues arise for other causal panel methods, whose identification assumptions cannot accommodate treatment selection that depends jointly on time-invariant characteristics and time-varying pretreatment covariates \citep{arkhangelsky2024causal}. Although \citet{arkhangelsky2023large} show that this issue is less relevant for synthetic control under simultaneous treatment, their argument does not directly extend to staggered adoption.
A common approach in the DiD literature to allow selection on observed pretreatment covariates is to impose the conditional parallel trends assumption \citep{abadie2005semiparametric, sant2020doubly, chang2020double, caetano2022difference,wooldridge2025two}. This requires that, before and after treatment time $T_0$, treated and untreated units should have parallel trends if they have the same covariates $X_{iT_0}$, i.e.,
\begin{align*}
\mathbb{E}\left( Y_{iT_0 }(0) - Y_{iT_0 - 1}(0)\mid D_i = 1, X_{iT_0}\right) = \mathbb{E}\left( Y_{iT_0}(0) - Y_{iT_0 - 1}(0)\mid D_i = 0, X_{iT_0}\right).
\end{align*}
However, conditional parallel trends may still fail when treatment selection depends jointly on time-invariant latent factors and pretreatment covariates, even if untreated potential outcomes admit an additive two-way fixed-effects form.
Figure~\ref{fig:motivation} presents a simple simulation to illustrate how violations of standard and conditional parallel trends lead to biased estimators even when $Y_{it}(0)$ admits an additive two-way fixed-effects form. When treatment selection, $D_i \sim \mathrm{Logit}(\alpha_i + Y_{iT_0-1})$, depends jointly on fixed effects and lagged outcomes, the distribution of the standard DiD estimator is centered well below the true ATT (Figure~\ref{fig:motivation}(A)). In addition, the distribution of the doubly robust DiD estimator proposed by \citet{sant2020doubly}, which relies on conditional parallel trends, is also heavily biased (Figure~\ref{fig:motivation}(B)). The substantial bias is distinct from the negative-weight problems emphasized in the recent literature\footnote{
The simulation is conducted under block assignment with homogeneous treatment effects, making standard DiD equivalent to the recently proposed DiD estimators \citep{callaway2021difference,sun2021estimating,borusyak2024revisiting}.
}. Instead, it arises because joint selection on latent heterogeneity and lagged outcomes violates both standard and conditional parallel trends.
\begin{figure}[H]
\centering
\includegraphics[width=1.0\textwidth]
{figures/motivation.pdf}
\caption{Finite-sample distributions of the proposed estimator, the standard DiD, and the doubly robust DiD estimator, with $N = 1000$ and $T_0 = 20$. }
\label{fig:motivation}
\end{figure}
I now introduce the following assumption to formalize the selection mechanism: treatment is as good as randomly assigned conditional on both time-invariant heterogeneity and time-varying covariates.
\begin{assumption}[Selection]\label{assumption:selection_simultaneous}
For each $t = T_0, \ldots, T$,
\begin{align}
Y_{it}(0) \perp D_i
\mid
X_{iT_0}, \alpha_i, \Gamma_T.
\end{align}
\end{assumption}
Assumption~\ref{assumption:selection_simultaneous} requires that, after conditioning on pretreatment covariates $X_{iT_0}$ and time-invariant heterogeneity $\alpha_i$, treatment assignment is independent of the untreated potential outcome $Y_{it}(0)$ for each $t=T_0,\ldots,T$. I also condition explicitly on $\Gamma_T$ to maintain the fixed-effect interpretation of the time factors. The conditional independence is imposed only on untreated potential outcomes, rather than on treated potential outcomes, because the estimand of interest is the ATT and only untreated counterfactual outcomes for treated individuals need to be imputed.
Assumption~\ref{assumption:selection_simultaneous} extends the conventional selection-on-fixed-effects assumption
by additionally allowing treatment selection to depend on time-varying pretreatment covariates. In addition, this assumption nests the sequential unconfoundedness assumption \citep{robins2000marginal,viviano2026dynamic,marx2024heterogeneous}, $Y_{it}(0) \perp D_i \mid X_{iT_0}$, as a special case in which treatment selection does not depend on unobserved fixed effects. Lastly, Assumption~\ref{assumption:selection_simultaneous} arises naturally in dynamic economic models \citep{heckman2007dynamic}, in which a forward-looking agent chooses treatment to maximize expected utility based on the state variables $(X_{iT_0},\alpha_i)$.
\begin{remark*}
Because $X_{iT_0}$ is realized before treatment assignment, it cannot include potential outcomes at or after treatment assignment, such as $(Y_{iT_0}(0),Y_{iT_0}(1))$. This timing restriction therefore rules out anticipation effects (see the discussion in~\citet{roth2023s}) and Roy-type selection. In addition, Assumption~\ref{assumption:selection_simultaneous} rules out unobserved time-varying confounders that jointly affect treatment assignment and untreated potential outcomes. Addressing such confounding typically requires additional information, such as an external instrument. I refer readers to \citet{gulek2025synthetic}, who combine synthetic control and instrumental variables to address such confounding.
\end{remark*}
For any $t = T_0, \ldots, T$, define the regression outcome and propensity score as:
\begin{equation}\label{eq:expected_outcome_ps_simultaneous}
\begin{gathered}
m_{t}(x , \alpha ) := \mathbb{E}_T(Y_{it}(0) \mid X_{iT_0} = x, \alpha_i = \alpha ), \quad
\pi(x, \alpha ) := \mathbb{P}_T(D_i = 1 \mid X_{iT_0} = x, \alpha_i = \alpha ).
\end{gathered}
\end{equation}
Here, the function $m_t:\mathbb{R}^{d_{X}}\times\mathbb{R}^{d_{\alpha}}\mapsto\mathbb{R}$ denotes the conditional expectation of the untreated outcome at time $t$ given pretreatment information, and the function $\pi:\mathbb{R}^{d_X}\times\mathbb{R}^{d_{\alpha}}\mapsto [0,1]$ denotes the propensity score given individuals' pretreatment information. I write
\begin{align*}
Y_{it}(0) = m_{t}(X_{iT_0} , \alpha_i ) + u_{it}, \quad D_{i} = \pi (X_{iT_0} , \alpha_i ) + e_{i},
\end{align*}
where $u_{it}$ and $e_i$ are projection errors satisfying $\mathbb{E}_T(u_{it} \mid X_{iT_0}, \alpha_i ) = \mathbb{E}_T(e_{i} \mid X_{iT_0}, \alpha_i ) = 0$. It is straightforward to check that $\mathbb{E}_T(u_{it}e_i\mid X_{iT_0},\alpha_i)=0$ by selection assumption. The following assumption imposes a mean-independent Markov restriction, under which earlier histories have no additional predictive power for future untreated outcomes or treatment assignment when conditioning on $(X_{iT_0},\alpha_i)$.
\begin{assumption}[Mean-independent Markov]\label{assumption:surrogacy_simultaneous}
Let $H_{it} := \{(Y_{i\tau}, X_{i\tau })\}_{\tau \leq t}$ denote the history up to $t$. Then, for any $t = T_0, \ldots, T$,
\begin{enumerate}[label=(\roman*)]
\item \label{item:assumption_surrogacy_simultaneous_outcome}$\mathbb{E}_T (u_{it} \mid H_{iT_0-1}, X_{iT_0}, \alpha_i ) = \mathbb{E}_T (u_{it} \mid X_{iT_0}, \alpha_i ) = 0$;
\item \label{item:assumption_surrogacy_simultaneous_propensity_score} $\mathbb{E}_T (e_i \mid H_{iT_0-1}, X_{iT_0}, \alpha_i ) = \mathbb{E}_T (e_i \mid X_{iT_0}, \alpha_i ) = 0$.
\item \label{item:assumption_surrogacy_simultaneous_product} $\mathbb{E}_T (u_{it}e_i \mid H_{iT_0-1}, X_{iT_0}, \alpha_i ) = \mathbb{E}_T (u_{it}e_i \mid X_{iT_0}, \alpha_i ) = 0$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:surrogacy_simultaneous}\ref{item:assumption_surrogacy_simultaneous_outcome} and \ref{item:assumption_surrogacy_simultaneous_propensity_score} state that, conditional on the time-invariant factor $\alpha_i$ and the covariates $X_{iT_0}$, the future untreated path and treatment assignment are mean independent of earlier histories $H_{iT_0-1}$. Specifically, when $X_{iT_0}:=(Y_{iT_0-1}(0),\ldots,Y_{iT_0-p}(0))$ consists of a finite number of lagged outcomes, histories prior to $T_0-p$ have no additional predictive power for future untreated outcomes or treatment. Assumption~\ref{assumption:surrogacy_simultaneous}\ref{item:assumption_surrogacy_simultaneous_product} further requires projection errors to remain uncorrelated even after conditioning on earlier histories $H_{iT_0-1}$.
Importantly, Assumption~\ref{assumption:surrogacy_simultaneous} is weaker than a standard Markov assumption on the untreated outcome process, which restricts the conditional distribution rather than only the conditional mean.
\subsection{Potential outcomes}
I assume that the untreated outcome $Y_{it}(0)$ is generated according to
\begin{align}\label{eq:potential_outcome}
Y_{it}(0) = f (X_{it}, \alpha_i, \gamma_t ) + \epsilon_{it}, \quad \mathbb{E}_T(\epsilon_{it} \mid H_{it-1}, X_{it}, \alpha_i) = 0.
\end{align}
Here, $f(\cdot)$ is a data-generating function, $\alpha_i$ and $\gamma_t$ are unobserved individual and time fixed effects, respectively. The exogeneity condition $\mathbb{E}_T(\epsilon_{it} \mid H_{it-1}, X_{it}, \alpha_i) = 0$ corresponds to the weak exogeneity condition in structural panel model literature that allows for feedbacks~\citep{bonhomme2025back}. The data-generating process in~\eqref{eq:potential_outcome} is flexible. I impose no parametric restrictions on the function $f$, allowing for rich interactions between individual-specific and time-specific heterogeneity. In addition, I do not specify the joint distribution of $X_{it}$ and the latent factors.
The model accommodates dynamic dependence in untreated potential outcomes when $X_{it}:=(Y_{it-1}(0),\ldots,Y_{it-p}(0))$. Dynamic dependence is ubiquitous in economics. Aggregate outcomes such as local employment and fertility rates often exhibit dynamic dependence because shocks are persistent over time, while individual outcomes may also exhibit state dependence, i.e., past employment may affect present employment \citep{heckman1981heterogeneity}. Such dynamic dependence is particularly relevant in staggered-adoption settings, as it allows researchers to disentangle treatment effects from dependence on lagged outcomes.
Model~\eqref{eq:potential_outcome} incorporates a number of important data-generating processes that are commonly used in the structural and causal panel literature. Here are some examples.
\begin{example}[Linear panel with IFE]\label{example:interactive_fixed_effects}
Consider the following linear panel model with covariates $X_{it}$ and interactive fixed effects:
\begin{equation}\label{eq:interactive_fixed_effects}
\begin{aligned}
Y_{it}(0) = \beta' X_{it} + \gamma_{t}'\alpha_i + \epsilon_{it}, \quad \mathbb{E}_{T} (\epsilon_{it} \mid H_{it-1}, X_{it} , \alpha_i) = 0.
\end{aligned}
\end{equation}
Here, $\gamma_t$ is a vector of macroeconomic shocks, e.g., technology shocks and financial crises, and $\alpha_i$ captures heterogeneous responses to these shocks. For example, in labor economics, $Y_{it}$ represents the wage rate, $X_{it}$ represents experience, $\alpha_i$ is a vector of unobserved skills, and $\gamma_t$ contains the time-varying market returns to these skills. Changes in technology or labor-market conditions therefore affect workers differently according to their latent skill composition.
Model~\eqref{eq:interactive_fixed_effects} is studied in~\citet{bai2009panel} and nests many commonly used untreated-outcome specifications in causal panel methods. For example, setting $\beta=0$ and imposing additive factor structure gives $Y_{it}(0)=\alpha_i+\gamma_t+\epsilon_{it}$, which is the untreated-outcome structure in standard DiD designs. Setting $\beta=0$ while allowing unrestricted interactive fixed effects gives $Y_{it}(0)=\alpha_i'\gamma_t+\epsilon_{it}$, which is typically assumed to support recently developed synthetic control and related factor-based causal panel methods \citep{arkhangelsky2021synthetic,athey2021matrix, bai2021matrix, chernozhukov2023inference}.
\end{example}
\begin{example}[Dynamic linear panel with IFE]\label{example:dynamic_interactive_fixed_effects}
Consider the following dynamic linear panel model with $X_{it}: = Y_{it-1}(0)$ and interactive fixed effects:
\begin{align}
\label{eq:dynamic_interactive_fixed_effects}
Y_{it}(0) = \rho Y_{it-1}(0) + \alpha_i'\gamma_t + \epsilon_{it}.
\end{align}
Here, the error terms are sequentially exogenous, i.e., $\mathbb{E}_{T}(\epsilon_{it} \mid \{Y_{i\tau}\}_{\tau \leq t-1}, \alpha_i ) = 0$. Model~\eqref{eq:dynamic_interactive_fixed_effects} captures outcome persistence through the lagged outcome $Y_{it-1}(0)$ and allows individuals to respond heterogeneously to time-varying aggregate shocks through the interactive term $\alpha_i'\gamma_t$.
This model is studied in \citet{moon2017dynamic}. A special case with additive fixed effects is $Y_{it}(0)=\rho Y_{it-1}(0)+ \alpha_i+\gamma_t+\epsilon_{it}$, which is a benchmark specification in the dynamic panel literature \citep{arellano1991some}.
Dynamic dependence is crucial in many economic settings. For example, workers' current earnings may depend on their earnings histories because income shocks are persistent over time, and firms' investment decisions may depend on their past investment because of adjustment costs.
\end{example}
\begin{example}[Dynamic binary panel with IFE]\label{example:dynamic_nonlinear_interactive_fixed_effects}
Consider the following binary response panel model,
\begin{align}
\label{eq:dynamic_nonlinear}
Y_{it}(0) = \boldsymbol{1}(\beta Y_{it-1}(0) + \alpha_i'\gamma_t- u_{it} \geq 0),
\end{align}
where the exogenous errors $\{u_{it}\}$ are independent (across $i$ and $t$) draws from a distribution with cumulative distribution function $F(\cdot)$.\footnote{E.g., $F(\cdot)$ can stand for the logistic or standard normal CDF in Logit and Probit models, respectively.}
The binary panel \eqref{eq:dynamic_nonlinear} allows for heterogeneous responses of units to aggregate shocks and state dependence, e.g., women's labor-force participation decisions may depend on their employment histories, as those who have been out of the labor market for a long time may be less likely to re-enter.
In this example, when $F(\cdot)$ can stand for the logistic distribution, the data generating process $f(\cdot)$ as in~\eqref{eq:potential_outcome} takes the form
\begin{align*}
f (Y_{it-1}(0), \alpha_i, \gamma_t ) = \frac{\exp\left(\beta Y_{it-1}(0) + \alpha_i'\gamma_t\right)}{1 + \exp\left(\beta Y_{it-1}(0) + \alpha_i'\gamma_t\right)},
\end{align*}
and $\epsilon_{it}$ is a centered Bernoulli error with $\mathbb{E}_T (\epsilon_{it} \mid \{Y_{i\tau}(0)\}_{\tau\leq t-1}, \alpha_i) = 0$.
\end{example}
\begin{assumption}[Potential Outcome]\label{assumption:potential_outcome_simultaneous}
$Y_{it}(0)$ evolves according to~\eqref{eq:potential_outcome} and satisfies:
\begin{enumerate}[label=(\roman*)]
\item \label{item:potential_outcome_simultaneous_time_factor_exogeneity} \textbf{(Time factor exogeneity)} conditional on the individual latent factors $\alpha_i$ and the time-factor path $\Gamma_t$, $Y_{it}(0)$ is mean independent of future time factors, i.e.,
\begin{align*}
\mathbb{E}\left(Y_{it}(0) \mid \alpha_i, \Gamma_T\right) = \mathbb{E}\left(Y_{it}(0) \mid \alpha_i, \Gamma_t\right), \quad \forall t.
\end{align*}
\item \label{item:potential_outcome_simultaneous_weak_dependence} \textbf{(Weak dependence)} $(\ldots, Y_{iT_0-1}(0), Y_{iT_0}(0), \ldots )$ are conditionally weakly dependent, such that $\sum_{s =-\infty}^{\infty} \left|\mathrm{Cov} \left(Y_{it}(0), Y_{is}(0)\mid \alpha_{i} = \alpha, \Gamma_T = \Gamma \right)\right| <\infty$ uniformly for any supported $\alpha$, any possible time path $\Gamma$, and any $t \in \mathbb{Z}$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_time_factor_exogeneity} requires that future time factors provide no additional information for predicting $Y_{it}(0)$ when conditioning on $(\alpha_i,\Gamma_t)$. This restriction is mild and has natural interpretations. First, it reflects the principle that the future cannot affect the present. This assumption does not preclude forward-looking behavior based on expectations formed using the macroeconomic information set $\Gamma_t$.
Second, since $(\gamma_{t+1}, \gamma_{t + 2}, \ldots)$ capture aggregate conditions, they should not be affected by any single individual's behavior.
Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_weak_dependence} states that, conditional on the fixed effects $(\alpha_i, \Gamma_T)$, the serial dependence of the untreated outcome vanishes sufficiently fast such that the sum of the absolute conditional autocovariances is finite. This condition holds uniformly over any supported individual fixed effects $\alpha_i$ and any possible time path $\Gamma_T$. This conditional weak dependence requirement is also mild since persistence captured by the individual and time fixed effects, $\alpha_i$ and $\Gamma_T$, is unrestricted.
The conditional weak dependence assumption is crucial to the proposed method. I discuss sufficient conditions under which it holds in Examples~\ref{example:interactive_fixed_effects}-\ref{example:dynamic_nonlinear_interactive_fixed_effects}.
\setcounter{example}{0}
\begin{example}[Linear panel with IFE - Continued]
It is common in the literature to assume that $X_{it}$ also admits a factor structure \citep{bai2009panel}. For simplicity, set $d_X = 1$. Then,
\begin{gather*}
Y_{it}(0) = \beta X_{it} + \gamma_{Y, t}'\alpha_{Y, i} + \epsilon_{Y, it}, \quad \mathbb{E}_{T} (\epsilon_{Y, it} \mid H_{it-1}, X_{it} , \alpha_i) = 0; \\
X_{it} = \gamma_{X, t}'\alpha_{X, i} + \epsilon_{X, it}, \quad \mathbb{E}_{T} (\epsilon_{X, it} \mid H_{it-1}, \alpha_i) = 0.
\end{gather*}
Here, $\alpha_{Y,i}$ and $\alpha_{X,i}$ are (possibly overlapping) subvectors of $\alpha_i$, and $\gamma_t := (\gamma_{Y,t}',\gamma_{X,t}')'$.
Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_weak_dependence} holds if I additionally assume that the idiosyncratic terms $(\epsilon_{X,it}, \epsilon_{Y,it})$ are not persistent over time, i.e., for all $\alpha\in \mathcal{A}$, all possible time paths $\Gamma$, and any $t\in \mathbb{Z}$,
\begin{align*}
\sum_{s=-\infty}^{\infty} \left|\mathrm{Cov}\left(\epsilon_{Y,it}, \epsilon_{Y,is}\mid \alpha_i = \alpha, \Gamma_T = \Gamma\right)\right| <\infty,
\quad
\sum_{s=-\infty}^{\infty} \left|\mathrm{Cov}\left(\epsilon_{X,it}, \epsilon_{X,is}\mid \alpha_i = \alpha, \Gamma_T = \Gamma\right)\right| <\infty,
\end{align*}
which are standard in the factor-model literature. The result extends to multidimensional $X_{it}$ by imposing the same weak-dependence condition on each component of $\epsilon_{X,it}$.
\end{example}
\begin{example}[Dynamic linear panel with IFE - Continued]
For the dynamic linear model with IFE, \eqref{eq:dynamic_interactive_fixed_effects} in Example~\ref{example:dynamic_interactive_fixed_effects}, it is straightforward to verify that Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_weak_dependence} holds if I additionally assume that (i) $|\rho| < 1$, and (ii) the idiosyncratic term $(\epsilon_{it})$ are not persistent over time, i.e., for all $\alpha\in \mathcal{A}$, all possible time paths $\Gamma$, and any $t\in \mathbb{Z}$,
\begin{align*}
\sum_{s=-\infty}^{\infty} \left|\mathrm{Cov}\left(\epsilon_{it}, \epsilon_{is}\mid \alpha_i = \alpha, \Gamma_T = \Gamma\right)\right| <\infty.
\end{align*}
The condition $|\rho|<1$ is standard in the time-series literature to ensure that the serial dependence decays geometrically over time. The result extends to higher-order lags as long as all roots of the autoregressive characteristic polynomial lie outside the unit circle.
\end{example}
\begin{example}[Dynamic nonlinear panel - Continued]
For the dynamic model \eqref{eq:dynamic_nonlinear} in Example~\ref{example:dynamic_nonlinear_interactive_fixed_effects}, since $\{u_{it}\}_{t\in\mathbb{Z}}$ is independent across $t$, $\{Y_{it}(0)\}_{t\in\mathbb{Z}}$ is a time-inhomogeneous Markov process conditional on the fixed effects $(\alpha_i,\Gamma_T)$.
If I further assume that (i) the cumulative distribution function $F(\cdot)$ is continuous and its derivative satisfies $F^{(1)}(\cdot)>0$, and (ii) the supports of $\alpha_i$ and $\gamma_t$ are uniformly bounded across $i$ and $t$, then each entry of the transition matrix is bounded away from $0$ and $1$ uniformly over all periods. This ensures that Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_weak_dependence} holds. The commonly used logit and probit specifications satisfy these restrictions on $F(\cdot)$.
\end{example}
\section{Method}\label{sec:overview_method}
This section provides an overview of the proposed method, which uses pretreatment histories to recover information about latent heterogeneity, and introduces a doubly robust estimator for the dynamic ATT. I focus on simultaneous treatment settings and formally extend the analysis to staggered adoption in Section~\ref{sec:extension}.
\subsection{Recovering latent heterogeneity}
The key technical challenge is that the individual fixed effects $\alpha_i$ are unobserved by econometricians. If $\alpha_i$ were observed, the problem would fall under selection on observables, and one could directly apply existing methods to estimate the dynamic ATT. Therefore, the problem of estimating dynamic $\mathrm{ATT}$ reduces to recovering information about the latent factors $\alpha_i$ from the observed information.
To see the basic idea, first consider a pure factor model without dynamics:
\begin{align*}
Y_{it}(0) = f(\alpha_i,\gamma_t) + \epsilon_{it}, \quad \mathbb{E}(\epsilon_{it}\mid \alpha_i, \Gamma_T) = 0.
\end{align*}
This is a nonlinear factor model in which the untreated outcome depends on an individual latent factor $\alpha_i$ and a time-varying aggregate factor $\gamma_t$, and $\{\epsilon_{it}\}$ are weakly dependent across time. If two individuals $i$ and $j$ have similar latent factors, then, for each $t$, the smoothness of $f(\cdot)$ implies that $f(\alpha_i,\gamma_t) \approx f(\alpha_j,\gamma_t)$. Thus, individuals with similar latent heterogeneity should have similar untreated histories up to idiosyncratic innovations. This observation provides the intuition for the reverse direction: similarity in pretreatment histories can be informative about similarity in latent factors. This is the logic behind the pseudo-distance approach studied in \citet{feng2023optimal} and \citet{deaner2025inferring}.
However, this idea does not directly extend to models with dynamic effects. To see this, even if two individuals $i$ and $j$ have similar latent factors, the terms $f(X_{it}, \alpha_i,\gamma_t)$ and $f(X_{jt}, \alpha_j,\gamma_t)$ can still differ when their covariates differ.
Despite this complication, I show that, under certain informativeness conditions, even in the presence of dynamic effects, the latent factors can still be recovered from pretreatment histories as long as the sequence is conditionally weakly dependent as in Assumption~\ref{assumption:potential_outcome_simultaneous}\ref{item:potential_outcome_simultaneous_weak_dependence}. To motivate the general treatment, consider first the linear model in Example~\ref{example:dynamic_interactive_fixed_effects}:
\begin{align*}
Y_{it}(0) = \rho Y_{it-1}(0) + \alpha_i'\gamma_t + \epsilon_{it}, \quad \mathbb{E}(\epsilon_{it} \mid \{Y_{i\tau}(0)\}_{\tau\leq t-1}, \alpha_i, \Gamma_T) = 0.
\end{align*}
I impose $|\rho|<1$ to ensure that the sequence is conditionally weakly dependent. Iterating the model backward gives, for any lag $L >0$,
\begin{align*}
Y_{it}(0) = \rho^{L} Y_{it-L}(0) + \alpha_i' \sum_{\tau = 0}^{L-1}\rho^{\tau }\gamma_{t - \tau} + \sum_{\tau = 0}^{L-1}\rho^{\tau } \epsilon_{it-\tau}.
\end{align*}
Letting $L\to\infty$, the first term vanishes under $|\rho|<1$. Hence,
\begin{align*}
Y_{it}(0) = \underbrace{\alpha_i' \sum_{\tau = 0}^{\infty}\rho^{\tau }\gamma_{t - \tau}}_{:= g(\alpha_i, \Gamma_t)} + \underbrace{\sum_{\tau = 0}^{\infty}\rho^{\tau} \epsilon_{it-\tau}}_{\tilde{\epsilon}_{it}} = g(\alpha_i, \Gamma_t) + \tilde{\epsilon}_{it}.
\end{align*}
Thus, the dynamic linear model admits a transformed pure factor representation. The transformed innovation $\tilde{\epsilon}_{it}$ is weakly dependent over time conditional on $(\alpha_i,\Gamma_T)$, and satisfies $\mathbb{E}\left(\tilde{\epsilon}_{it} \mid \alpha_i,\Gamma_T\right) = 0$ by sequential exogeneity.
This implies that the pseudo-distance method based on pure factor models can also be applied to dynamic linear models. The intuition is that the dynamic effect of the initial condition vanishes, i.e., $\lim_{L\to\infty}\rho^L Y_{it-L}(0)=0$, and the behavior of the path is eventually driven by latent heterogeneity $\alpha_i$.
I now extend the idea from the linear dynamic model to nonlinear cases. Motivated by the idea of backward iteration, I project $Y_{it}(0)$ onto the space of functions of $(\alpha_i, \Gamma_T)$:
\begin{align}\label{eq:potential_outcome_factor_representation}
Y_{it}(0) = \underbrace{\mathbb{E}\left(Y_{it}(0) \mid \alpha_i, \Gamma_T \right)}_{:=g(\alpha_i, \Gamma_t)} + \underbrace{Y_{it}(0) - \mathbb{E}\left(Y_{it}(0) \mid \alpha_i, \Gamma_T \right)}_{\tilde{\epsilon}_{it}} = g(\alpha_i, \Gamma_t) + \tilde{\epsilon}_{it}.
\end{align}
Here, since future time factors provide no additional information about $Y_{it}(0)$ conditional on $(\alpha_i,\Gamma_t)$, the transformed pure factor representation $g(\alpha_i, \Gamma_t)$ depends only on time factors up to period $t$. In addition, the projection error $\tilde{\epsilon}_{it}$ is weakly dependent over time conditional on $(\alpha_i,\Gamma_T)$ and satisfies $\mathbb{E}\left(\tilde{\epsilon}_{it} \mid \alpha_i,\Gamma_T\right) = 0$.
Thus, I show that the dynamic nonlinear model also admits a transformed pure factor representation.
Example~\ref{example:dynamic_nonlinear_interactive_fixed_effects} provides a concrete example of this transformed pure factor representation under nonlinear models.
\setcounter{example}{2}
\begin{example}[Dynamic nonlinear panel - Continued]
Suppose that $Y_{it}(0)$ evolves as in Example~\ref{example:dynamic_nonlinear_interactive_fixed_effects} and let $F(\cdot)$ denote the cumulative distribution function of $u_{it}$. For notational simplicity, define $\rho_{it}: = F\left(\beta + \gamma_t' \alpha_i \right) - F \left( \gamma_t' \alpha_i \right)$ as the time-inhomogeneous discount factor, which satisfies $\sup_{i, t} |\rho_{it}| < 1$. It is straightforward to verify that one-step backward iteration gives
$g(\alpha_i, \Gamma_t) = F \left(\gamma_t' \alpha_i \right) + \rho_{it}g(\alpha_i, \Gamma_{t-1})$.
Then, for any $L>0$, $L$-step backward iteration gives
\begin{align*}
g(\alpha_i, \Gamma_t)
= & \sum_{\tau = 0}^{L-1} \left(\prod_{s = 0}^{\tau-1} \rho_{i, t-s} \right) F \left(\gamma_{t-\tau}' \alpha_i \right) + \left(\prod_{s = 0}^{L-1} \rho_{i, t-s}\right) g(\alpha_i, \Gamma_{t-L}),
\end{align*}
where the empty product is equal to one. Letting $L\to\infty$, the last term vanishes as $\sup_{i,t}|\rho_{it}|<1$. Hence,
\begin{align*}
g(\alpha_i, \Gamma_t) =\sum_{\tau = 0}^{\infty} \left(\prod_{s = 0}^{\tau-1} \rho_{i, t-s} \right) F \left(\gamma_{t-\tau}'\alpha_i \right).
\end{align*}
A rigorous treatment is provided in the appendix.
\end{example}
The pure factor representation of the dynamic model allows us to compare individuals through their pretreatment histories. Intuitively, when pretreatment histories are informative about latent factors, closeness between the pretreatment histories of individuals $i$ and $j$ implies similarity between their latent factors, i.e., $\|\alpha_i - \alpha_j\|\approx 0$.
A natural measure of the difference between pretreatment histories is the squared $L_2$ distance. Specifically, for any $i, j$,
\begin{align*}
\widehat{d}_{2, ij}^2 := \frac{1}{T_0} \sum_{t = 0}^{T_0-1}(Y_{it} - Y_{jt})^2.
\end{align*}
However, $\widehat{d}_{2, ij}^2$ cannot serve as a good distance for comparing histories, because it captures not only differences in the factor component of the history but also the conditional variance of idiosyncratic innovations (see the discussion in Appendix~\ref{appendix_sub:informativeness_comparison}).
As an alternative, I employ the pseudo-distance approach from the nonlinear factor literature \citep{feng2023optimal,deaner2025inferring}. Specifically, for any $i, j$, define the pseudo-distance as:
\begin{align}\label{eq:definition_pseudo_distance}
\widehat{d}_{ij} = \max_{\substack{k_1, k_2 =1, \ldots, n, \\ k_1, k_2 \neq i, j}} \left|\frac{1}{T_0} \sum_{t = 0}^{T_0-1} (Y_{k_1t} - Y_{k_2t})(Y_{it} - Y_{jt})\right|.
\end{align}
Unlike the $L_2$ distance, the pseudo-distance effectively filters the idiosyncratic innovations by comparing the history difference between $i$ and $j$ against history differences between other individuals. It therefore provides a cleaner measure of similarity in the denoised histories.
I provide sufficient conditions in the following sections to ensure that the pseudo-distance $\widehat{d}_{ij}$ can serve as a proxy for latent similarity, $\|\alpha_i-\alpha_j\|$. The key requirement is \emph{informativeness}: individuals that are close in terms of the pseudo-distance should also be close in their underlying latent factors, which means that one can find individuals with similar latent heterogeneity using the pseudo-distance. This condition is related to completeness conditions widely used in nonparametric identification problems. I formally define the informativeness condition and discuss sufficient conditions under which it holds in Sections~\ref{sec:identification} and~\ref{sec:estimation}.
\begin{remark*}\textbf{(Common time trend)}
The pseudo-distance remains well-defined even when $Y_{it}(0)$ contains a common time trend, in which case the factor representation of the data-generating process takes the form
\begin{align*}
Y_{it}(0)
= \underbrace{\lambda_t + \tilde{g}(\alpha_i,\Gamma_t)}_{=g(\alpha_i,\Gamma_t)}
+ \tilde{\epsilon}_{it}
= g(\alpha_i,\Gamma_t)+\tilde{\epsilon}_{it},
\end{align*}
where $\lambda_t$ is common across units and may be nonstationary and unbounded. The pseudo-distance proposed in this paper, unlike the pseudo-distances in~\citet{feng2023optimal}, \citet{feng2020causal}, and~\citet{deaner2025inferring}, accommodates such common trends without modification because these trends cancel in the cross-sectional differences entering construction~\eqref{eq:definition_pseudo_distance}.
Allowing for common time trends is important in many empirical applications. For example, aggregate outcomes such as local GDP often exhibit long-run trends. In addition, long outcome histories may span major aggregate structural breaks, such as World War II and the baby boom. Such breaks do not affect the analysis as long as they are captured by additive time effects.
\end{remark*}
Although $\alpha_i$ is unobserved, when the pseudo-distance can serve as a proxy for the latent distance $\|\alpha_i-\alpha_j\|$, I can identify and estimate the dynamic $\mathrm{ATT}$. For identification, the pseudo-distance allows us to identify units with identical latent heterogeneity, so that the problem reduces to identifying the $\mathrm{ATT}$ under selection on observables. For estimation, standard nonparametric methods, such as kernel regression or $k$-nearest-neighbor methods, can then be used to impute missing potential outcomes and estimate propensity scores.
\begin{remark*}\textbf{(Alternatives to long history)}
Although this paper uses long pretreatment histories to recover time-invariant heterogeneity, this is not the only way to do so. For example, when the panel is short but rich cross-sectional covariates are available, these covariates can instead be used to infer latent heterogeneity. This idea is closely related to factor-augmented regression \citep{stock2002forecasting} and \citet{feng2020causal}. This paper focuses on long pretreatment histories to maintain a setting comparable to synthetic control and matrix completion methods.
\end{remark*}
\subsection{Estimator}
When $\widehat{d}_{ij}$ can serve as a proxy for the latent distance, I construct a doubly robust estimator based on the observed covariates and the pseudo-distance.
For notational convenience, let $m_{it} := m_{t}(X_{iT_0}, \alpha_i)$ and $\pi_{i} := \pi (X_{iT_0}, \alpha_i)$. Under the selection mechanism in Assumption~\ref{assumption:selection_simultaneous} and the mean-independent Markov condition in Assumption~\ref{assumption:surrogacy_simultaneous}, an oracle estimator for $\mathrm{ATT}(t)$ takes the following doubly robust form:
\begin{align*}
\widehat{\mathrm{ATT}(t)}^{\mathrm{oracle}} := \frac{1}{n_{1}} \sum_{i = 1}^{n} \left(D_{i}Y_{it} - \frac{(1 - D_i)\pi_iY_{it} + (D_i - \pi_i)m_{it}}{1 - \pi_i }\right),
\end{align*}
where $n_1:=\sum_{i=1}^{n} \boldsymbol{1}(D_i = 1)$ is the number of treated individuals. The estimator is oracle in the sense that it uses the true propensity score $\pi_i$ and outcome regression $m_{it}$, which are unknown in practice. Let $\widehat{\pi}_i$ and $\widehat{m}_{it}$ be estimates of $\pi_i$ and $m_{it}$ respectively. The corresponding doubly robust estimator takes the form:
\begin{equation}\label{eq:ATT_DR_simultaneous}
\begin{aligned}
\widehat{\mathrm{ATT}(t)} := \frac{1}{n_{1}} \sum_{i = 1}^{n} \left(D_{i}Y_{it} - \frac{(1 - D_i)\widehat{\pi}_iY_{it} + (D_i - \widehat{\pi}_i)\widehat{m}_{it}}{1 - \widehat{\pi}_i }\right).
\end{aligned}
\end{equation}
\paragraph{Nadaraya-Watson estimator} If latent heterogeneity $\alpha_i$ were observable, researchers could apply the Nadaraya-Watson (NW) estimator to nonparametrically estimate the functions $m_t(\cdot,\cdot)$ and $\pi(\cdot,\cdot)$, and thus $m_{it}$ and $\pi_i$.
Since the pseudo-distance $\widehat{d}_{ij}$ serves as a proxy for the latent distance $\|\alpha_i-\alpha_j\|$, I therefore replace the latent distance in the NW estimator with the pseudo-distance $\widehat{d}_{ij}$ and obtain a feasible estimator. Specifically, for any $i,j=1,\ldots,n$, define the kernel weight
\begin{align*}
\widehat{K}_{h, ij} := K\left(\left(X_{jT_0} - X_{iT_0}\right)/ h \right) K\left(\widehat{d}_{ij}/ h\right).
\end{align*}
Here, $K(\cdot)$ denotes the (product) kernel function and $h\rightarrow 0$ is the bandwidth. The hat notation emphasizes that the kernel weight is constructed using the estimated pseudo-distance $\widehat{d}_{ij}$. The Nadaraya-Watson estimators based on the pseudo-distance are given by
\begin{align*}
\widehat{m}_{it} = \sum_{j =1}^{n} (1 - D_j)Y_{jt}\widehat{K}_{h, ij} \big / \sum_{j =1}^{n} (1 - D_j) \widehat{K}_{h, ij}, \quad \widehat{\pi}_i = \sum_{j =1}^{n} D_j \widehat{K}_{h, ij} \big / \sum_{j =1}^{n} \widehat{K}_{h, ij}.
\end{align*}
Plugging these estimates into \eqref{eq:ATT_DR_simultaneous} yields a feasible doubly robust estimator for $\mathrm{ATT}(t)$.
\paragraph{Double cross-fitting}
Instead of directly plugging the nuisance-function estimates into the doubly robust estimator, researchers employ cross-fitting to debias the estimator and facilitate valid inference. This approach follows the double machine learning literature \citep[e.g.,][]{chernozhukov2018double}, which uses separate samples to estimate the nuisance functions (the potential-outcome regressions and propensity scores) and evaluate the ATT, thereby reducing overfitting bias. I further employ a \emph{double} cross-fitting procedure \citep{newey2018cross, mcclean2026double} to estimate the outcome regressions and propensity scores on separate samples, which reduces the dependence between $\widehat{m}_{it}$ and $\widehat{\pi}_i$. I show that, compared with standard cross-fitting, combining double cross-fitting with undersmoothing yields a faster convergence rate and produces a root-$n$ consistent, asymptotically normal, and asymptotically unbiased ATT estimator in more general settings.
\begin{algorithm}[htbp]
\caption{Doubly robust estimator for dynamic ATT (Basic Idea)}\label{alg:simultaneous_basic}
\begin{algorithmic}
\Require Kernel $K(\cdot)$, bandwidths $h_m, h_\pi$, and confidence level $1 - \alpha$.
\Ensure Dynamic ATT estimates $\widehat{\mathrm{ATT}}(t)$, standard errors $\widehat{\mathrm{se}}_t$, and level $1-\alpha$ confidence interval.
\State
\State \textbf{Step 1: Sample-splitting}
\State Randomly partition the sample into $3$ disjoint folds $\{\mathcal{I}_1, \mathcal{I}_2, \mathcal{I}_3\}$ of equal sizes.
\State\State \textbf{Step 2: Calculate similarity}
\State For each $i, j = 1, \ldots, n$, calculate the pseudo-distance using pretreatment outcomes:
\begin{align*}
\widehat{d}_{ij} = \max_{\substack{k_1, k_2 =1, \ldots, n, \\ k_1, k_2 \neq i, j}} \left|\frac{1}{T_0} \sum_{t = 0}^{T_0-1} (Y_{k_1t} - Y_{k_2t})(Y_{it} - Y_{jt})\right|.
\end{align*}
\State \rule{0pt}{0.4em}
\State\State \textbf{Step 3: Compute expected outcomes and propensity scores}
\State The Nadaraya-Watson estimators under double cross-fitting are given by
\begin{align*}
\widehat{m}_{it} = \frac{\sum_{j\in \mathcal{I}_m(i)} (1 - D_j)Y_{jt}\widehat{K}_{h_m, ij}}{\sum_{j\in \mathcal{I}_m(i)} (1 - D_j) \widehat{K}_{h_m, ij}} , \quad \widehat{\pi}_i = \frac{\sum_{j\in \mathcal{I}_\pi(i)} D_j \widehat{K}_{h_\pi, ij} }{\sum_{j\in \mathcal{I}_\pi (i)} \widehat{K}_{h_\pi, ij}}.
\end{align*}
For any $i,j=1,\ldots,n$ and $h>0$, the kernel weight $\widehat{K}_{h, ij}$ is defined as:
\begin{align*}
\widehat{K}_{h, ij} := K\left(\left(X_{jT_0} - X_{iT_0}\right)/ h \right) K\left(\widehat{d}_{ij}/ h\right).
\end{align*}
\State\State \textbf{Step 4: Construct $\widehat{\mathrm{ATT}(t)}$, standard error $\widehat{\mathrm{se}}_t$, and confidence interval}
\State Let $n_1 :=\sum_{i=1}^{n} \boldsymbol{1}(D_i = 1)$, and compute
\begin{align*}
\widehat{\mathrm{ATT}(t)} = &\frac{1}{n_{1}} \sum_{i = 1}^{n} \left(D_{i}Y_{it} - \frac{(1 - D_i)\widehat{\pi}_iY_{it} + (D_i - \widehat{\pi}_i)\widehat{m}_{it}}{1 - \widehat{\pi}_i }\right), \\
\widehat{\mathrm{se}}_t = & \left(\frac{1}{n_1^2} \sum_{i=1}^{n} \left(D_{i}Y_{it} - \frac{(1 - D_i)\widehat{\pi}_iY_{it} + (D_i - \widehat{\pi}_i)\widehat{m}_{it}}{1 - \widehat{\pi}_i } - D_i \widehat{\mathrm{ATT}(t)} \right)^2\right)^{1/2}.
\end{align*}
The confidence interval is $\mathrm{CI}_t = [\widehat{\mathrm{ATT}}(t) \pm Z_{1 - \alpha/2} \cdot \widehat{\mathrm{se}}_t ]$.
\end{algorithmic}
\end{algorithm}
The double cross-fitting procedure is implemented as follows. I randomly partition the sample into three separate folds, $\{\mathcal{I}_1, \mathcal{I}_2, \mathcal{I}_3\}$. For each $i=1,\ldots,n$, let $\mathcal{I}(i)$ denote the subsample containing unit $i$, and let $\mathcal{I}_m(i)$ and $\mathcal{I}_\pi(i)$ denote the subsamples used to estimate the potential-outcome regression and the propensity score for unit $i$, respectively. For example, when $i\in\mathcal{I}_1$, I use $\mathcal{I}_m(i)=\mathcal{I}_2$ to estimate the potential-outcome regression and obtain $\widehat{m}_{it}$, and use $\mathcal{I}_\pi(i)=\mathcal{I}_3$ to estimate the propensity score and obtain $\widehat{\pi}_i$.\footnote{
More broadly, if $i\in\mathcal{I}_s$, then $\mathcal{I}(i)=\mathcal{I}_s$, $\mathcal{I}_m(i)=\mathcal{I}_{1+(s\bmod 3)}$, and $\mathcal{I}_\pi(i)=\mathcal{I}_{1+((s+1)\bmod 3)}$.
}
I use different bandwidths $\{h_m,h_\pi\}$ to estimate $m_{it}$ and $\pi_i$, respectively. The NW estimators under double cross-fitting are given by
\begin{align*}
\widehat{m}_{it} = \frac{\sum_{j\in \mathcal{I}_m(i)} (1 - D_j)Y_{jt}\widehat{K}_{h_m, ij}}{\sum_{j\in \mathcal{I}_m(i)} (1 - D_j) \widehat{K}_{h_m, ij}} , \quad \widehat{\pi}_i = \frac{\sum_{j\in \mathcal{I}_\pi(i)} D_j \widehat{K}_{h_\pi, ij}}{\sum_{j\in \mathcal{I}_\pi (i)} \widehat{K}_{h_\pi, ij}} .
\end{align*}
Plugging these estimates into \eqref{eq:ATT_DR_simultaneous} yields the proposed estimator for $\mathrm{ATT}(t)$. To highlight the main idea of the proposed algorithm, a simplified version that abstracts from cross-validation for selecting the bandwidths $\{h_m,h_\pi\}$ is summarized in Algorithm~\ref{alg:simultaneous_basic}.
\section{Theory for Identification}\label{sec:identification}
This section establishes identification of the dynamic ATT. If $\alpha_i$ were observed or identified, the selection mechanism would reduce to selection on observables, and the dynamic ATT could be identified directly under standard conditions. I show that, although $
\alpha$ is unobserved, the pseudo-distance is informative about latent heterogeneity, thereby enabling a matching strategy that identifies the dynamic ATT.
\paragraph{Pseudo-distance} I first discuss conditions under which the sample pseudo-distance $\widehat{d}_{ij}$ serves as a proxy for the latent distance. The argument proceeds in two steps: (i) I show that, as $n,T_0\to\infty$, $\widehat{d}_{ij}$ converges to a \emph{population} pseudo-distance $d(\alpha_i,\alpha_j)$ that depends only on the latent heterogeneity $\alpha_i$ and $\alpha_j$, and (ii) I introduce an informativeness condition under which the population pseudo-distance $d(\alpha_i,\alpha_j)$ reveals the latent distance $\|\alpha_i-\alpha_j\|$.
Establishing the population limit of $\widehat{d}_{ij}$ requires the underlying time averages to converge as $T_0\to\infty$, which in turn requires restrictions on the long-run behavior of the time fixed effects.
One way to formulate these restrictions is to condition on the path $\{\gamma_t\}_{t\in\mathbb Z}$ and impose restrictions directly on this path. This approach is consistent with the fixed-effect interpretation of $\gamma_t$ and is common in the factor model literature, but the resulting conditions are typically high-level and difficult to interpret, especially in nonlinear models.
To focus on the main idea and provide more primitive conditions, I instead embed the realized time effects in a superpopulation model. Specifically, $\Gamma_{T} $ is viewed as one realization of an underlying stochastic process, and restrictions are imposed on this process to characterize the limiting behavior of the realized path as $T_0\to\infty$\footnote{
This formulation does not change the fixed-effect interpretation in the analysis, because once the sample path is realized, the time effects are treated as fixed.
}. The superpopulation assumption serves only to provide transparent and easy-to-verify sufficient conditions for the required time-series limits. Importantly, all stochastic-process assumptions imposed below can be replaced by direct restrictions on the realized path that imply the same limiting results.
\begin{assumption}[Pseudo distance]\label{assumption:identification_d_simultaneous}
I assume that:
\begin{enumerate}[label=(\roman*)]
\item \label{item:identification_d_simultaneous_panel_iid} \textbf{(Sampling)} Conditional on $\{\gamma_t\}_{t \in \mathbb{Z}}$, the panel $\{(Y_{it}, X_{it}, D_{i}, \alpha_i)\}_{i=1, \ldots, n, t = 1, \ldots, T}$ is i.i.d. across $i$, and $\{\gamma_t\}_{t \in \mathbb{Z}}$ is stationary and ergodic over time. The panel is large with $n\rightarrow\infty$ and the number of pretreatment periods $T_0\rightarrow\infty$.
\item \label{item:identification_d_simultaneous_compact} \textbf{(Compact)} The support of $\alpha$ is compact.
\item \label{item:identification_d_simultaneous_finite_moment} \textbf{(Moment)} $\mathbb{E}\left( Y^2_{it}(0) \mid \alpha_i = \alpha, \Gamma_{t} = \Gamma , D_i = d \right) <\infty $ uniformly over all supported $\alpha$, all possible paths $\Gamma$, $d\in\{0, 1\}$, and $t = 1, \ldots, T$.
\item \label{item:identification_d_simultaneous_ULLN} \textbf{(Envelope)} The envelope of the factor structure $g(\alpha, \Gamma)$ (defined in~\eqref{eq:potential_outcome_factor_representation}) has a finite second moment, i.e., $\mathbb{E}\left(\sup_{\alpha \in \mathcal{A}} \left(g(\alpha, \Gamma) \right)^2\right) < \infty$.
\item \label{item:identification_d_simultaneous_smoothness} \textbf{(Smoothness)} The factor structure $g(\alpha, \Gamma)$ defined in~\eqref{eq:potential_outcome_factor_representation} is uniformly continuous in $\alpha$, i.e., for any $\epsilon >0$, there exists a $\delta >0$ such that
$$\sup_{\Gamma\in \mathrm{supp}(\Gamma_{T})} \sup_{\alpha, \alpha'\in \mathcal{A}, \| \alpha - \alpha' \| \leq \delta} | g(\alpha, \Gamma) - g(\alpha', \Gamma)| \leq \epsilon. $$
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_panel_iid} requires panel data with $n\rightarrow\infty$ individuals and $T_0\rightarrow\infty$ pretreatment periods, and imposes conditional independence across individuals $i$ given the realized time fixed effects $\Gamma_{T}$. The conditional independence condition is weak and does not require time-homogeneity condition \citep{chernozhukov2013average,chernozhukov2026linear}. I also assume that $\{\gamma_t\}_{t\in\mathbb{Z}}$ is stationary and ergodic over time, so that $\{\Gamma_t\}_{t\in\mathbb Z}$ is also stationary and ergodic. Consequently, the distribution of $\Gamma_t$ is invariant over $t$. These restrictions, similar to the stochastic-process restrictions imposed in \citet{feng2023optimal} and \citet{deaner2025inferring}, ensure the convergence of the time averages in the sample pseudo-distance. These stationarity and ergodicity assumptions are stronger than necessary and serve as primitive conditions for the identification argument in the main text. Appendix~\ref{appendix:informativeness} discusses alternative conditions that allow for nonstationarity.
Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_compact} is a standard compact-support condition on the latent heterogeneity $\alpha_i$.
Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_finite_moment} imposes moment conditions on untreated potential outcomes. It requires that the conditional second moment $\mathbb{E}\left(Y^2_{it}(0) \mid \alpha_i=\alpha, \Gamma_t=\Gamma, D_i=d\right)$ be uniformly bounded over the support of $\alpha$, all possible realizations of $\Gamma$, $d\in\{0,1\}$, and all time periods $t=1,\ldots,T$.
It is worth noting that, although the presence of unbounded common time trends in $Y_{it}(0)$ violates Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_finite_moment}, the analysis remains valid because the common trend cancels out in the cross-sectional differences used to construct the pseudo-distance, as discussed in Section~\ref{sec:overview_method}.
Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_ULLN} requires that the envelope of $g^2(\cdot, \cdot)$ is integrable over the time-factor path $\Gamma$. The condition is mild in that it requires only a finite second moment of the envelope and does not require $g(\alpha,\Gamma)$ to be uniformly bounded.
Assumption~\ref{assumption:identification_d_simultaneous}\ref{item:identification_d_simultaneous_smoothness} imposes that $g(\alpha,\Gamma)$ is uniformly continuous in $\alpha$ over all possible paths $\Gamma$. This assumption is necessary because it ensures that individuals with similar $\alpha$ behave similarly over a long time horizon.
\begin{proposition}[Identification of $d(\alpha_i, \alpha_j)$]\label{prop:identification_d_simultaneous}
Under Assumptions~\ref{assumption:potential_outcome_simultaneous} and~\ref{assumption:identification_d_simultaneous}, for each $i,j = 1, \ldots, n$, $\widehat d_{ij}\stackrel{p}{\longrightarrow} d(\alpha_i,\alpha_j)$, as $n,T_0\to\infty$, where $d(\alpha_i, \alpha_j)$ is given by
\begin{equation}\label{eq:definition_d}
\begin{aligned}
d(\alpha_i,\alpha_j) : = \sup_{\alpha_1, \alpha_2 \in \mathcal{A}} \left| \int (g(\alpha_{1}, \Gamma) -g(\alpha_{2}, \Gamma) )(g(\alpha_{i}, \Gamma ) -g(\alpha_{j}, \Gamma) )\mathrm{d}\mathbb{P}(\Gamma) \right|.
\end{aligned}
\end{equation}
\end{proposition}
The identification of $d(\alpha_i,\alpha_j)$ follows directly from establishing that $\widehat d_{ij}\stackrel{p}{\longrightarrow} d(\alpha_i,\alpha_j)$. This constitutes an identification-by-construction approach (see the discussion in \citet{lewbel2019identification}). The population pseudo-distance $d(\alpha_i,\alpha_j)$ depends only on $(\alpha_i,\alpha_j)$ and does not depend on the realized outcome paths or the realization of the time fixed effects. In addition, it is finite, symmetric, and uniformly continuous in $(\alpha_i,\alpha_j)$, with $d(\alpha_i,\alpha_j)=0$ whenever $\alpha_i=\alpha_j$.
The integral in~\eqref{eq:definition_d} is taken with respect to the stationary distribution of the time-factor path $\Gamma$ because, under stationarity and ergodicity, the ergodic theorem implies that time averages along the realized path $\{\Gamma_t\}_{t\in\mathbb Z}$ converge to expectations under this stationary distribution. The existence of the population pseudo-distance and the identification result, however, do not rely on stationarity or ergodicity. I provide the definition of $d(\alpha_i, \alpha_j)$ under more general conditions in Appendix~\ref{appendix_sub:definition_d_without_stationary_ergodicity}.
\begin{remark*}
Proposition~\ref{prop:identification_d_simultaneous} identifies the $n(n-1)/2$ pairwise pseudo-distances $d(\alpha_i,\alpha_j)$ in the observed sample. The proposition does not identify the latent heterogeneity $\alpha_i$ itself or the functional form of $d(\cdot,\cdot)$ over the entire latent space. Nevertheless, I show that this pairwise information is sufficient for the identification of the dynamic ATT, as it allows us to identify and match units with similar latent heterogeneity.
\end{remark*}
\paragraph{Identification of $\mathrm{ATT}(t)$}
To use the identified pseudo-distance to recover latent similarity, I impose the following informativeness condition.
\begin{assumption}[Informativeness - identification]\label{assumption:informativeness_identification_simultaneous}
Suppose the population pseudo-distance exists and for every $\varepsilon >0$, there exists a $\delta >0$, such that
$$\sup_{\alpha_1, \alpha_2\in \mathcal{A}, d(\alpha_1, \alpha_2)\leq \delta} \|\alpha_1 - \alpha_2 \| \leq \varepsilon. $$
\end{assumption}
This assumption is key to identifying dynamic ATT, which requires that, if two individuals are close in terms of $d(\alpha_i,\alpha_j)$, then their latent factors must also be close. Assumption~\ref{assumption:informativeness_identification_simultaneous}, analogous to the identification assumption in \citet[Assumption~3]{auerbach2022identification}, is mild and holds if distinct latent factors cannot generate exactly the same factor structure, i.e., for any $\alpha_1\neq\alpha_2$, $g(\alpha_1,\Gamma)$ and $g(\alpha_2,\Gamma)$ are not almost surely identical (see Lemma~\ref{lemma:sufficient_informativeness_identification_simultaneous} in Appendix~\ref{appendix_sub:informativeness_comparison} for a formal discussion). In addition, the dynamic panel models in Examples~\ref{example:dynamic_interactive_fixed_effects} and~\ref{example:dynamic_nonlinear_interactive_fixed_effects} satisfy the informativeness assumption under general conditions, and the discussion is postponed to Section~\ref{sec:estimation}. This assumption is called informativeness because, although $\alpha$ is unobserved, the pretreatment history contains sufficient information to recover similarity in latent heterogeneity through the population pseudo-distance $d(\alpha_i,\alpha_j)$.
\begin{assumption}[Identification of ATT]\label{assumption:identification_ATT_simultaneous}
I assume that:
\begin{enumerate}[label=(\roman*)]
\item \label{item:identification_ATT_simultaneous_finite_moment} \textbf{(Moment)} $\mathbb{E}_{T}\left( | Y_{it}(0)| \mid D_i = 1\right) , \mathbb{E}_{T}\left( | Y_{it}(1)| \mid D_i = 1\right)<\infty $ for each $t = T_0, \ldots, T$.
\item \label{item:identification_ATT_simultaneous_smoothness} \textbf{(Smoothness)} For each $t=T_0,\ldots,T$, the expected outcome $m_{t}(x,\alpha)$ defined in~\eqref{eq:expected_outcome_ps_simultaneous}, is uniformly continuous in $(x, \alpha)$, i.e., for any $\epsilon >0$, there exists $\delta >0$ such that for all $(x, \alpha) \in \mathrm{supp} (X_{iT_0}, \alpha_i )$,
\begin{align*}
\sup_{ \|(x', \alpha') - (x, \alpha) \| \leq \delta } \left| m_{t}(x, \alpha ) - m_{t}(x', \alpha' ) \right| \leq \epsilon.
\end{align*}
\item \label{item:identification_ATT_simultaneous_overlap} \textbf{(Overlap)} $\mathbb{P}_T(D_i=1\mid X_{iT_0},\alpha_i)\in(0,1)$ a.s.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:identification_ATT_simultaneous}\ref{item:identification_ATT_simultaneous_finite_moment} imposes finite first-moment conditions on both untreated and treated potential outcomes. Assumption~\ref{assumption:identification_ATT_simultaneous}\ref{item:identification_ATT_simultaneous_smoothness} imposes uniform continuity on expected outcome $m_t(\cdot, \cdot)$. Uniform continuity of $m_t$ is a mild and standard condition commonly imposed for nonparametric identification of conditional expectations.
Assumption~\ref{assumption:identification_ATT_simultaneous}\ref{item:identification_ATT_simultaneous_overlap} imposes the standard overlap condition, which is commonly required in causal inference to identify counterfactual outcomes.
The following theorem presents the identification result for the dynamic ATT. The key idea is that, although $\alpha$ is unobserved to econometricians, the identified population pseudo-distance reveals latent similarity. Therefore, the dynamic ATT is identified by matching units with similar latent heterogeneity and pretreatment history.
\begin{theorem}[Identification of ATT]\label{thm:identification_simultaneous}
Suppose that the population pseudo-distance $d(\alpha_i,\alpha_j)$ is identified for each $i,j=1,\ldots,n$. Then, under Assumptions~\ref{assumption:selection_simultaneous}-\ref{assumption:potential_outcome_simultaneous}, \ref{assumption:informativeness_identification_simultaneous}, and \ref{assumption:identification_ATT_simultaneous}, $\mathrm{ATT}(t)$ is identified by
\begin{align*}
\mathbb{E}_T(Y_{it} \mid D_i = 1) - \mathbb{E}_T\left( \lim_{\delta\downarrow 0}\mathbb{E}_T (Y_{jt} \mid X_{jT_0} = X_{iT_0}, d(\alpha_i, \alpha_j) \leq \delta, D_j = 0 ) \mid D_i = 1 \right).
\end{align*}
\end{theorem}
Theorem~\ref{thm:identification_simultaneous} shows that the dynamic ATT is identified using a matching strategy based on the pseudo-distance.
Identification of $\mathrm{ATT}(t)$ requires imputing the missing untreated counterfactual outcome for treated units, i.e., $\mathbb{E}_T\left(Y_{it}(0)\mid X_{iT_0},\alpha_i\right)$. I impute this counterfactual using matching. To elaborate, for each treated unit $i$, I match $i$ to untreated units $j$ with the same pretreatment covariates and a population pseudo-distance $d(\alpha_i,\alpha_j)$ no greater than $\delta$. For a fixed $\delta$, the imputed quantity is
\begin{align*}
\mathbb{E}_T(Y_{jt} \mid X_{jT_0} = X_{iT_0}, d(\alpha_i, \alpha_j) \leq \delta, D_j = 0 ).
\end{align*}
By the informativeness condition, a small pseudo-distance implies similar latent heterogeneity. Together with the smoothness and selection conditions, this allows $\delta$ to approach zero, yielding
\begin{align*}
\mathbb{E}_T\left(Y_{it}(0) \mid X_{iT_0}, \alpha_i\right) = \lim_{\delta\rightarrow 0} \mathbb{E}_T (Y_{jt} \mid X_{jT_0} = X_{iT_0}, d(\alpha_i, \alpha_j) \leq \delta, D_j = 0 ).
\end{align*}
Finally, averaging the imputed counterfactual over the distribution of treated units yields the identification formula in Theorem~\ref{thm:identification_simultaneous}.
My imputation strategy is closely related to the missing-outcome imputation approach in the panel-data causal inference literature, including DiD~\citep{borusyak2024revisiting}, synthetic control~\citep{abadie2003economic}, and factor-model imputation methods~\citep{bai2021matrix}.
I consider a more general nonlinear dynamic model, and as a result, the latent heterogeneity cannot be eliminated by differencing, as in standard DiD designs, and cannot be directly estimated using linear factor methods. The pseudo-distance identifies only pairwise latent similarity: it reveals whether two units have close latent types, but it does not identify the latent heterogeneity or its distribution. Nevertheless, Theorem~\ref{thm:identification_simultaneous} shows that even pairwise information is sufficient for identification.
The studies most closely related to this paper are \citet{feng2023optimal}, \citet{athey2025identification}, and \citet{deaner2025inferring}, which also consider imputation-based methods for nonlinear models. However, their frameworks allow selection to depend only on latent heterogeneity, whereas this paper extends the analysis to settings in which selection also depends on pretreatment covariates.
\section{Theory for Estimation and Inference}\label{sec:estimation}
This section establishes the asymptotic properties of the proposed Nadaraya-Watson estimators using the sample pseudo-distance to impute untreated outcomes and estimate propensity scores. I then provide sufficient conditions under which the doubly robust ATT estimator with double cross-fitting is $\sqrt{n}$-consistent, asymptotically unbiased, and asymptotically normal, thereby enabling valid inference for the dynamic ATT.
In this section, I focus on the case in which both the latent heterogeneity and the pretreatment covariates are continuously distributed. The analysis extends straightforwardly to settings where either $\alpha$ or $X$ has finite support.
\subsection{Asymptotic analysis of pseudo distance}
I impose two assumptions to ensure that a sufficiently small pseudo-distance $\widehat{d}_{ij}$ implies that $\alpha_i$ and $\alpha_j$ are close in the latent space. The first assumption characterizes the rate of convergence of $\widehat{d}_{ij}$ to its population counterpart $d(\alpha_i,\alpha_j)$.
\begin{assumption}[Pseudo-distance approximation]\label{assumption:estimation_consistency_d}
There exist constants $\lambda_1, \lambda_2 >0$ such that the following inequality holds wpa1:
\begin{align*}
\max_{i,j \in \{1, \ldots, n\}} |\widehat{d}_{ij} - d(\alpha_i, \alpha_j) | \leq \lambda_1 n^{-1/d_\alpha}\log (n) + \lambda_2 T_0^{-1/2}\sqrt{\log(nT_0)},
\end{align*}
\end{assumption}
Assumption~\ref{assumption:estimation_consistency_d} strengthens the pointwise convergence result used for identification by imposing a uniform convergence rate for $\widehat{d}_{ij}$. This assumption states that the estimation error between $\widehat{d}_{ij}$ and $d(\alpha_i,\alpha_j)$ is uniformly bounded by two terms in an asymptotic sense. The first term, $\lambda_1 n^{-1/d_\alpha}\sqrt{\log n}$, depends only on the dimension of the latent heterogeneity $d_\alpha$ and not on the length of the pretreatment history $p$. It captures the matching discrepancy in the latent characteristics $\alpha$, and its rate deteriorates as $d_\alpha$ increases. The second term, $\lambda_2 T_0^{-1/2}\sqrt{\log(nT_0)}$, captures the sampling error that arises because the pseudo-distance is constructed from noisy observed outcomes.
The convergence rate in Assumption~\ref{assumption:estimation_consistency_d} is consistent with the indirect matching rates in \citet[Theorem~4.1]{feng2023optimal} and \citet[Lemma~A.1]{deaner2025inferring} for pure factor models, while the setting considered here is more general in allowing for dynamic dependence. I show that this condition is not restrictive even in the dynamic panel setting and provide sufficient conditions under which it holds (see Lemma~\ref{lemma:sufficient_estimation_consistency_d} in Appendix~\ref{appendix_sub:sufficient_consistency_d_without_proof}). Lastly, for notational simplicity, let $\delta_{n, T_0} := \lambda_1 n^{-1/d_\alpha}\log (n) + \lambda_2 T_0^{-1/2}\sqrt{\log(nT_0)}$ denote a uniform upper bound on these estimation errors.
\bigskip
The second assumption is critical and requires that the population pseudo-distance $d (\alpha_i, \alpha_j)$ be
informative to reveal the latent distance $\|\alpha_i - \alpha_j\|$.
\begin{assumption}[Informativeness-estimation]\label{assumption:informativeness_estimation}
There exists a constant $\eta > 0$ such that for any $\alpha_1, \alpha_2 \in \mathcal{A}$, $\|\alpha_1 - \alpha_{2}\|\leq \eta d (\alpha_1, \alpha_{2})$.
\end{assumption}
Assumption~\ref{assumption:informativeness_estimation} strengthens the informativeness condition used for identification (Assumption~\ref{assumption:informativeness_identification_simultaneous}) by requiring the latent distance to be linearly bounded by the population pseudo-distance. Assumption~\ref{assumption:informativeness_estimation} is analogous to the conditions in \citet[Assumption~4.1]{feng2023optimal} and \citet[Assumption~4(vi)]{deaner2025inferring} and is conceptually related to completeness conditions in nonparametric identification problems.
Verifying Assumption~\ref{assumption:informativeness_estimation} is challenging, especially in nonlinear models with dynamic effects and unobserved heterogeneity. I provide sufficient conditions under which many commonly used dynamic panel models, including Examples~\ref{example:interactive_fixed_effects} and~\ref{example:dynamic_nonlinear_interactive_fixed_effects}, satisfy the informativeness assumption under general conditions.
\setcounter{example}{1}
\begingroup
\begin{example}[Linear model - Continued]
Here I focus on Example~\ref{example:dynamic_interactive_fixed_effects}, since the verification for Example~\ref{example:interactive_fixed_effects} can be obtained by setting $\rho=0$ in the argument below. Suppose that $Y_{it}(0)$ evolves according to $Y_{it}(0) = \rho Y_{it-1}(0) + \gamma_t'\alpha_i + \epsilon_{it}$, where $\{\epsilon_{it}\}_{i=1, \ldots, n, t = 1, \ldots, T}$ are sequentially mean-independent idiosyncratic errors, i.e., $\mathbb{E}_T\left(\epsilon_{it}\mid \{Y_{i\tau}(0)\}_{\tau\leq t-1}, \alpha_i\right) = 0$. Then, under certain regularity conditions on $\mathcal{A}$, Assumption~\ref{assumption:informativeness_estimation} holds if (i) $|\rho|<1$, (ii) $\{\gamma_t\}_{t\in\mathbb{Z}}$ are stationary and ergodic over time, and (iii) the largest eigenvalue of $\mathbb{E}\left(\gamma_t\gamma_t'\right)$ is finite, and its smallest eigenvalue is strictly positive. Here, the condition $|\rho|<1$ ensures that the effects of past idiosyncratic shocks decay over time. In addition, the positive smallest eigenvalue of $\mathbb{E}\left(\gamma_t\gamma_t'\right)$, which can be viewed as analogous to a strong-factor condition in the factor-model literature, ensures that each dimension of $\alpha$ makes an independent contribution to the factor component that cannot be perfectly explained by the remaining dimensions. A formal statement of this result is provided in Lemma~\ref{lemma:sufficient_informativeness_linear_IFE} in Appendix~\ref{appendix_sub:sufficient_informativeness_without_proof}.
\end{example}
\endgroup
\begin{example}[Dynamic nonlinear model - Continued]
Suppose that $Y_{it}(0)$ evolves according to a dynamic Logit model, i.e., $Y_{it}(0) = \boldsymbol{1}(\beta Y_{it-1}(0) + \alpha_i'\gamma_t - u_{it} \geq 0)$, where the idiosyncratic errors $\{u_{it}\}$ are independent draws from the standard logistic distribution. Then, under certain regularity conditions on $\mathcal{A}$, Assumption~\ref{assumption:informativeness_estimation} holds if (i) $\{\gamma_t\}_{t\in\mathbb{Z}}$ are i.i.d. over time, (ii) the supports of $\alpha$ and $\gamma_t$ are bounded, and (iii) the largest eigenvalue of $\mathbb{E}\left(\gamma_t\gamma_t'\right)$ is finite and its smallest eigenvalue is strictly positive.
Similarly, the positive smallest eigenvalue of $\mathbb{E}\left(\gamma_t\gamma_t'\right)$ ensures that each dimension of $\alpha$ makes an independent contribution to the factor component that cannot be perfectly explained by the remaining dimensions. The same result also holds for dynamic Probit models. A formal statement is provided in Lemma~\ref{lemma:sufficient_informativeness_nonlinear_IFE} in Appendix~\ref{appendix_sub:sufficient_informativeness_without_proof}.
\end{example}
\subsection{Asymptotic analysis of $\widehat{\mathrm{ATT}(t)}$}
I impose the following regularity conditions to establish the asymptotic properties of $\widehat{m}_{it}$, $\widehat{\pi}_{i}$, and $\widehat{\mathrm{ATT}(t)}$.
\begin{assumption}[Estimation]\label{assumption:estimation_simultaneous} I assume that
\begin{enumerate}[label=(\roman*)]
\item \label{item:estimation_simultaneous_panel_weak_dependent} \textbf{(Sampling)} Conditional on $\Gamma_T$, the panel $\{(Y_{it},X_{it}, D_{i}, \alpha_i)\}_{i=1, \ldots, n, t = 1, \ldots, T}$ is i.i.d. across $i$. The panel is large with $n\rightarrow\infty$ and the number of pretreatment periods $T_0\rightarrow\infty$.
\item \label{item:estimation_simultaneous_finite} \textbf{(Bounded)} $Y_{it}(0)$ is uniformly bounded for all $i=1, \ldots, n$ and all $t= 1, \ldots, T$. In addition, for each $t=T_0,\ldots,T$, $\mathbb{E}_T\left(|Y_{it}(1)|^{2+\delta}\right)<\infty$ for some $\delta >0$.
\item \label{item:estimation_simultaneous_compact} \textbf{(Compact)} The support of $\alpha$, $\mathcal{A}$, is compact. Let $\rho_1$ and $\rho_2$ denote the radius of $\mathcal{A}$ and $X_{iT_0}$'s support, respectively. There exist constants $\underline{c}_1, \overline{c}_1 >0$, such that for any $\alpha_0 \in \mathcal{A}$ and any $ r \in (0, \rho_1]$, $\underline{c}_1 r^{d_{\alpha}} \leq \mathbb{P}\left(\|\alpha - \alpha_0\| \leq r \right) \leq \overline{c}_1 r^{d_{\alpha}}$. In addition, there exist constants $0<\underline{c}_2\leq\overline {c}_2<\infty$ such that, for any $(x,\alpha)\in\mathrm{supp}(X_{iT_0},\alpha_i) $ and any $ r \in (0, \max\{\rho_1, \rho_2)\}]$, $\underline{c}_2 r^{d_X + d_{\alpha}} \leq \mathbb{P}_T \left(\|X_{iT_0}-x \| \leq r, \|\alpha_i-\alpha\|\leq r \right)\leq \overline c_2r^{d_X + d_{\alpha}}$.
\item \label{item:estimation_simultaneous_L_continuous} \textbf{(Lipschitz continuity)} For each $t = T_0, \ldots, T$, the functions $m_{t}(x, \alpha)$ and $\pi (x, \alpha)$ are Lipschitz continuous in $(x, \alpha) \in \mathrm{supp} (X_{iT_0}, \alpha_i)$. \item \label{item:estimation_simultaneous_overlap} \textbf{(Common overlap)} There exist constants $\underline{p}, \overline{p}\in (0, 1)$ such that $\underline{p} \leq \pi( X_{iT_0}, \alpha_i)\leq \overline{p}$ a.s.
\item \label{item:estimation_simultaneous_kernel} \textbf{(Kernel function)} The kernel $K: \mathbb{R}\mapsto \mathbb{R}_{+}$ is bounded by $\overline{K} >0$ and supported on $[-1, 1]$. In addition, $K(0)>0$ and $K(\cdot)$ is Lipschitz continuous with constant $L_K>0$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_panel_weak_dependent} imposes conditional independence across units given the time fixed effects $\Gamma_T$.
Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_finite} requires that the untreated potential outcome $Y_{it}(0)$ be uniformly bounded over individuals and time periods, while the treated potential outcome $Y_{it}(1)$ satisfies a finite $(2+\delta)$-th moment condition, which is standard for applying the central limit theorem. Since pretreatment outcomes enter as regressors in the nonparametric estimation, the boundedness of $Y_{it}(0)$ is imposed to avoid technical complications. Notably, the boundedness condition on $Y_{it}(0)$ is stronger than necessary. For example, although an unbounded common time trend in untreated potential outcomes violates this condition, all analysis remains valid because those trends are canceled out by cross-sectional differencing (see the discussion in Section~\ref{sec:overview_method}).
Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_compact} first imposes compactness on the support of $\alpha$ and requires that the distribution of $\alpha$ be neither locally too sparse nor too concentrated. That is, for any supported $\alpha$ and any small radius $r$, the probability that $\alpha_i$ falls within a ball of radius $r$ centered at $\alpha$ is of the same order as the volume of that ball. In addition, this assumption imposes the same restriction on the joint distribution of $(X_{iT_0}, \alpha)$ and requires that the joint distribution of $(X_{iT_0}, \alpha_i)$ be neither locally too sparse nor too concentrated. Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_compact} is common in establishing the uniform asymptotic properties of nonparametric kernel estimators and is automatically satisfied if (i) the support of $(X_{iT_0}, \alpha)$ is compact and convex, and (ii) the density of $\alpha$ and the joint density of $(X_{iT_0}, \alpha)$ are uniformly bounded above and bounded away from zero.
Lastly, although this assumption is primarily applicable to continuously distributed $\alpha$ and $X$, the density condition can be replaced by the corresponding condition on point-mass probabilities to accommodate the case in which either $\alpha$ or $X$ has finite support.
Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_L_continuous} requires that both the expected outcomes and propensity scores be Lipschitz continuous functions of pretreatment histories and individual heterogeneity. This condition is mild and does not require differentiability.
Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_overlap} is a commonly adopted overlap condition in the causal inference literature. Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_kernel} imposes standard regularity conditions on the kernel function. These conditions are satisfied by many commonly used kernels, including the Epanechnikov kernel.
\begin{proposition}[Convergence rates for counterfactual imputation]\label{prop:estimation_entry_simultaneous}
Under Assumptions~\ref{assumption:selection_simultaneous}-\ref{assumption:potential_outcome_simultaneous} and Assumptions~\ref{assumption:estimation_consistency_d}-\ref{assumption:estimation_simultaneous}, for any $h\in\{h_{\pi},h_m\}$ such that $h\rightarrow 0$, $nh^{d_\alpha+d_{X}}/\log n\rightarrow\infty$, and $\delta_{n,T_0}/h\rightarrow 0$, the following result holds for each $t = T_0, \ldots, T$,
\begin{align*}
\max_{i = 1, \ldots, n}\left|\widehat{m}_{it} - m_{it} \right| = O_P\left(h_m + \left(nh_m^{d_\alpha+d_{X}}\right)^{-1/2}\sqrt{\log (n)} \right).
\end{align*}
In addition,
\begin{align*}
\max_{i = 1, \ldots, n}\left|\widehat{\pi}_{i} - \pi_{i} \right| = O_P\left(h_{\pi} + \left(nh_{\pi}^{d_\alpha+d_{X}}\right)^{-1/2}\sqrt{\log (n)} \right).
\end{align*}
\end{proposition}
Proposition~\ref{prop:estimation_entry_simultaneous} establishes the uniform convergence rates for counterfactual outcome estimators and propensity score estimators over all individuals. These convergence rates imply that, when the estimation error of the sample pseudo-distance is asymptotically negligible relative to the bandwidth $h$, the convergence rate of the estimator using $\widehat{d}_{ij}$ is identical (up to a logarithmic term) to that of the Nadaraya-Watson estimator if $\alpha$ were observed.
In standard practice, researchers set $h_m, h_\pi \asymp n^{-1/(d_\alpha+d_{X}+2)}$ to balance the bias and variance terms. The uniform convergence rates in this case are of order $n^{-1/(d_\alpha+d_{X}+2)}\log(n)$ and achieve Stone's optimal rate under Lipschitz continuity (up to a logarithmic term). Then, one typically employs standard single cross-fitting and plugs these nuisance-function estimators into the doubly robust estimator of the ATT. However, I show that the bandwidth choice that is rate-optimal for imputing counterfactual outcomes and estimating propensity scores is not necessarily optimal for estimating the ATT. By combining double cross-fitting with a different choice of bandwidths, I obtain a faster convergence rate for the ATT estimator and establish root-$n$ consistency, asymptotic normality, and asymptotic unbiasedness under a broader range of settings. The following Theorem formalizes this idea.
\begin{theorem}[Asymptotic theory for ATT estimator]\label{thm:estimation_ATT_simultaneous}
Under conditions in Proposition~\ref{prop:estimation_entry_simultaneous}, for any $h\in\{h_{\pi},h_m\}$ such that $h\rightarrow 0$, $nh^{d_\alpha+d_{X}}/\log n\rightarrow\infty$, and $\delta_{n,T_0}/h\rightarrow 0$, the following result holds for each $t=T_0, \ldots, T$,
\begin{align*}
\widehat{\mathrm{ATT}(t)} - \mathrm{ATT}( t) = \frac{1}{\sqrt{n}} \mathcal{N}(0, V_t) + O_P\left(h_m \left(h_{\pi} + (nh_{\pi}^{d_{\alpha} + d_{X}})^{-1/2}\sqrt{\log(n)}\right)\right) + o_P(n^{-1/2}).
\end{align*}
Here, $\mathcal{N}(\cdot,\cdot)$ denotes the normal distribution, and $V_t$ is the asymptotic variance of the oracle estimator.
\end{theorem}
Theorem~\ref{thm:estimation_ATT_simultaneous} is the main result of this paper. Ignoring logarithmic terms for simplicity, the theorem shows that the difference between the oracle estimator and the proposed doubly robust estimator using double cross-fitting is $O_P\left(h_m\left(h_{\pi}+(nh_{\pi}^{d_{\alpha}+d_{X}})^{-1/2}\right)\right)$.
Therefore, when this remainder term is $o_P(1/\sqrt n)$, the proposed estimator is asymptotically equivalent to the oracle estimator and achieves root-$n$ asymptotic normality and asymptotic unbiasedness.
Double cross-fitting achieves a faster convergence rate than standard single cross-fitting. Under standard single cross-fitting, the remainder term is the product of nuisance estimation errors, and is of order $O_P\left(\left(h_m + (nh_{m}^{d_{\alpha}+d_{X}})^{-1/2}\right)\left(h_{\pi}+(nh_{\pi}^{d_{\alpha}+d_{X}})^{-1/2}\right)\right)$. Double cross-fitting eliminates the variance term $(nh_m^{d_{\alpha}+d_{X}})^{-1/2}$ from the remainder, which enables researchers to adopt an \emph{undersmoothing} choice of $h_m$ rather than the bandwidth that is optimal for estimating $m_{it}$ and achieves a faster convergence rate for the remainder term.
The intuition behind this improvement is that double cross-fitting reduces the dependence between the outcome regression and propensity score estimators. Under single cross-fitting, $\widehat{m}_{it}$ and $\widehat{\pi}_{i}$ are constructed from the same sample and hence have correlated estimation errors. As a result, the product error is controlled by the Cauchy-Schwarz inequality, which depends on the convergence rates of both $\widehat{m}_{it}$ and $\widehat{\pi}_{i}$ in the $L_2$ norm. Double cross-fitting removes this dependence by estimating the two nuisance functions using separate samples. Therefore, the remainder term depends on the average bias of $\widehat{m}_{it}$ (which is of order $h_m$) rather than the $L_2$ norm of the estimation errors, allowing for undersmoothing the outcome regression.
\begin{corollary}[Inference]\label{corollary:estimation_ATT_simultaneous}
Under conditions in Proposition~\ref{prop:estimation_entry_simultaneous}, set $h_\pi \asymp n^{-1/(d_{\alpha} + d_{X} + 2)}$ and $h_m \asymp n^{-1/(d_{\alpha} + d_{X})}(\log (n))^{2/(d_{\alpha} + d_{X})} $. In addition, suppose $\delta_{n,T_0}/h_m \rightarrow 0$. Then, for each $t=T_0, \ldots, T$ and $d_{\alpha}+ d_{X} \leq 3$,
\begin{align*}
\sqrt n\left(\widehat{\mathrm{ATT}(t)}-\mathrm{ATT}(t)\right) \stackrel{d}{\longrightarrow} \mathcal{N}(0, V_t).
\end{align*}
Here, $V_t$ is the asymptotic variance of the oracle estimator.
\end{corollary}
Corollary~\ref{corollary:estimation_ATT_simultaneous} shows that, when $d_{\alpha}+d_{X}\leq 3$, combining double cross-fitting with undersmoothing yields a root-$n$ asymptotically normal and asymptotically unbiased ATT estimator, enabling valid inference on the dynamic ATT. I choose $h_\pi \asymp n^{-1/(d_{\alpha} + d_{X} + 2)}$, which is optimal for the MSE of the propensity score estimator, while aggressively undersmoothing $h_m$ because the outcome regression contributes to the remainder through its bias. Corollary~\ref{corollary:estimation_ATT_simultaneous} expands the range of settings in which root-$n$ inference can be obtained relative to the standard single cross-fitting, which establishes root-$n$ inference for the ATT estimator only when $d_{\alpha} + d_{X}=1$ (see \citet{deaner2025inferring}).
\begin{remark*}\textbf{(Other choices of bandwidths)}
Although Corollary~\ref{corollary:estimation_ATT_simultaneous} requires $h_m$ to be undersmoothed as aggressively as possible to achieve the fastest convergence rate of the remainder term, the bandwidth choice is not unique. It is straightforward to verify that there exists a range of bandwidth choices that make the remainder term $o_P(1/\sqrt{n})$ and therefore yield root-$n$ valid inference. I adopt this most aggressive undersmoothing choice because it leads to a concise theoretical characterization of the remainder term while also providing practical guidance for selecting bandwidths, which will be discussed later.
\end{remark*}
\begin{remark*}\textbf{(Discrete Covariates)}
Proposition~\ref{prop:estimation_entry_simultaneous}, Theorem~\ref{thm:estimation_ATT_simultaneous}, and Corollary~\ref{corollary:estimation_ATT_simultaneous} are primarily applicable to continuously distributed covariates. When $X_{it}$ takes discrete values, for example, in the binary response model as in Example~\ref{example:dynamic_nonlinear_interactive_fixed_effects}, the proposed estimator and the associated analysis remain valid, and the only difference is that $d_X$ does not enter the convergence rate asymptotically. Specifically, one only needs to replace all terms of the form $nh^{d_\alpha+d_X}$ with $nh^{d_\alpha}$ in the previous discussion to accommodate the discrete case. In addition, combining double cross-fitting with undersmoothing yields a root-$n$ asymptotically normal and asymptotically unbiased ATT estimator when $d_{\alpha}\leq 3$.
\end{remark*}
\begin{remark*}\textbf{(Asymptotics of $T_0$)}
Although Corollary~\ref{corollary:estimation_ATT_simultaneous} does not impose an explicit restriction on $T_0$, $T_0$ needs to grow sufficiently fast to achieve root-$n$ valid inference, as it affects the accuracy of $\widehat{d}_{ij}$ through $\delta_{n,T_0}$. A sufficient condition is that $T_0 \gg n^{2/(d_{\alpha}+d_{X})}$.
\end{remark*}
\begin{remark*}\textbf{(Functional form)}
The proposed approach is nonparametric and does not impose functional form restrictions on the potential outcome or propensity score. If researchers are willing to impose parametric restrictions, for example, assuming that the potential outcome follows a linear dynamic model with interactive fixed effects or that the propensity score follows a logit model with a linear index, the dimensionality restrictions can be substantially relaxed. In addition, how to incorporate high-dimensional covariates into the outcome regression (see, e.g., \citet{viviano2026dynamic}), while controlling for individual fixed effects remains a promising direction for future research.
\end{remark*}
\subsection{Bandwidth selection in practice}
The theoretical bandwidth choice in Corollary~\ref{corollary:estimation_ATT_simultaneous} provides practical guidance for selecting $h_m$ and $h_\pi$ using data-driven methods. Since $h_\pi$ is set to minimize MSE of the propensity score estimator, I propose to use standard leave-one-out cross-validation. To elaborate, given $h >0$, for each $i =1, \ldots, n$, compute the within-sample leave-one-out estimator,
\begin{equation*}
\widehat{\pi}_{i}^{(-i)}(h) = \sum_{i' \in \mathcal{I}(i)\setminus \{i\} } D_{i'} \widehat{K}_{h, ii'} \big / \sum_{i' \in \mathcal{I}(i)\setminus \{i\} } \widehat{K}_{h, ii'},
\end{equation*}
and select the optimal bandwidth $h_{\pi}^*$ as $h_{\pi}^* = \operatorname*{argmin}_{h > 0 } \frac{1}{n}\sum_{i=1}^{n}(D_i - \widehat{\pi}_{i}^{(-i)}(h) )^2$.
Since Corollary~\ref{corollary:estimation_ATT_simultaneous} requires $h_m\asymp n^{-1/(d_{\alpha}+d_{X})}(\log(n))^{2/(d_{X} + d_{\alpha})} $ to be undersmoothed as aggressively as possible, this immediately means that the neighbors for each unit is of order $\sim (\log (n))^2$. Given any bandwidth $h>0$, for each $i$, I define the effective sample size as
\begin{gather*}
n_{i, \mathrm{eff}}(h) = \big(\sum_{i' \in \mathcal{I}(i)\setminus \{i\} } (1 - D_{i'}) \widehat{K}_{h, ii'}\big)^2 \big /\sum_{i' \in \mathcal{I}(i)\setminus \{i\} } (1 - D_{i'}) \widehat{K}^2_{h, ii'},
\end{gather*}
Following the theoretical guidance, a proper choice of $h_m$ should be as small as possible while maintaining a sufficient number of effective matches of order $\sim (\log (n))^2$. Therefore, in practice, I choose predetermined constants $\kappa>0$ and $\underline q\in(0,1]$ and select the smallest bandwidth that ensures at least a fraction $\underline q$ of observations have an effective sample size larger than $(\kappa \log(n))^2$:
\begin{align*}
h_m^* = \min \left\{h >0: \frac{1}{n}\sum_{i=1}^{n}\boldsymbol{1}\left(n_{i,\mathrm{eff}}(h)>(\kappa \log(n))^2\right)\geq \underline q\right\}.
\end{align*}
Numerical simulations in Section~\ref{sec:simulation} show that $(\kappa=0.2, \underline{q} = 0.8 )$ performs well under various settings, and I recommend this value for empirical applications.
\section{Extension: Staggered Adoption}\label{sec:extension}
Section~\ref{sec:extension} extends the analysis to staggered adoptions, which are ubiquitous in empirical research. I adopt the idea of sequential unconfoundedness from \citet{robins2000marginal} and combine it with the pseudo-distance approach to identify and estimate dynamic ATT. The identification, estimation, and inference results developed for simultaneous treatment continue to hold under this extension with appropriate modifications.
\subsection{Setup and dynamic selection}
,
Consider the panel data $\{(Y_{it}, X_{it}, D_{it}, \alpha_i)\}$ on $n$ units, indexed by $i=1, \ldots, n$, over $T$ periods, indexed by $t = 1, \ldots, T$.
The outcome variable is $Y_{it}\in\mathbb{R}$. The pretreatment covariates are $X_{it} \in \mathbb{R}^{d_X}$ and can include lagged outcomes, for example, $X_{it} = Y_{it-1}$. The treatment status is denoted by $D_i\in\{0,1\}$. I assume that treatment is \emph{absorbing}, i.e., $D_{it+1}\geq D_{it}$ for any $t$, in the sense that once a unit becomes treated, it remains treated in all subsequent periods. Let $G_i$ denote the first period in which unit $i$ receives the treatment, i.e., $D_{it} = \boldsymbol{1}\{t \ge G_i\}$. For those individuals never treated in the observed time period, we set $G_i = \infty$. $\alpha_i$ is the unobservable individual fixed effect.
To allow for dynamic treatment effects, potential outcomes are indexed by the entire treatment sequence $(d_1,...,d_T) \in \{0,1\}^T$, $Y_{it}(d_1,...,d_T)$. Since treatment is an absorbing state, the potential outcomes can be indexed by the first treatment period $G_i$ only. Here, we define $Y_{it}(g): = Y_{it}(\boldsymbol{0}_{g-1}, \boldsymbol{1}_{T-g+1})$, and $Y_{it}(\infty): = Y_{it}(\boldsymbol{0}_T)$. I follow the standard non-anticipation assumption such that
\begin{align*}
Y_{it} = \left\{
\begin{array}{lr}
Y_{it}(\infty), & t < G_i\\
Y_{it}(G_i), & t\geq G_i
\end{array}
\right.
\end{align*}
Similarly, I use $\mathbb{E}_T(\cdot) := \mathbb{E}\left(\cdot \mid \Gamma_{T}\right)$ and $\mathbb{P}_T(\cdot) := \mathbb{P}\left(\cdot \mid \Gamma_{T}\right)$ to denote expectation and probability conditional on the realized time fixed effects $\Gamma_{T}$, respectively.
For each $t \geq g \geq T_0$, the estimand of interest is the dynamic ATT:
\begin{align*}
\mathrm{ATT}(g, t) := \mathbb{E}_T \left( Y_{it}(g) - Y_{it}(\infty) \mid G_i = g \right).
\end{align*}
It denotes the ATT at time $t$ for units that first receive the treatment in period $g$. As discussed in Section~\ref{sec:selection}, although the definition of the ATT here appears to differ from the standard definition, $\mathbb{E}(Y_{it}(g)-Y_{it}(\infty)\mid G_i=g)$, the two definitions are equivalent, and the difference is purely notational.
\begin{assumption}[Selection]\label{assumption:selection_extension} Assume that for each $T_0\leq g\leq t\leq T$,
\begin{align*}
Y_{it}(\infty) \perp \{G_{i} = g \} \mid X_{ig}, \alpha_i, \Gamma_{T}, G_i > g-1.
\end{align*}
\end{assumption}
The assumption extends Assumption~\ref{assumption:selection_extension} to staggered adoption, and is closely related to the sequential unconfoundedness assumption in the dynamic treatment effect literature. Assumption~\ref{assumption:selection_extension} requires that, for individuals who have not yet received the treatment before period $g$, treatment adoption at period $g$ is independent of future untreated potential outcomes, conditional on pretreatment covariates $X_{ig}$, individual fixed effects, and time factors. The conditional orthogonality is imposed only on untreated potential outcomes because the estimand of interest is the ATT, and only untreated outcomes need to be imputed.
Similar to~\eqref{eq:expected_outcome_ps_simultaneous}, for each $T_0 \leq t'\leq t\leq T$ and $\alpha \in \mathcal{A}$, define
\begin{equation}\label{eq:expected_outcome_ps_extension}
\begin{aligned}
m_{t \mid t'}(x, \alpha ) := &\mathbb{E}_T(Y_{it}(\infty) \mid X_{it'} = x, \alpha_i = \alpha, G_i > t'-1), \\
\pi_{t'}(x, \alpha ) := &\mathbb{P}(G_i=t' \mid X_{it'} = x, \alpha_i=\alpha, G_i>t'-1).
\end{aligned}
\end{equation}
The function $m_{t\mid t'}(x,\alpha)$ is the conditional expectation of the untreated outcome at period $t$ for units that remain untreated at time $t'-1$, given their pretreatment information summarized by $X_{it'} = x$ and latent heterogeneity $\alpha_i=\alpha$.
The function $\pi_{t'}(x,\alpha)$ is the propensity score of treatment adoption at time $t'$ for units that have not been treated by time $t'-1$, given their observed covariates $X_{it'}$, and individual fixed effects $\alpha_i=\alpha$. I write
\begin{align*}
Y_{it}(\infty) = & m_{t \mid t'}(X_{it'} , \alpha_i ) + u_{i, t \mid t'}, \quad &\forall\text{ } T_0 \leq t'\leq t\leq T, \text{ and } G_i >t'-1, \\
\quad D_{it'} = & \pi_{t'} (X_{it'} , \alpha_i ) + e_{it'}, \quad &\forall\text{ } T_0 \leq t' \leq T, \text{ and } G_i >t' -1.
\end{align*}
Here, $u_{i, t\mid t'}$ and $e_{it'}$ are projection errors satisfying $\mathbb{E}_T(u_{i, t\mid t'} \mid X_{it'}, \alpha_i, G_i >t' -1 ) = \mathbb{E}_T(e_{it'} \mid X_{it'}, \alpha_i, G_i >t'-1 ) = 0$. It is straightforward to check that $\mathbb{E}_T(u_{i, t\mid t'}e_{it'}\mid X_{it'},\alpha_i, G_i>t'-1)=0$ by the selection assumption. The following assumption imposes a mean-independent Markov restriction.
\begin{assumption}[Mean-independent Markov]\label{assumption:surrogacy_extension}
Let $H_{it'} := \{(Y_{i\tau}, X_{i\tau }, D_{i\tau})\}_{\tau = -\infty}^{t'}$ denote the history up to $t'$. For each $T_0\leq t' \leq t\leq T$,
\begin{enumerate}[label=(\roman*)]
\item \label{item:assumption_surrogacy_extension_outcome}$\mathbb{E}_T(u_{i, t\mid t'} \mid H_{it'-1}, X_{it'}, \alpha_i, G_i >t' -1 ) = \mathbb{E}_T(u_{i, t\mid t'} \mid X_{it'}, \alpha_i, G_i >t' -1 ) = 0$;
\item \label{item:assumption_surrogacy_extension_propensity_score} $\mathbb{E}_T(e_{it'} \mid H_{it'-1}, X_{it'}, \alpha_i, G_i >t'-1 ) = \mathbb{E}_T(e_{it'} \mid X_{it'}, \alpha_i, G_i >t'-1 ) = 0$;
\item \label{item:assumption_surrogacy_extension_product} $\mathbb{E}_T(u_{i, t\mid t'}e_{it'}\mid H_{it'-1}, X_{it'},\alpha_i, G_i>t'-1) = \mathbb{E}_T(u_{i, t\mid t'}e_{it'}\mid X_{it'},\alpha_i, G_i>t'-1) = 0$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:surrogacy_extension} extends Assumption~\ref{assumption:surrogacy_simultaneous} to staggered adoption. In addition, this assumption is weaker than a standard Markov assumption, which restricts the conditional distribution rather than only the moment condition.
\begin{assumption}[Potential Outcome]\label{assumption:potential_outcome_extension}
$Y_{it}(\infty)$ evolves according to~\eqref{eq:potential_outcome} and satisfies:
\begin{enumerate}[label=(\roman*)]
\item \label{item:potential_outcome_extension_time_factor_exogeneity} \textbf{(Time factor exogeneity)} conditional on the individual latent factors $\alpha_i$ and the time-factor path $\Gamma_t$, $Y_{it}(\infty)$ is mean independent of future time factors, i.e.,
\begin{align*}
\mathbb{E}\left(Y_{it}(\infty) \mid \alpha_i, \Gamma_T\right) = \mathbb{E}\left(Y_{it}(\infty) \mid \alpha_i, \Gamma_t\right), \quad \forall t.
\end{align*}
\item \label{item:potential_outcome_extension_weak_dependence} \textbf{(Weak dependence)} $(\ldots, Y_{it-1}(\infty), Y_{it }(\infty), \ldots )$ are conditionally weakly dependent, such that $\sum_{s =-\infty}^{\infty} \left|\mathrm{Cov} \left(Y_{it}(\infty), Y_{is}(\infty)\mid \alpha_{i} = \alpha, \Gamma_{T} = \Gamma \right)\right| <\infty$ uniformly for any supported $\alpha$, any possible time path $\Gamma$, and any $t \in \mathbb{Z}$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:potential_outcome_extension} is identical to Assumption~\ref{assumption:potential_outcome_simultaneous}, except that $Y_{it}(0)$ is replaced by $Y_{it}(\infty)$ to denote the untreated potential outcomes under staggered adoption.
\subsection{Identification}
The identification of the pseudo-distance and its informativeness under staggered adoption is identical to that under simultaneous treatment (see Assumptions~\ref{assumption:identification_d_simultaneous} and \ref{assumption:informativeness_identification_simultaneous} and Proposition~\ref{prop:identification_d_simultaneous}), except that $Y_{it}(0)$ needs to be replaced by $Y_{it}(\infty)$ to accommodate staggered adoption. Therefore, I omit the discussion of pseudo-distance identification and informativeness and focus directly on the identification of the dynamic ATT.
\begin{assumption}[Identification of ATT]\label{assumption:identification_ATT_extension}
I assume that:
\begin{enumerate}[label=(\roman*)]
\item \label{item:identification_ATT_extension_finite_moment} \textbf{(Moment)} For each $t\geq g \geq T_0$, $\mathbb{E}_{T}\left( | Y_{it}(\infty)| \mid G_i = g \right), \mathbb{E}_{T}\left( | Y_{it}(g)| \mid G_i = g \right)<\infty $.
\item \label{item:identification_ATT_extension_smoothness} \textbf{(Smoothness)} For each $T_0\leq t'\leq t\leq T$, $m_{t\mid t'}(x,\alpha)$ is uniformly continuous in $(x, \alpha)$, i.e., for any $\epsilon >0$, there exists $\delta >0$ such that for all $(x, \alpha) \in \mathrm{supp} (X_{it'}, \alpha_i )$
\begin{align*}
\sup_{ \|(x', \alpha') - (x, \alpha) \| \leq \delta } \left| m_{t\mid t'}(x, \alpha ) - m_{t\mid t'}(x', \alpha' ) \right| \leq \epsilon.
\end{align*}
\item \label{item:identification_ATT_extension_overlap} \textbf{(Overlap)} $\mathbb{P}_T(G_i = t\mid X_{it},\alpha_i, G_i >t-1)\in(0,1)$ a.s.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:identification_ATT_extension}\ref{item:identification_ATT_extension_finite_moment} is standard and imposes finite first-moment conditions on both untreated and treated potential outcomes. Assumption~\ref{assumption:identification_ATT_extension}\ref{item:identification_ATT_extension_smoothness} extends Assumption~\ref{assumption:identification_ATT_simultaneous}\ref{item:identification_ATT_simultaneous_smoothness} by imposing uniform continuity on the expected untreated outcomes $m_{t\mid t'}(\cdot, \cdot)$. Uniform continuity of $m_{t\mid t'}(\cdot, \cdot)$ is a mild and standard condition commonly imposed for nonparametric identification of conditional expectations.
Assumption~\ref{assumption:identification_ATT_extension}\ref{item:identification_ATT_extension_overlap} strengthens Assumption~\ref{assumption:identification_ATT_simultaneous}\ref{item:identification_ATT_simultaneous_overlap} and is commonly required in causal inference to identify counterfactual outcomes.
\begin{theorem}[Identification]\label{thm:identification_extension}
Suppose that the population pseudo-distance $d(\alpha_i,\alpha_j)$ is identified for each $i,j=1,\ldots,n$. Then, under Assumptions~\ref{assumption:selection_extension}-\ref{assumption:identification_ATT_extension} and Assumption~\ref{assumption:informativeness_identification_simultaneous}, the following identification results hold: (i) For each $t\geq T_0$ and each unit $i$ such that $G_i>t-1$, the contemporaneous counterfactual is identified by
\begin{align*}
m_{t\mid t}(X_{it}, \alpha_i)= \lim_{\delta\downarrow 0}\mathbb{E}_{T}(Y_{jt} \mid X_{jt} = X_{it}, d(\alpha_j, \alpha_i)\leq \delta, G_j >t);
\end{align*}
(ii) for any $t' = t-1, \ldots, g$ and each unit $i$ such that $G_i>t' -1 $, the dynamic counterfactual expectation is identified recursively backward by
\begin{align*}
m_{t\mid t'}(X_{it'}, \alpha_i) = \lim_{\delta \downarrow 0} \mathbb{E}_{T}\left( m_{t\mid t'+1 } (X_{jt' + 1}, \alpha_j)\mid X_{jt'} = X_{it'}, d(\alpha_j , \alpha_i)\leq \delta, G_j >t' \right);
\end{align*}
and (iii) for any $t\geq g\geq T_0$, the dynamic ATT is identified by
\begin{align*}
\mathrm{ATT}(g, t) = \mathbb{E}_{T}\left( Y_{it}\mid G_i = g\right) - \mathbb{E}_{T}\left( m_{t\mid g}(X_{ig}, \alpha_i) \mid G_i = g\right).
\end{align*}
\end{theorem}
Theorem~\ref{thm:identification_extension} extends the identification result to staggered-adoption settings and establishes the identification of $\mathrm{ATT}(g,t)$. It uses a pseudo-distance-based matching strategy with a backward recursive procedure to accommodate dynamic selection.
The theorem proceeds in three steps. First, it identifies the contemporaneous counterfactual expectation $m_{t\mid t}(\cdot,\cdot)$. Second, it identifies dynamic counterfactual expectations $m_{t\mid t'}(\cdot,\cdot)$ through a backward recursive procedure over $t'=t,t-1,\ldots,g$. Finally, averaging the imputed counterfactual $m_{t\mid g}(\cdot,\cdot)$ over the distribution of treated units with $G_i=g$ identifies $\mathrm{ATT}(g,t)$.
\subsection{Estimator}
For notational simplicity, for each $t \geq t' \geq T_0$ and each $i$ who has not been treated at time $t'-1$, define its expected (counterfactual) potential outcome as $m_{i, t\mid t'} := m_{t \mid t'}(X_{it'}, \alpha_i )$. In addition, for each $t = T_0, \ldots, T$ and each $i$ such that $G_i>t-1$, define its propensity score as $\pi_{it} = \pi_t(X_{it}, \alpha_i )$.
An oracle estimator for $\mathrm{ATT}(g,t)$ under staggered adoption takes the following doubly robust form:
\begin{equation}\label{eq:ATT_DR_extension}
\begin{aligned}
\widehat{ATT(g, t)}^{\mathrm{oracle}} = \frac{1}{n_g} \sum_{i = 1}^{n} \Bigg[ & \boldsymbol{1}\left(G_i = g\right) \left(Y_{it} - {m}_{i, t\mid g}\right) \\
& - \sum_{t' = g}^{t}\boldsymbol{1}\left(G_i > t'\right) \frac{{\pi}_{ig}}{1- {\pi}_{ig}} \left(\prod_{r = g+1}^{t'} \frac{1}{1 - {\pi}_{ir}}\right)\left( {m}_{i, t\mid t'+1} - {m}_{i, t\mid t'} \right) \Bigg].
\end{aligned}
\end{equation}
Here, (i) $n_g:=\sum_{i=1}^{n}\boldsymbol{1}(G_i=g)$ denotes the number of units first treated in period $g$; (ii) $m_{i,t\mid t+1}:=Y_{it}$ is defined for notational convenience; and (iii) an empty product is defined as one, so that $\prod_{r=g+1}^{t'}(1-\pi_{ir})^{-1}=1$ when $t'=g$. This is an oracle estimator in the sense that it uses the true propensity scores and outcome regressions, which are unknown in practice. As in the simultaneous treatment case, I use Nadaraya-Watson estimators based on the pseudo-distance to nonparametrically estimate $m_{i,t\mid t'}$ and $\pi_{it'}$, and then employ a \emph{double} cross-fitting procedure to achieve a faster convergence rate.
The full algorithm, similar to the identification strategy, involves a backward procedure and is tedious to present in the main text. Therefore, a simplified version that abstracts from cross-validation for selecting the bandwidths ${h_m,h_\pi}$ is summarized in Algorithm~\ref{alg:extension_basic} to highlight the main idea. The cross-validation procedure for bandwidth selection is discussed in the following text, and the full algorithm is presented in Algorithm~\ref{alg:extension_cv}.
\subsection{Estimation}
The convergence rate of the pseudo-distance and its informativeness under staggered adoption are identical to those under simultaneous treatment (see Assumptions~\ref{assumption:estimation_consistency_d} and \ref{assumption:informativeness_estimation}), except that $Y_{it}(0)$ needs to be replaced by $Y_{it}(\infty)$ to accommodate dynamic selection. Therefore, I omit the discussion of the accuracy of the sampled pseudo-distance and informativeness and focus directly on the inference on the dynamic ATT.
\begin{assumption}[Estimation]\label{assumption:estimation_extension} I assume that
\begin{enumerate}[label=(\roman*)]
\item \label{item:estimation_extension_panel_weak_dependent} \textbf{(Sampling)} Conditional on $\{\gamma_t\}_{t \in \mathbb{Z}}$, the panel $\{(Y_{it}, X_{it}, D_{it}, \alpha_i)\}_{i=1, \ldots, n, t = 1, \ldots, T}$ is i.i.d. across $i$. The panel is large with $n\rightarrow\infty$ and the number of pretreatment periods $T_0\rightarrow\infty$.
\item \label{item:estimation_extension_finite} \textbf{(Bounded)} $Y_{it}(\infty)$ is uniformly bounded for all $i=1, \ldots, n$ and all $t= 1, \ldots, T$. In addition, for each $t \geq g \geq T_0$, $\mathbb{E}_T\left(|Y_{it}(g)|^{2+\delta}\right)<\infty$ for some $\delta >0$.
\item \label{item:estimation_extension_compact} \textbf{(Compact)} The support of $\alpha$, $\mathcal{A}$, is compact. Let $\rho_1$ denote the radius of $\mathcal{A}$. There exist constants $\underline{c}_1, \overline{c}_1 >0$, such that for any $\alpha_0 \in \mathcal{A}$ and any $ r \in (0, \rho_1]$, $\underline{c}_1 r^{d_{\alpha}} \leq \mathbb{P}\left(\|\alpha - \alpha_0\| \leq r \right) \leq \overline{c}_1 r^{d_{\alpha}}$.
In addition, for each $t = T_0, \ldots, T$, let $\mathcal S_{t-1}$ be the support of $(X_{it},\alpha_i)$ conditional on $\{G_i > t-1\}$, and let $\rho_2$ to denote the radius of $\mathcal{S}_{t-1}$.
There exist constants $0<\underline{c}_2\leq\overline {c}_2<\infty$ such that, for any $t = T_0, \ldots, T$ and any $(x,\alpha_0)\in\mathcal{S}_{t-1} $,
\begin{align*}
\underline{c}_2 r^{d_X + d_{\alpha}} \leq \mathbb{P}\left(\|X_{it}-x \| \leq r, \|\alpha_i-\alpha_0\|\leq r, G_i > t-1 \right)\leq \overline c_2 r^{d_X + d_{\alpha}}.
\end{align*}
\item \label{item:estimation_extension_L_continuous} \textbf{(Lipschitz continuity)} For each $T_0 \leq t'\leq t\leq T$, let $\mathcal S_{t'-1}$ be the support of $(X_{it'},\alpha_i)$ conditional on $\{G_i >t'-1\}$. The functions $m_{t \mid t'}(x, \alpha)$ and $\pi_{t'}(x, \alpha)$ are Lipschitz continuous in $(x, \alpha) \in \mathcal{S}_{t'-1}$.
\item \label{item:estimation_extension_overlap} \textbf{(Common overlap)} There exist constants $\underline{p}, \overline{p}\in (0, 1)$ such that for any $t = T_0, \ldots, T$, $\underline{p} \leq \pi_{t}(X_{it}, \alpha_i)\leq \overline{p}$ a.s.
\item \label{item:estimation_extension_kernel} \textbf{(Kernel function)} The kernel $K: \mathbb{R}\mapsto \mathbb{R}_{+}$ is bounded by $\overline{K} >0$ and supported on $[-1, 1]$. In addition, $K(0)>0$ and $K(\cdot)$ is Lipschitz continuous with constant $L_K>0$.
\end{enumerate}
\end{assumption}
Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_panel_weak_dependent} and Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_kernel} are identical to those conditions in the simultaneous treatment case. Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_finite} impose a finite $(2 + \delta)$-moment condition on $Y_{it}(g)$ to accommodate dynamic selection. In addition, Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_compact} extends Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_compact} to staggered adoption.
Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_L_continuous} is similar to Assumption~\ref{assumption:estimation_simultaneous}\ref{item:estimation_simultaneous_L_continuous}, but extends
continuity condition to the dynamic counterfactual functions $m_{t\mid t'}(\cdot,\cdot)$ and the propensity scores at different times. Lastly, Assumption~\ref{assumption:estimation_extension}\ref{item:estimation_extension_overlap} is a common overlap condition imposed for treatment adoption at different periods.
\begin{proposition}[Convergence rates for counterfactual imputation]\label{prop:estimation_entry_extension}
Under Assumptions~\ref{assumption:estimation_consistency_d}-\ref{assumption:estimation_simultaneous}, \ref{assumption:selection_extension}-\ref{assumption:potential_outcome_extension}, and~\ref{assumption:estimation_extension}, for any $h\in\{h_{\pi},h_m\}$ such that $h\rightarrow 0$, $nh^{d_\alpha+d_{X}}/\log n\rightarrow\infty$, and $\delta_{n,T_0}/h\rightarrow 0$, we have that for each $T_0 \leq t'\leq t\leq T$,
\begin{align*}
\max_{i: G_i >t'-1}\left|\widehat{m}_{i, t\mid t'} - m_{i, t\mid t'} \right| = O_p\left(h_m + \left(nh_m^{d_\alpha+d_{X}}\right)^{-1/2}\sqrt{\log (n)} \right).
\end{align*}
In addition, for each $t = T_0, \ldots, T$,
\begin{align*}
\max_{i: G_i >t-1}\left|\widehat{\pi}_{it} - \pi_{it} \right| = O_p\left(h_{\pi} + \left(nh_{\pi}^{d_\alpha+d_{X}}\right)^{-1/2}\sqrt{\log (n)} \right).
\end{align*}
\end{proposition}
Proposition~\ref{prop:estimation_entry_extension} establishes the convergence rates for counterfactual outcome estimators and propensity score estimators. The convergence rates are identical to those in the simultaneous treatment case.
\begin{theorem}[Asymptotic normality]\label{thm:estimation_ATT_extension}
Under conditions in Proposition~\ref{prop:estimation_entry_extension}, for each $T_0 \leq g \leq t\leq T$,
\begin{align*}
\widehat{\mathrm{ATT}}(g, t) - \mathrm{ATT}(g, t) = \frac{1}{\sqrt{n}} \mathcal{N}(0, V_{g, t}) + O_p\left(h_m \left(h_{\pi} + (nh_{\pi}^{d_{\alpha} + d_{X}})^{-1/2}\sqrt{\log(n)}\right)\right).
\end{align*}
Here, $\mathcal{N}(\cdot,\cdot)$ denotes the normal distribution, and $V_{g, t}$ is the asymptotic variance of the oracle estimator.
\end{theorem}
Theorem~\ref{thm:estimation_ATT_extension} shows that, similar to the simultaneous treatment case, combining double cross-fitting with undersmoothing of $h_m$ yields a faster convergence rate than standard single cross-fitting.
\begin{corollary}[Inference]\label{corollary:estimation_ATT_extension}
Under conditions in Proposition~\ref{prop:estimation_entry_extension}, set $h_\pi \asymp n^{-1/(d_{\alpha} + d_{X} + 2)}$ and $h_m \asymp n^{-1/(d_{\alpha} + d_{X})}(\log (n))^{2/(d_X + d_{\alpha})} $. In addition, suppose $\delta_{n,T_0}/h_m \rightarrow 0$. Then, for each $T_0\leq g \leq t\leq T$ and $d_{\alpha}+d_{X} \leq 3$,
\begin{align*}
\sqrt n\left(\widehat{\mathrm{ATT}(g, t)}-\mathrm{ATT}(g, t)\right) \stackrel{d}{\longrightarrow} \mathcal{N}(0, V_{g, t}).
\end{align*}
Here, $\mathcal{N}(\cdot,\cdot)$ denotes the normal distribution, and $V_{g, t}$ is the asymptotic variance of the oracle estimator.
\end{corollary}
Corollary~\ref{corollary:estimation_ATT_extension} shows that, when $d_{\alpha}+d_{X}\leq 3$, combining double cross-fitting with undersmoothing yields a root-$n$ asymptotically normal and asymptotically unbiased estimator of the dynamic ATT, enabling valid inference.
\subsection{Bandwidth selection in practice}
The theoretical bandwidth choice in Corollary~\ref{corollary:estimation_ATT_extension} provides practical guidance for selecting $h_m$ and $h_\pi$ using data-driven methods. Similar to the simultaneous treatment case, the choice of $h_\pi$ optimizes the MSE of the propensity score estimator, while $h_m$ is aggressively undersmoothed.
I propose to use standard leave-one-out cross-validation for selecting $h_\pi$ and the effective number of matched neighbors for selecting $h_m$. The bandwidth selection procedure is almost identical to that in the simultaneous treatment case, except that it is implemented separately for each period. The full procedure is lengthy and is therefore omitted from the main text. The complete algorithm, including the bandwidth selection procedure, is provided in Algorithm~\ref{alg:extension_cv}.
\section{Simulation}\label{sec:simulation}
In this section, I illustrate the finite-sample properties of the proposed method under different DGPs, including dynamic models with additive fixed effects, dynamic models with interactive fixed effects, and nonlinear dynamic models.
I report three estimators. (i) \textbf{DID}: the canonical difference-in-differences estimator. (ii) $\boldsymbol{\mathrm{DR2}}$: the doubly robust estimator based on a standard two-way sample split, with $m$ and $\pi$ estimated using the same training sample. The bandwidths for $m$ and $\pi$ are selected by leave-one-out cross-validation to minimize mean squared errors. I also use $\mathrm{DR2}^*$ to denote the infeasible estimator that uses the true distance between individual latent factors rather than the pseudo-distance. (iii) $\boldsymbol{\mathrm{DR3}}$: the proposed estimator using a three-way sample split, with $m$ and $\pi$ estimated on separate training samples. I also use $\mathrm{DR3}^*$ to denote the infeasible estimator that uses the true distance between individual latent factors rather than the pseudo-distance. I vary $\kappa$ to examine how different choices affect the estimator's performance.
\subsection{Dynamic panel with additive fixed effects}
Consider the following data-generating process:
\begin{equation}\label{eq:sim_TWFE_Y_simultaneous}
\begin{aligned}
Y_{it} =
\left\{
\begin{array}{lr}
\rho Y_{it-1} + \alpha_i + \gamma_t + \epsilon_{Y, it}, & t < T_0 \\
\rho Y_{it-1} + \alpha_i + \gamma_t + 0.5 \cdot D_i + \epsilon_{Y, it}, & t \geq T_0 \\
\end{array}
\right..
\end{aligned}
\end{equation}
The individual latent factors $\{\alpha_i\}_{i=1}^{n}$ and the time latent factors $\{\gamma_t\}_{t\leq T}$ are independent random variables drawn from the uniform distribution on $[-1/4,1/4]$. The error terms $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$ are independent of the latent factors and are i.i.d. across both dimensions, following the uniform distribution on $[-1/2,1/2]$. I set $\rho=0.8$ to capture the dynamic effect. The contemporaneous ATT is $0.5$, and for each $t>T_0$, the dynamic treatment effect is $0.5\left(1+\rho+\ldots+\rho^{t-T_0}\right)$. Treatment assignment depends jointly on latent heterogeneity and the lagged outcome:
\begin{align}\label{eq:sim_TWFE_selection_simultaneous}
D_i = \boldsymbol{1}\left( \alpha_i/2 + Y_{i T_0 - 1}/2 + \epsilon_{D, i} \geq 0\right).
\end{align}
Here, the error terms $\{\epsilon_{D,i}\}_{i=1}^{n}$ are i.i.d. across individuals and are independent of the latent factors and $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$. In the numerical designs, $\epsilon_{D,i}$ follows a standard logistic distribution. I vary the sample sizes $(N,T_0)\in\left\{(1000,20),(200,100),(1000,100),(2000,100)\right\}$ and perform $5{,}000$ replications for each design.
\begin{table}[h]
\centering
\caption{Simulation Results: Dynamic Panel with Additive Fixed Effects}\label{tab:TWFE_simultaneous}
\begin{tabular}{cccccccccc}
\toprule
& DID & $\mathrm{DR2}^*$ & DR2 & \multicolumn{3}{c}{$\mathrm{DR3}^*$ } & \multicolumn{3}{c}{DR3} \\
\cmidrule(lr){5-7} \cmidrule(lr){8-10} & & & &$\kappa=.15$ & {$\kappa=.2$} & {$\kappa=.25 $} &$\kappa=.15$ & {$\kappa=.2$} & {$\kappa=.25 $} \\
\midrule
$N= 1000, T_0= 20$ & & & & & & & & & \\
Bias & 2.16 & 0.54 & 0.84 & 0.15 & 0.21 & 0.32 & 0.47 & 0.54 & 0.64 \\
SD & 1.94 & 2.07 & 2.15 & 2.51 & 2.41 & 2.30 & 2.54 & 2.44 & 2.37 \\
Coverage & 79.8 & 94.0 & 92.1 & 93.7 & 93.7 & 93.8 & 93.7 & 93.7 & 93.1 \\
\addlinespace[5pt]
$N= 200, T_0= 100$ & & & & & & & & & \\
Bias & 2.16 & 0.95 & 1.20 & 0.60 & 0.64 & 0.79 & 0.80 & 0.83 & 1.00 \\
SD & 1.94 & 5.38 & 5.45 & 7.05 & 7.03 & 6.94 & 7.18 & 7.15 & 7.01 \\
Coverage & 79.8 & 93.9 & 93.7 & 93.3 & 93.2 & 92.7 & 92.7 & 92.8 & 92.4 \\
\addlinespace[5pt]
$N= 1000, T_0= 100$ & & & & & & & & & \\
Bias & 2.16 & 0.57 & 0.78 & 0.22 & 0.26 & 0.36 & 0.41 & 0.49 & 0.59 \\
SD & 1.94 & 2.05 & 2.08 & 2.49 & 2.39 & 2.32 & 2.51 & 2.41 & 2.33 \\
Coverage & 79.8 & 93.6 & 93.5 & 94.2 & 93.9 & 93.5 & 93.4 & 93.3 & 93.5 \\
\addlinespace[5pt]
$N= 2000, T_0= 100$ & & & & & & & & & \\
Bias & 2.16 & 0.38 & 0.57 & 0.10 & 0.14 & 0.18 & 0.30 & 0.37 & 0.43 \\
SD & 1.94 & 1.43 & 1.45 & 1.68 & 1.62 & 1.58 & 1.71 & 1.62 & 1.59 \\
Coverage & 79.8 & 93.8 & 92.8 & 93.2 & 94.1 & 93.8 & 94.4 & 93.7 & 93.2 \\
\bottomrule
\end{tabular}
\vspace{0.3cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} The table reports Monte Carlo results for the contemporaneous ATT based on $5000$ replications of the dynamic panel model with additive fixed effects specified in~\eqref{eq:sim_TWFE_Y_simultaneous} and the treatment selection mechanism specified in~\eqref{eq:sim_TWFE_selection_simultaneous}. I fix $\underline{q}=80\%$ and vary $\kappa \in \{.15,.2,.25\}$ when selecting the bandwidth $h_m^*$ according to Algorithm~\ref{alg:simultaneous_cv}. I report the bias, standard deviation, and coverage rate for the difference-in-differences estimator ($\mathrm{DID}$), the infeasible doubly robust estimator based on a standard two-way sample split and the true latent distance ($\mathrm{DR2}^*$), and its feasible counterpart based on the estimated pseudo-distance ($\mathrm{DR2}$). I also report the bias, standard deviation, and coverage rate for the infeasible estimator based on a three-way sample split and the true latent distance ($\mathrm{DR3}^*$), as well as the proposed estimator based on a three-way sample split and the estimated pseudo-distance ($\mathrm{DR3}$).
Biases and standard deviations are reported in units of $\boldsymbol{0.01}$, and coverage rates are reported in percentage points $\boldsymbol{1\%}$.
\end{minipage}
\end{table}
Table~\ref{tab:TWFE_simultaneous} reports the simulation results for the contemporaneous ATT and shows that the proposed estimator performs well across all designs. Relative to DID and the standard doubly robust estimator, the proposed estimator substantially reduces bias and has coverage close to the nominal 95\% level. In addition, the proposed estimator using the pseudo-distance performs similarly to the infeasible DR3$^*$ using the true latent distance, indicating that the estimated pseudo-distance provides a good approximation to the true latent distance even when $T_0=20$. It is worth noting that, when $N=200$, DR2 has a smaller standard deviation and coverage closer to the nominal level, although the proposed estimator has a smaller absolute bias. This is because the double cross-fitting procedure leaves relatively few observations in each subsample. Finally, the performance of the proposed estimator is also stable across different values of $\kappa$, and I recommend $(\kappa,\underline{q})=(0.2,0.8)$ as the default choice.
\subsection{Dynamic panel with interactive fixed effects}
Consider the following data-generating process:
\begin{equation}\label{eq:sim_IFE_Y_simultaneous}
\begin{aligned}
Y_{it} =
\left\{
\begin{array}{lr}
\rho Y_{it-1} + \alpha_i \gamma_t + \epsilon_{Y, it}, & t < T_0 \\
\rho Y_{it-1} + \alpha_i \gamma_t + 0.5 \cdot D_i + \epsilon_{Y, it}, & t \geq T_0 \\
\end{array}
\right..
\end{aligned}
\end{equation}
The individual latent factors $\{\alpha_i\}_{i=1}^{n}$ and the time latent factors $\{\gamma_t\}_{t\leq T}$ are independent random variables drawn from uniform distributions on $[-1,1]$ and $[-2,2]$, respectively. The error terms $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$ are independent of the latent factors and are i.i.d. across both dimensions, following a uniform distribution on $[-1/2,1/2]$. I set $\rho=0.8$ to capture the dynamic effect. The contemporaneous ATT is $0.5$, and for each $t>T_0$, the dynamic treatment effect is $0.5\left(1+\rho+\ldots+\rho^{t-T_0}\right)$. Treatment assignment depends jointly on latent heterogeneity and the lagged outcome:
\begin{align}\label{eq:sim_IFE_selection_simultaneous}
D_i = \boldsymbol{1}\left( \alpha_i/2 + Y_{i T_0 - 1}/2 + \epsilon_{D, i} \geq 0\right).
\end{align}
Here, the error terms $\{\epsilon_{D,i}\}_{i=1}^{n}$ are i.i.d. across individuals and are independent of the latent factors and $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$. In the numerical designs, $\epsilon_{D,i}$ follows a standard logistic distribution. I vary the sample sizes $(N,T_0)\in\left\{(1000,20),(200,100),(1000,100),(2000,100)\right\}$ and perform $5000$ replications for each design.
\begin{table}[h]
\centering
\caption{Simulation Results: Dynamic Panel with Interactive Fixed Effects}\label{tab:IFE_simultaneous}
\begin{tabular}{ccccccccc}
\toprule
& $\mathrm{DR2}^*$ & DR2 & \multicolumn{3}{c}{$\mathrm{DR3}^*$ } & \multicolumn{3}{c}{DR3} \\
\cmidrule(lr){4-6} \cmidrule(lr){7-9} & & &$\kappa=.15$ & {$\kappa=.2$} & {$\kappa=.25$} & {$\kappa=.15$} & {$\kappa=.2$} & {$\kappa=.25$} \\
\midrule
$N= 1000, T_0= 20$ & & & & & & & & \\
Bias & 0.60 & 0.82 & 0.24 & 0.34 & 0.46 & 0.46 & 0.57 & 0.69 \\
SD & 2.27 & 2.56 & 2.70 & 2.61 & 2.53 & 3.00 & 2.90 & 2.82 \\
Coverage & 93.5 & 92.1 & 93.6 & 93.6 & 93.5 & 93.1 & 92.8 & 92.3 \\
\addlinespace[5pt]
$N= 200, T_0= 100$ & & & & & & & & \\
Bias & 1.06 & 1.16 & 0.86 & 0.90 & 1.07 & 0.91 & 0.94 & 1.18 \\
SD & 6.11 & 6.29 & 8.12 & 8.09 & 7.97 & 8.29 & 8.25 & 8.17 \\
Coverage & 93.9 & 93.8 & 92.9 & 92.6 & 92.4 & 92.8 & 92.8 & 92.7 \\
\addlinespace[5pt]
$N= 1000, T_0= 100$ & & & & & & & & \\
Bias & 0.62 & 0.75 & 0.30 & 0.36 & 0.50 & 0.41 & 0.50 & 0.62 \\
SD & 2.29 & 2.35 & 2.74 & 2.65 & 2.58 & 2.83 & 2.74 & 2.68 \\
Coverage & 93.3 & 92.7 & 93.7 & 93.5 & 93.1 & 93.8 & 93.4 & 92.6 \\
\addlinespace[5pt]
$N= 2000, T_0= 100$ & & & & & & & & \\
Bias & 0.40 & 0.52 & 0.13 & 0.20 & 0.27 & 0.25 & 0.32 & 0.40 \\
SD & 1.56 & 1.62 & 1.83 & 1.76 & 1.71 & 1.92 & 1.84 & 1.80 \\
Coverage & 93.7 & 92.6 & 94.2 & 93.9 & 93.7 & 93.8 & 93.3 & 93.2 \\
\bottomrule
\end{tabular}
\vspace{0.3cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} The table reports Monte Carlo results for the contemporaneous ATT based on $5000$ replications of the dynamic panel model with interactive fixed effects specified in~\eqref{eq:sim_IFE_Y_simultaneous} and the treatment selection mechanism specified in~\eqref{eq:sim_IFE_selection_simultaneous}. I fix $\underline{q}=80\%$ and vary $\kappa \in \{.15,.2,.25\}$ when selecting the bandwidth $h_m^*$ according to Algorithm~\ref{alg:simultaneous_cv}. I report the bias, standard deviation, and coverage rate for the infeasible doubly robust estimator based on a standard two-way sample split and the true latent distance ($\mathrm{DR2}^*$), and its feasible counterpart based on the estimated pseudo-distance ($\mathrm{DR2}$). I also report the bias, standard deviation, and coverage rate for the infeasible estimator based on a three-way sample split and the true latent distance ($\mathrm{DR3}^*$), as well as the proposed estimator based on a three-way sample split and the estimated pseudo-distance ($\mathrm{DR3}$).
Biases and standard deviations are reported in units of $\boldsymbol{0.01}$, and coverage rates are reported in percentage points $\boldsymbol{1\%}$.
\end{minipage}
\end{table}
Table~\ref{tab:IFE_simultaneous} reports the simulation results for the contemporaneous ATT under interactive fixed effects (the DiD estimator is omitted because the DGP does not admit an additive two-way fixed-effects representation). The table shows that the proposed estimator reduces bias and achieves coverage close to the nominal 95\% level. Also, the proposed estimator performs similarly to the infeasible DR3$^*$ estimator, indicating that the estimated pseudo-distance provides a good approximation to the true latent distance even when $T_0=20$ under the interactive fixed effects model. When $N=200$, the performance of the proposed estimator is exceeded by that of the standard doubly robust estimator because the double cross-fitting procedure leaves relatively few observations in each subsample. Still, I recommend $(\kappa,\underline{q})=(0.2,0.8)$ as the default choice.
\subsection{Nonlinear dynamic panel with fixed effects}
Consider the following data-generating process:
\begin{equation}\label{eq:sim_NL_Y_simultaneous}
\begin{aligned}
Y_{it} =
\left\{
\begin{array}{lr}
\boldsymbol{1}\left(-2 + \rho Y_{it-1} + 4 \alpha_i + \epsilon_{Y, it} \geq 0 \right), & t < T_0 \\
\boldsymbol{1}\left(-2 + \rho Y_{it-1} + 4 \alpha_i + D_i + \epsilon_{Y, it} \geq 0 \right), & t \geq T_0 \\
\end{array}
\right..
\end{aligned}
\end{equation}
The individual latent factors $\{\alpha_i\}_{i=1}^{n}$ consist of independent random variables drawn from the uniform distribution on $[-1,1]$. The error terms $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$ are independent of the latent factors and are i.i.d. across both dimensions, following the standard logistic distribution. I set $\rho=0.5$ to capture dynamic dependence. The contemporaneous ATT, computed using numerical simulation, is approximately $0.12$. Dynamic treatment effects for $t>T_0$ can also be calculated numerically. Treatment assignment depends jointly on latent heterogeneity and the lagged outcome:
\begin{align}\label{eq:sim_NL_selection_simultaneous}
D_i = \boldsymbol{1}\left( \alpha_i/2 + Y_{i T_0 - 1}/2 + \epsilon_{D, i} \geq 0\right).
\end{align}
Here, the error terms $\{\epsilon_{D,i}\}_{i=1}^{n}$ are i.i.d. across individuals and are independent of the latent factors and $\{\epsilon_{Y,it}\}_{i=1,\ldots,n,\;t=0,\ldots,T}$. In the numerical designs, $\epsilon_{D,i}$ follows the standard logistic distribution. I vary the sample sizes $(N,T_0)\in\left\{(1000,20),(200,100),(1000,100),(2000,100)\right\}$ and perform $5000$ replications for each design.
\begin{table}[htbp]
\centering
\caption{Simulation Results: Nonlinear Dynamic Panel with Fixed Effects}\label{tab:NL_simultaneous}
\begin{tabular}{ccccccccc}
\toprule
& $\mathrm{DR2}^*$ & DR2 & \multicolumn{3}{c}{$\mathrm{DR3}^*$ } & \multicolumn{3}{c}{DR3} \\
\cmidrule(lr){4-6} \cmidrule(lr){7-9} & & &$\kappa=.15$ & {$\kappa=.2$} & {$\kappa=.25$} & {$\kappa=.15$} & {$\kappa=.2$} & {$\kappa=.25$} \\
\midrule
$N= 1000, T_0= 20$ & & & & & & & & \\
Bias & 0.20 & 1.10 & -0.31 & -0.24 & -0.02 & -0.25 & -0.05 & 0.38 \\
SD & 2.51 & 2.67 & 3.21 & 3.17 & 2.97 & 3.08 & 3.07 & 3.04 \\
Coverage & 95.1 & 93.1 & 94.2 & 94.3 & 94.9 & 94.0 & 94.1 & 94.1 \\
\addlinespace[5pt]
$N= 200, T_0= 100$ & & & & & & & & \\
Bias & 0.38 & 0.77 & -0.48 & -0.47 & -0.38 & -0.54 & -0.50 & -0.33 \\
SD & 6.72 & 6.96 & 9.90 & 9.89 & 9.83 & 10.01 & 10.01 & 10.01 \\
Coverage & 94.3 & 94.5 & 93.6 & 93.4 & 93.2 & 92.4 & 92.4 & 92.6 \\
\addlinespace[5pt]
$N= 1000, T_0= 100$ & & & & & & & & \\
Bias & 0.21 & 0.63 & -0.31 & -0.25 & -0.06 & -0.09 & -0.04 & 0.06 \\
SD & 2.56 & 2.58 & 3.19 & 3.14 & 2.99 & 3.08 & 3.06 & 3.03 \\
Coverage & 95.1 & 94.5 & 94.5 & 94.6 & 94.2 & 93.9 & 94.1 & 94.1 \\
\addlinespace[5pt]
$N= 2000, T_0= 100$ & & & & & & & & \\
Bias & 0.12 & 0.55 & -0.06 & -0.06 & -0.06 & -0.05 & -0.04 & 0.00 \\
SD & 1.76 & 1.77 & 2.06 & 2.06 & 2.06 & 2.06 & 2.05 & 2.05 \\
Coverage & 95.1 & 94.5 & 94.8 & 94.8 & 94.6 & 94.2 & 94.3 & 94.2 \\
\bottomrule
\end{tabular}
\vspace{0.3cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} The table reports Monte Carlo results for the contemporaneous ATT based on $5000$ replications of the nonlinear dynamic panel model with fixed effects specified in~\eqref{eq:sim_NL_Y_simultaneous} and the treatment selection mechanism specified in~\eqref{eq:sim_NL_selection_simultaneous}. I fix $\underline{q}=80\%$ and vary $\kappa \in \{.15,.2,.25\}$ when selecting the bandwidth $h_m^*$ according to Algorithm~\ref{alg:simultaneous_cv}. I report the bias, standard deviation, and coverage rate for the infeasible doubly robust estimator based on a standard two-way sample split and the true latent distance ($\mathrm{DR2}^*$), and its feasible counterpart based on the estimated pseudo-distance ($\mathrm{DR2}$). I also report the bias, standard deviation, and coverage rate for the infeasible estimator based on a three-way sample split and the true latent distance ($\mathrm{DR3}^*$), as well as the proposed estimator based on a three-way sample split and the estimated pseudo-distance ($\mathrm{DR3}$).
Biases and standard deviations are reported in units of $\boldsymbol{0.01}$, and coverage rates are reported in percentage points $\boldsymbol{1\%}$.
\end{minipage}
\end{table}
Table~\ref{tab:NL_simultaneous} reports the simulation results for the contemporaneous ATT under the nonlinear dynamic panel model (the DiD estimator is omitted because the DGP does not admit an additive two-way fixed-effects representation). The table shows that both the standard doubly robust estimator and the proposed estimator have small biases and achieve coverage close to the nominal 95\% level. The standard doubly robust estimator performs particularly well in this design because $Y_{it}$ is discrete, reducing the effective dimension of the nonparametric regression to one and allowing the estimator to be root-$N$ consistent and asymptotically unbiased \citep{deaner2025inferring}. The proposed estimator performs similarly to the infeasible DR3$^*$ estimator, indicating that the estimated pseudo-distance provides a good approximation to the true latent distance even when $T_0=20$. Finally, I recommend $(\kappa,\underline{q})=(0.2,0.8)$ as the default choice.
\bigskip
Overall, the proposed estimator substantially reduces bias and achieves coverage rates close to the nominal level. It performs well across different DGPs and sample sizes. In particular, its performance remains strong when the number of pretreatment periods is moderate, such as $T_0=20$. This result is especially relevant because panels with a very large time dimension are uncommon in empirical applications. The simulations show that the performance of the proposed estimator is more sensitive to the cross-sectional sample size $N$ and deteriorates when $N$ is small.
\section{Empirical Application}\label{sec:empirical}
I apply the proposed method to the setting studied by \citet{bailey2012reexamining}, estimating the long-run effects of U.S. family planning programs on childbearing, allowing treatment assignment to depend on county fixed effects and pretreatment fertility rates.
The original data is a county-year panel covers $3,037$ counties from $1959$ to $1988$. The outcome, $Y_{it}$, is the general fertility rate (GFR), defined as the number of live births per $1000$ women aged $15-44$ in county $i$ and year $t$.
Treatment timing is summarized in Table~\ref{tab:treatment_timing} and is measured by the fiscal year in which a county first received a recorded federal family planning grant. Counties were treated sequentially between 1965 and 1973, and those without a recorded grant during this period are coded as untreated. Following \citet{bailey2012reexamining}, I pool adjacent treatment years into three groups: $1965-1967$, $1968-1969$, and $1970-1973$\footnote{
This is because estimating the ATT separately for each $g = 1965, \ldots, 1973$ is infeasible in practice: the number of counties first treated in some years (e.g., 1965 and 1973) is very small, making the corresponding nonparametric estimates unstable and imprecise.
In addition, since the analysis focuses on medium- and long-run effects, pooling adjacent cohorts sacrifices little timing variation while preserving the staggered rollout of the program.
}.
\begin{table}[H]
\centering
\caption{Treatment Timing}
\label{tab:treatment_timing}
\begin{tabular*}{0.85\textwidth}{@{\extracolsep{\fill}}lccc@{}}
\toprule
& Years & No. treated by year & Total \\
\midrule
Group 1 & $1965-1967$ & 6, \; 43, \; 74 & 123 \\
Group 2 & $1968-1969$ & 52, \; 278 & 330 \\
Group 3 & $1970-1973$ & 63, \;75, \;53, \;10 & 201 \\
\bottomrule
\end{tabular*}
\end{table}
\paragraph{Pseudo distance}
To obtain a longer pretreatment outcome history, I match the Bailey
replication sample to the \textit{U.S. County-Level Natality and Mortality
Data, 1915--2007}\footnote{
The data is publicly available at \href{https://www.icpsr.umich.edu/web/ICPSR/studies/36603}{https://www.icpsr.umich.edu/web/ICPSR/studies/36603}.
} assembled by \citet{BaileyEtAl2016Data}. This database provides
county-year live births by the mother's county of residence and estimates of
the female population aged 15-44 since $1937$. This enables me to construct the historical GFR from $1937$ to $1958$\footnote{
For the overlapping period $1959$-$1988$, I just retain the GFR reported in the
\citet{bailey2012reexamining} replication files rather than replacing it with the reconstructed
series.
}.
The final panel contains $3,017$ county units observed from $1937$ to $1988$ with $28$ years of pretreatment outcomes for each county\footnote{
Counties that could not be matched or lacked the required outcome history were excluded from the final sample. For example, Los Alamos was excluded because its outcome data are unavailable before $1950$.
}.
I use the pretreatment outcomes from $1937-1964$ to calculate pseudo-distance between each pair of counties. Although someone may argue that there are structural changes in the series because the period span the Second World War and the baby boom, this does not affect the construction when the structural change can be captured by an additive common time trend, as the common trend is differenced out by cross-sectional differencing. The results are similar when I use data from $1946-1964$ to calculate the pseudo-distance.
\paragraph{Selection} Selection into the federal family planning program may depend on both persistent
county-level heterogeneity and lagged outcomes. Although the early federal grant-making process operated under substantial pressure, with funds disbursed on a first-come, first-served basis (\citet[pp.~70--71]{bailey2012reexamining}), whether and when a local organization entered the applicant pool remained endogenous. In fact, \citet{bailey2012reexamining} documents that funded counties were larger and more urban and, on average, more educated and affluent than unfunded counties, which supports the selection on county fixed effects.
In addition, selection may depend on pretreatment fertility rates. For example, local authorities in areas with higher fertility rates may have had greater demand for subsidized family planning services and, consequently, stronger incentives to apply. On the other hand, local authorities that were already more supportive of contraception before treatment may have applied earlier, and some applicants may have initiated local family planning programs before receiving federal funding, both of which could generate a negative relationship between pretreatment fertility rates and selection into the program.
To investigate this concern, \citet{bailey2012reexamining} uses linear regressions of treatment year on lagged outcomes and finds small, statistically insignificant coefficients. However, these results do not fully rule out selection based on pretreatment fertility rates, as the regressions remain subject to omitted-variable bias and functional-form misspecification.
\begin{table}[h]
\centering
\caption{Selection on lagged outcomes}
\begin{tabular}{lcccccc}
\toprule
& \multicolumn{2}{c}{ Treated $1965-67$}
& \multicolumn{2}{c}{Treated $1968-69$}
& \multicolumn{2}{c}{Treated $1970-73$} \\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}
& (1) & (2) & (1) & (2) & (1) & (2) \\
\bottomrule
& & & & & & \\
State fixed effects & $0.24$ & $0.20$ & $1.07^{***}$ & $1.16^{***}$ & $0.35$ & $0.23$ \\
& $(0.45)$ & $(0.49)$ & $(0.34)$ & $(0.36)$ & $(0.48)$ & $(0.53)$ \\
& & & & & & \\
$+$ Urban & $0.71$ & $0.65$ & $1.33^{***}$ & $1.50^{***}$ & $0.92^*$ & $0.80$\\
& $(0.51)$ & $(0.53)$ & $(0.37)$ & $(0.38)$ & $(0.54)$ & $(0.57)$ \\
& & & & & & \\
$+$ Education & $0.50$ & $0.43$ & $1.18^{***}$ & $1.35^{***}$ & $1.07^{**}$ & $1.00^{*}$ \\
& $(0.53)$ & $(0.55)$ & $(0.39)$ & $(0.40)$ & $(0.55)$ & $(0.59)$ \\
& & & & & & \\
$+$ Income & $0.33$ & $0.27$ & $1.08^{***}$ & $1.29^{***}$ & $1.07^{**}$ & $1.01^{*}$ \\
& $(0.54)$ & $(0.56)$ & $(0.39)$ & $(0.40)$ & $(0.55)$ & $(0.59)$\\
& & & & & & \\
$+$ Race and population & $0.34$ & $0.30$ & $0.86^{**}$ & $0.98^{**}$ & $0.92^*$ & $0.82$ \\
& $(0.54)$ & $(0.56)$ & $(0.39)$ & $(0.41)$ & $(0.56)$ & $(0.61)$ \\
& & & & & & \\
Observations & $3037$ & $3037$ & $2914$ & $2914$ & $2584$ & $2584$ \\
\bottomrule
\end{tabular}\label{tab:empirical_selection_lag_outcomes}
\vspace{0.3cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} This table reports the estimated logistic regression coefficients on lagged outcomes for counties first funded in 1965-67, 1968-69, and 1970-73. Columns (1) and (2) use the one-period lag and the average of the three preceding lags, respectively. Controls are added sequentially across rows. The baseline specification includes state fixed effects, followed by controls for urbanization (the percentage of the population living in urban areas in 1960), education (the percentage of the population with more than 12 years of education in 1960), income (the percentage of the population with an annual income below 3,000 US dollars), race (the percentage of the population that was nonwhite in 1960), and population (in 1960). Standard errors are reported in parentheses. $***$, $**$, and $*$ denote statistical significance at the 1\%, 5\%, and 10\% levels, respectively.
\end{minipage}
\end{table}
I find strong evidence between lagged fertility rates and selection into treatment. In Table~\ref{tab:empirical_selection_lag_outcomes}, for each treatment group, I report logistic regressions of selection into treatment on either the one-period lag of the outcome or the average of its three most recent lags. Starting from state fixed effects, I sequentially add controls for urbanization, education, income, race, and population (only the coefficient on the lagged outcome is reported). The results show that lagged outcomes are strongly associated with treatment for counties treated in 1968--69, with positive and statistically significant coefficients across all specifications. The relationship is weaker for counties treated in 1970-73. Although the coefficients are not statistically significant for the 1965--67 group, this does not rule out selection based on pretreatment fertility rates because the two opposing selection mechanisms described above can offset each other. This evidence motivates the use of the proposed method to allow treatment assignment to depend both on county fixed effects and on lagged fertility rates.
\paragraph{Estimation}
For estimation, I control for unobserved time-invariant heterogeneity using the pseudo-distance and for observed lagged fertility rates. Specifically, I use
\begin{align*}
X_{it}
= \left(Y_{it-1}+Y_{it-2}+Y_{it-3}\right) / 3
\end{align*}
to capture the possibility that treatment assignment depends on a county's lagged fertility history\footnote{
The average is calculated using fertility rates from the three years preceding the earliest treatment year in that group. For example, for the $1965-1967$ group, I use the average fertility rate over $1962-1964$.}.
I set $\kappa=0.2$, as suggested by the simulations. The estimates in \citet{bailey2012reexamining} are based on event-study regressions and are therefore potentially affected by negative weights, so I also report estimates obtained using the DiD estimator of~\citet{callaway2021difference}.
For each treatment group, I report ATT estimates beginning in the year after the latest treatment year in that group. For example, for the $1965-1967$ group, the reported ATT estimates begin in 1968.
\begin{figure}[H]
\centering
\includegraphics[width=1.0\textwidth]
{figures/weighted_ATT_comparison.pdf}
\caption{Weighted dynamic ATT.}
\label{fig:att}
\vspace{0.3cm}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} Following \citet{bailey2012reexamining}, the results are weighted by the number of women aged $15-44$ in $1970$.
\end{minipage}
\end{figure}
\paragraph{Results}
Figure~\ref{fig:att} presents the dynamic ATT estimates and 95\% confidence intervals obtained using the proposed method and the DiD estimator of~\citet{callaway2021difference}. The first row displays the dynamic ATT estimates from the proposed method. The estimates show that the family planning program has a statistically significant negative effect on fertility. In addition, the long-run effect is statistically significant for the group treated in 1968-69. The DiD estimates in the second row show similar results.
The last row of Figure~\ref{fig:att} compares estimates from the proposed method and the DiD estimator over the first $15$ post-treatment periods. The proposed method yields larger estimated reductions in fertility for the first two treatment groups, whereas the estimates are similar for the third group. For the first two groups, the gap becomes larger at longer post-treatment horizons. These results highlight the importance of accounting jointly for selection on lagged outcomes and their dynamic effects when evaluating policies.
\begin{table}[h]
\centering
\caption{ATT by Post-Treatment Horizon}
\begin{tabular}{lccc}
\toprule
& \multicolumn{1}{l}{ATT (1 - 5)} & \multicolumn{1}{l}{ATT (6 - 10)} & \multicolumn{1}{l}{ATT (11 - 15)} \\
\midrule
Proposed estimator & $-1.70^{**}$ & $-3.39^{***}$ & $-3.03^{**}$ \\
& (0.67) & (1.01) & (1.39) \\
DiD Callaway-Sant'anna & $-1.37^{**}$ & $-2.63^{***}$ & $-1.42$ \\
& (0.54) & (0.80) & (0.94) \\
Difference & -0.32 & -0.76 & -1.61 \\
& (0.76) & (1.11) & (1.28) \\
\bottomrule
\end{tabular}
\label{tab:long_term_ATT}
\begin{minipage}{\textwidth}
\footnotesize
\textbf{Note:} Following \citet{bailey2012reexamining}, the results are weighted by the number of women aged $15-44$ in $1970$.
\end{minipage}
\end{table}
Table~\ref{tab:long_term_ATT} reports treatment effects averaged over three five-year windows and across the three treatment groups. The Callaway-Sant'Anna DiD estimates are similar to those reported by \citet{bailey2012reexamining}, but the proposed estimator produces more negative estimates than the Callaway-Sant'Anna DiD estimator in all three windows. The magnitude of the difference increases from 0.32 in years 1-5 to 0.76 in years 6-10 and then to 1.61 in years 11-15. Consistent with the dynamic-selection mechanism discussed above, this widening gap shows that accounting for selection based on lagged fertility rates is particularly important for estimating longer-run treatment effects.
\bibliographystyle{aea}
\bibliography{ref}
\clearpage