Inferring Treatment Effects in Large Panels by Uncovering Latent Similarities

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

118,031 characters · 11 sections · 118 citation commands

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

abstractThe presence of unobserved confounders is one of the main challenges in identifying treatment effects. In this paper, we propose a new approach to causal inference using panel data with large $N$ and $T$. Our approach imputes the untreated potential outcomes for treated units using the outcomes for untreated individuals with similar values of the latent confounders. In order to find units with similar latent characteristics, we utilize long pre-treatment histories of the outcomes. Our analysis is based on a nonparametric, nonlinear, and nonseparable factor model for untreated potential outcomes and treatments. The model satisfies minimal smoothness requirements. We impute both missing counterfactual outcomes and propensity scores using kernel smoothing based on the constructed measure of latent similarity between units, and demonstrate that our estimates can achieve the optimal nonparametric rate of convergence up to log terms. Using these estimates, we construct a doubly robust estimator of the period-specifc average treatment effect on the treated (ATT), and provide conditions, under which this estimator is $\sqrt{N}$-consistent, and asymptotically normal and unbiased. Our simulation study demonstrates that our method provides accurate inference for a wide range of data generating processes.

Introduction

The potential presence of unobserved confounding factors presents a challenging obstacle for casual inference. However, the availability of rich panel data can help solve this problem. Panel data contain multiple observations for each individual and thus it provides a possible means of controlling for time-invariant confounders (fixed effects). Moreover, panel data can enable the imputation of individual counterfactual outcomes. Methods for causal inference in panels have become dominant in empirical work, and their development is an active and rapidly growing area of research (see, e.g., arkhangelsky2024causal).

Many popular methods for causal inference using panel data assume that untreated potential outcomes have a linear/affine structure. For example, difference-in-differences (DiD) and two-way fixed effects (TWFE) methods rely on parallel trend assumptions, whereas synthetic control (SC) and most matrix completion methods require that untreated potential outcomes obey a linear factor model. The assumption of a linear model can simplify estimation and inference, and may be necessary given limited data availability. However, the assumption of linearity may be overly restrictive.

This paper presents a new method for estimation of, and inference on, the causal effects of a binary treatment $w_{i,t}$ in large panels. We assume that the untreated potential outcome of unit $i$ in period $t$ follows the possibly nonlinear and non-separable factor model below

align[align omitted — 96 chars of source]

where $\alpha_i$ and $\lambda_t$ are unobserved unit and time effects, $\mu^{(0)}(\cdot, \cdot)$ is an unknown function, and $u_{i,t}$ is an (exogenous) residual.\footnote{It is also straightforward to allow for a linear covariates adjustment in (ref). We abstract from this extension for clarity of exposition.} We assume that conditional on the latent factors, untreated potential outcomes are independent of treatment. Importantly, we do not assume that $\mu^{(0)}(\cdot, \cdot)$ has any particular parametric form. Thus, the model (ref) is substantially more general than the commonly used linear formulations $\mu^{(0)}(\alpha_i, \lambda_t) = \alpha_i + \lambda_t$ and $\mu^{(0)}(\alpha_i, \lambda_t) = \alpha_i' \lambda_t$ imposed by DiD/TWFE and SC methods, respectively.

While the generality of formulation (ref) is appealing, it obscures identification and estimation of treatment effects because $\alpha_i$ cannot be identified without strong additional restrictions are on the form of $\mu^{(0)}(\cdot, \cdot)$. Instead of trying to model $\mu^{(0)}(\cdot, \cdot)$ and to estimate $\alpha_i$ directly, we follow an approach first developed by zhang2017estimating to estimate the graphon in the context of network data. We use a long history of pre-treated outcomes to construct a pseudo-distance $\hat{d}_{i,j}$ that is informative about the closeness of the latent characteristics of individuals $i$ and $j$. That is, the closeness of $\alpha_i$ and $\alpha_j$. This allows us to impute the values of conditional mean potential outcomes and propensity scores (conditional on the latent factors) in a given post-treatment period $t$ for all units by kernel smoothing based on $\hat{d}_{ij}$. This is despite the fact that $\alpha_i$ is neither observed nor identified.

We construct estimates of $\mu^{(0)}_{i,t}$ and $p_{i,t}$ and establish their rates of convergence (uniformly over $i$) when both the number of units $N$ and the number of pre-treatment periods $T_0$ go to infinity. Notably, if $T_0$ goes to infinity sufficiently quickly, then our estimators achieve the optimal nonparametric rate of convergence in stone1980optimal (up to a log term), as if we observed the latent $\alpha_i$.

Once the values of the conditional means and the propensity scores are imputed, we construct a doubly robust estimate of the average treatment effect on the treated (ATT) for a given post-treatment period $t$ (e.g., robins1994estimation,chernozhukov2018double). We propose a cross-fitting scheme which, together with double robustness, helps to ensure that our treatment effect estimates are $\sqrt{N}$-consistent, and asymptotically normal and unbiased. Remarkably, under certain regularity conditions and if $T_0$ grows sufficiently fast, our estimator achieves the semiparametric efficiency bound as if $\alpha_i$ was observed.

We provide extensive simulation evidence of the efficacy of our methods. We show reliable confidence interval coverage over a range of DGPs.

Related Literature

This paper contributes to the vast and rapidly growing literature on causal inference in panels. We refer the reader to arkhangelsky2024causal for an excellent recent overview of this field. DiD and TWFE are very widely used still growing in popularity in applied work, perhaps due to their practicality and perceived transparency (goldsmith2024tracking); see de2023two and roth2023s for extensive reviews of the DiD and TWFE literature. A focus of the recent econometrics literature concerning these methods is the accomodation of heterogeneous treatment effects (e.g., de2020two,callaway2021difference,sun2021estimating,wooldridge2021two,borusyak2024revisiting). However, the validity of DiD and TWFE methods crucially relies on the parallel trends assumption. A number of recent papers re-evaluate the restrictiveness of parallel trends in DiD and TWFE and in some cases propose relaxing the assumption at the expense of point identification (e.g., manski2018right,rambachan2023more,ghanem2022selection). In contrast, our framework relaxes the parallel trends assumption and allows for rich heterogeneity in treatment effects and selection mechanisms, but maintains point identification. However, in order to achieve this we require a sufficiently long pre-treatment history in contrast to these other methods.

Alternative popular approaches to causal inference in panels include SC and matrix completion methods; see abadie2021using for a recent review of SC methods. Some recent developments in the rapidly growing SC and causal matrix completion literatures include, among others, arkhangelsky2021synthetic,cattaneo2021prediction,chernozhukov2021exact,chernozhukov2018t and athey2021matrix,bai2021matrix,agarwal2023causal,abadie2024doubly, respectively. Similarly to the interactive fixed effects panel literature (e.g., pesaran2006estimation,bai2009panel), SC and matrix completion methods generalize the TWFE framework by allowing the (untreated) potential outcomes to follow a factor model. Similarly to this paper, these methods impute the missing counterfactual outcomes by leveraging the factor structure of a long history of pre-treatment outcomes. However, to establish statistical guarantees and validity of inference, most methods assume a linear factor model, i.e, $\mu^{(0)} (\alpha_i,\lambda_t) = \alpha_i' \lambda_t$, whereas we allow for a general nonparametric and nonlinear factor model. While some matrix completion methods consider extensions to nonlinear factor models, these papers either lack inferential theory and/or impose strong smoothness requirements on $\mu^{(0)}(\cdot,\cdot)$, e.g., assume that $\mu^{(0)}(\cdot,\cdot)$ belongs to a H\"older class (e.g., agarwal2020synthetic,fernandez2021low,athey2025identificationaveragetreatmenteffects). In contrast, our imputation method (nearly) achieves the optimal nonparametric rate of convergence and allows us to provide a semiparametrically efficient estimator of the ATTs under a substantially weaker Lipschitz continuity assumption.

We measure latent similarity using a pseudo-distance suggested in zhang2017estimating, allowing us to identify individuals with similar latent characteristics, despite the inherent non-identifiability of the latent factors $\alpha$. Various versions of this pseudo-distance have been utilized in the recent literature in numerous applications, including non-parametric graphon estimation zhang2017estimating,zeleneev2020identification,nowakowicz2024nonparametric, controlling for unobservables using network data auerbach2022identification, and estimation of nonlinear factor models feng2024optimal. Other recent applications to estimation of treatment effects in network and panel models include wang2022linking, hoshino2024estimating, and athey2025identificationaveragetreatmenteffects but, unlike this paper, none of these works provides inferential theory.

Perhaps most closely related to the present work are feng2024causal and abadie2024doubly. feng2024causal studies causal inference in a cross-sectional setting. He uses a large number of auxiliary variables following a nonlinear factor model to control for unobservables. These auxiliary variables play effectively the same role as the pre-treatment history in our setting. Similarly to our work, feng2024causal imputes the counterfactual means and the propensity scores and then uses the imputed values to construct doubly robust estimators and confidence intervals for causal estimands of interest. However, his imputation method is different from ours. Specifically, feng2024causal combines the pseudo-distance of zhang2017estimating with local principal component analysis (PCA) to estimate latent factors and loadings and then employs quasi-maximum likelihood using these estimates for imputation. On the other hand, our approach imputes the counterfactual means and propensity score directly without estimating the factors nor factor loadings. The simplicity of our approach results in both theoretical and practical advantages. First, we are able to establish consistency of our estimators and obtain the desired rates of convergence under weaker smoothness requirements. Second, we find that our approach also performs better than feng2024causal's method in small and moderate sample sizes typical for microeconometric applications.

abadie2024doubly is another recent paper that provides a doubly robust approach to inference in latent factor models. While the standard matrix completion methods exclusively focus on imputing the counterfactual means, abadie2024doubly suggest applying these methods to estimate the matrix of propensity scores as well. The latter imputation approach is fundamentally different than the one proposed in this paper. Specifically, abadie2024doubly rely on denoising a large matrix of treatment assignments which have an underlying low-rank structure, whereas our approach imputes the propensity scores based on the pre-treatment outcomes. Another important difference is that abadie2024doubly consider linear factor models, whereas we focus on and provide formal statistical guarantees for nonlinear factor models while imposing minimal smoothness requirements.

Our estimation approach is based on a doubly robust estimate of the ATT. As shown in stone1980optimal, the optimal rate of convergence for nonparametric estimators is slow under weak smoothness assumptions, which can complicate inference. The use of doubly robust/Neyman orthogonal estimation (e.g., robins1994estimation, robins1995semiparametric, hahn1998role, scharfstein1999adjusting,chernozhukov2018double,chernozhukov2022debiased) can help ensure centered asymptotic normality for low-dimensional estimands in the presence of nonparametric nuisance parameters. Some recent literature on panel data leverages double robustness to achieve valid inference under relatively weak conditions (e.g., arkhangelsky2022doubly,arkhangelsky2024double,sant2020doubly).

The rest of the paper is organized as follows. Section (ref) introduces the framework and describes our estimator. Section (ref) provides formal statistical guarantees. Section (ref) presents numerical evidence. All auxiliary lemmas and proofs are provided in Appendix (ref).

comment\subsection{Previous intro and literature review} \subsubsection{Introduction from previous draft} Treatment effects estimation is important for policy evaluation. From the workhorse difference-in-differences methods with two-way fixed effects approach to recently developed event study methods, researchers are able to identify and estimate the causal effects with large panels, e.g., borusyak2021revisiting, callaway2021difference, sun2021estimating, wooldridge2021two, etc. These methods often rely on the assumption of unconditional or conditional parallel trends for outcomes. Consider the potential outcome framework with two treatment statuses: $Y_{i,t}(1)$ and $Y_{i,t}(0)$ are the counterfactual outcomes for individual $i$ at time $t$ if it gets treated or not. The unconditional parallel trend assumption requires that the outcome differences between different time periods should stay the same across different treatment statuses, that is, \begin{align*} \mathbb{E}[Y_{i,t}(0) - Y_{i,\ell}(0) \vert i in treatment group] = \mathbb{E}[Y_{i,t}(0) - Y_{i,\ell}(0) \vert i in control group] \end{align*} Despite these assumptions providing a more accessible identification and estimation strategy, these assumptions are invalid in a lot of scenarios with unobserved heterogeneity. Therefore, policy evaluation is even more challenging with complex treatment assignment mechanisms. Consider the following model with unobserved latent factors: \begin{align*} Y_{i,t}(0) &= \mu^{(0)}(\alpha_i, \lambda_t) + u_{i,t} \\ w_{i,t} &= p(\alpha_i, \lambda_t) + \epsilon_{i,t} \end{align*} where $Y_{i,t}(0)$ and $Y_{i,t}(1)$ are the outcome for individual $i$ at time $t$ for two treatment statuses, $\mu^{(0)}$ and $\mu^{(1)}$ are the unknown and unobserved functions, $\alpha_i$ and $\lambda_t$ are the unobserved individual and time latent factors, $w_{i,t}$ is the treatment status of individual $i$ at time $t$ that depends on $\alpha_i$ and $\lambda_t$, and $u_{i,t}, v_{i,t}, \epsilon_{i,t}$ are the unobserved error terms with mean zero conditional on the latent factors. Under the selection on unobserved latent factors setting, matrix completion methods tackle this model by imputing counterfactual outcomes on different treatment statuses as they are built with the factor model in mind and generally allow different missing mechanisms. And the matrix completion methods for causal inference, such as agarwal2021causal, athey2021matrix, bai2021matrix, fernandez2021low, require the smoothness assumption on the underlying function and the panel (or matrix) that it generated. Commonly used assumptions in the matrix completion literature allow that the counterfactual outcomes can be approximated by a low-rank matrix. The first set of assumptions is imposed on the panel or matrix itself, such as that the matrix needs to be generated with a specific factor structure, e.g., $\mu^{(0)}(\alpha_i, \lambda_t) = \alpha_i + \lambda_t$ or $\mu^{(0)}(\alpha_i, \lambda_t) = \alpha_i^{\intercal} \lambda_t$. This leads to the low-rank assumption that is heavily used, that is, $\mbox{rank}\left( \left[ \mu^{(0)}(\alpha_i, \lambda_t) \right]_{i\in\{1,...,N\}, t\in\{1,...,T\}} \right)$ is low. Another set of assumptions is the smoothness condition directly imposed on the function $g$, e.g., $g$ is infinitely differentiable or has fast singular value decay. However, these low-rank or smoothness assumptions that allow the matrix to be well approximated by the finite number of factors are strict and may not be satisfied in complex data generating processes. Hence, we aim to relax the assumption and impose a weaker smoothness condition. The proposed method imputes unobserved potential outcomes through matrix completion and accesses treatment effects estimates with weaker smoothness assumptions. Following zhang2017estimating, our approach is to, first, construct the identifiable and estimable pseudo distance measurement on observable to allow matching individuals with similar unobserved latent characteristics. Then we estimate counterfactuals with the corresponding neighborhood for each treated individual. Finally, we estimate the $\mathrm{ATT}_t$ with the previously acquired counterfactual estimations. The research contributes to the literature on causal inference for treatment evaluation and matrix completion. The proposed method is built on a general model specification with latent factors without assuming the low-rank structure and the reliance on the low-rank approximation, which is in contrast to most of the existing matrix completion literature, for the nonparametric estimation. It also allows the selection of latent factors and flexible missing mechanisms. We establish a consistency result with an acceptable rate of convergence for the treatment effect estimators. \subsubsection{Summary of other literature} \begin{itemize} • Proposed method estimates treatment effects with weaker smoothness assumption (Lipschitz condition) \begin{itemize} • Parallel trend assumption used in difference-in-differences and event study is not required: ben2021synthetic; borusyak2021revisiting; callaway2021difference; sun2021estimating • Strong smoothness assumption used in synthetic control and matrix completion is not required\\ abadie2010synthetic; agarwal2021causal; arkhangelsky2021synthetic; athey2021matrix; bai2021matrix; fernandez2021low \end{itemize} • Asymptotic inference through doubly robust AIPW estimation with latent factors \begin{itemize} • stone1980optimal shows that the optimal rate of convergence for nonparametric estimator is slow under weak smoothness assumption, results in unaccessible inference theorem. • Double robustness of AIPW estimator and double/debiased machine learning\\ abadie2024doubly; chernozhukov2018double; chernozhukov2022debiased; hahn1998role; robins1994estimation; robins1995semiparametric; scharfstein1999adjusting; singh2024double; smucler2019unifying • As mentioned in smucler2019unifying, there are two approaches to achieve rate double robustness: either imposing strong Donsker/smoothness assumptions on both nuisance functions as in robins2008higher, or performing sample splitting to regain independence as in chernozhukov2018double and abadie2024doubly. \end{itemize} • Pseudo distance related literature \begin{itemize} • zhang2017estimating; zeleneev2020identification • feng2024causal studies nonlinear factor model by utilizes zhang2017estimating's distance measurement to construct local neighborhood with similar latent characteristics. With such neighborhood, they perform local PCA to estimate latent factors and loadings for outcomes and propensity scores imputation. • hoshino2024estimating uses zhang2017estimating's distance measurement to estimate both the outcome model and propensity score with kernel smoothing, but their context is on dyadic data. Also, although they estimated the propensity score as we did, they didn't provide the inference method with a confidence interval. Alternatively, they provide a permutation test for size control. • wang2022linking study the effect of relationships (linking effect or peer effect) in the network context and utilizes zhang2017estimating's distance measurement to estimate the propensity score where the treatment is the assignment of links between individual. \end{itemize} \end{itemize}

Notation and Proposed Method

Model and Notation

Our sample consists of individuals indexed by $i=1,...,n$ and time periods indexed by $t=1,...,T$. For each individual $i$ and period $t$, we observe an outcome $Y_{i,t}$ and a binary indicator $w_{i,t}$. The indicator $w_{i,t}$ is equal to one if individual $i$ is treated at or prior to time $t$ and zero otherwise. We assume throughout that no individual in the population is treated prior to a period $T_0+1$ with $T_0<T$. We let $Y_{i,s:t}$ denote the vector of outcomes for individual $i$ from periods $s$ to $t$ inclusive and similarly for other variables.

Let $Y_{i,t}(0)$ denote the potential outcome of individual $i$ at time $t$ under a counterfactual in which the individual has not yet received treatment at period $t$. We assume that if $w_{i,t}=0$ then $Y_{i,t}=Y_{i,t}(0)$. This precludes the possibility that individuals anticipate future treatment and that this impacts their outcomes (e.g., Abbring2003, borusyak2021revisiting, sun2021estimating). We state this formally in Assumption (ref) below.

assu{1}[No Anticipation] If $w_{i,t}=0$ then $Y_{i,t}=Y_{i,t}(0)$.

Central to our analysis is the assumption of a non-linear factor model for the untreated potential outcomes. Let $\alpha_i$ be some latent and time-invariant characteristics of individual $i$, and let $\lambda_t$ be period-specific factors. $\alpha_i$ may capture say, unobserved demographic characteristics, individual preferences, or innate ability. $\lambda_t$ may capture e.g., unobserved macro-economic conditions, government policy, or environmental factors. We model the potential outcomes and realized treatments as follows.

align[align omitted — 258 chars of source]

The residuals $u_{i,t}$ and $\epsilon_{i,t}$ are unobserved. The functions $\mu^{(0)}$ and $p_t$ are unknown and we do not assume that they have any particular functional form. Note that the model above implies that the mean of $Y_{i,t}(0)$ given the latent factors is time-invariant and does not depend on the individual $i$. However, the time-subscript on $p_t$ allows the conditional mean of $w_{i,t}$ to vary over time, which reflects the fact that this variable is increasing and that it is therefore non-stationary.

Knowledge of the latent factors $\alpha_i$ and $\lambda_t$ is insufficient to identify causal quantities of interest. This is because the residual in the treatment model ((ref)) may be correlated with the residual in the outcome equation ((ref)). In particular, there may be confounding factors that are not included in $\alpha_i$ and $\lambda_t$ and which influence both the outcome and treatment status. In order to achieve identification we make the key assumption that the latent factors $\alpha_i$ and $\lambda_t$ together account for all confounding between the outcome and treatment. Formally, we assume that after controlling for these latent factors, there is no residual dependence between the untreated potential outcome and treatment.

assu{2}[Latent Unconfoundedness] $Y_{i,t}(0)\perp \!\!\! \perp w_{i,t}|\alpha_i,\lambda_t$.

The assumption of latent unconfoundedness is common in the literature (e.g. abadie2024doubly, agarwal2021causal, arkhangelsky2022doubly, athey2021matrix, fernandez2021low). Note that the assumption effectively requires that the latent factors are sufficiently rich.

Exposure to macroeconomic conditions and other shared aggregate factors may induce dependence between the potential outcomes of different individuals. We assume that the time-specific factors $\lambda_t$ are sufficiently rich that after controlling for these factors and the individual-specific factors, the untreated potential outcomes of any two individuals are independent.

assu{3}[Latent Independence] For any individuals $i\neq j$, $Y_{i,t}(0)\perp \!\!\! \perp Y_{j,t}(0)|\{\alpha_k\}_{k=1}^N,\lambda_t$.

Under Assumption (ref), key counterfactual quantities of interest can be written in terms of the latent factors and functions $\mu^{(0)}$ and $p$. For notational convenience, let $p_{i,t}:=p_t(\alpha_i) $ and $\mu^{(0)}_{i,t}:=\mu^{(0)}(\alpha_i, \lambda_t)$. Under Assumption (ref), the average effect of treatment on the treated at time $t$, which is defined as $\mathbb{E}_t[Y_{i,t}-Y_{i,t}(0)|w_{i,t}=1]$, can be written in doubly-robust form as follows. \[ \mathrm{ATT}_t=\frac{1}{P(w_{i,t}=1)}\mathbb{E}_t\bigg[Y_{i,t}w_{i,t}-\frac{(1-w_{i,t})Y_{i,t}p_{i,t}+(w_{i,t}-p_{i,t})\mu^{(0)}_{i,t}}{1-p_{i,t}}\bigg] \] The time subscript on the expectation above indicates that it is taken with respect to the period $t$-specific distribution of the observables (i.e., conditional on $\lambda_t$). The doubly-robust form of ATT above is proposed and used in panel data settings by sant2020doubly. Additionally, the proposed method can also adopt other causal estimands with their respective doubly robust score functions (e.g., arkhangelsky2022doubly).

In this paper we provide new methods for estimating $p_{i,t}$ and $\mu^{(0)}_{i,t}$. Given these estimates, one can construct a doubly-robust estimate of $\mathrm{ATT}_t$ and perform inference on this object. Let $\hat{p}_{i,t}$ and $\hat{\mu}^{(0)}_{i,t}$ be estimates of $p_{i,t}$ and $\mu^{(0)}_{i,t}$ respectively. Then a corresponding doubly-robust estimate of $\mathrm{ATT}_t$ is given below, where $N$ is the sample size, and $N_{1,t}$ the number of individuals treated by time $t$.

equation[equation omitted — 215 chars of source]

comment\begin{itemize}[nosep] • $Y_{i,t}(0)$ and $Y_{i,t}(1)$: Counterfactual outcomes for individual $i$ and time $t$ • $w_{i,t}$: Indicator if individual $i$ is treated at time $t$ • $\alpha_i$ and $\lambda_t$: Unobserved individual and time latent factors • $\mu^{(0)}$, $\mu^{(1)}$, $p$: Unknown functions, where \begin{itemize}[nosep] • $\mu^{(0)}_{i,t} := \mu^{(0)}(\alpha_i, \lambda_t)$ and $\mu^{(1)}_{i,t} := \mu^{(1)}(\alpha_i, \lambda_t)$ are the expected counterfactual outcomes • $p(\alpha_i, \lambda_t)$ is the propensity score, $p_i(\alpha_i, \lambda_t) = p_{i,t} = \mathbb{E}\{ \mathbf{1}\{w_{i,t} = 1\} | \alpha_i, \lambda_t \} = \mathbb{E}\{ w_{i,t} | \alpha_i, \lambda_t \}$ \end{itemize} • $u_{i,t}$, $v_{i,t}$, $\epsilon_{i,t}$: Unobserved error terms with $\mathbb{E}\left(u_{i,t}\vert\alpha, \lambda\right)=\mathbb{E}\left(v_{i,t}\vert\alpha, \lambda\right)=\mathbb{E}\left(\epsilon_{i,t}\vert\alpha, \lambda\right)=0$ \end{itemize} Under the potential outcome framework, researchers are only be able to observe one of the possible intervention. For example, for individuals in the treatment group, we are not able to observe their counterfactual outcomes if they happen to be in the control group. In order to impute and estimate the counterfactual outcomes, we will need to rely on other auxilary information: \begin{align*} Z_{i,t} = g(\alpha_i, \lambda_t) + \eta_{i,t} \end{align*} where both dimensions of $i$ and $t$ are large. In the panel data setting, the pre-treatment period for all groups are used as such auxilary information. \begin{align*} Z_{i,t} = g(\alpha_i, \lambda_t) + \eta_{i,t} \Rightarrow Y_{i,t}(0) &= \mu^{(0)}(\alpha_i, \lambda_t) + u_{i,t} \end{align*} Here, for simplicity, we only consider simultaneous treatment assignment mechanism, where \begin{itemize}[nosep] • $N$: Number of total individuals • $T$: Number of total time periods \begin{itemize}[nosep] • $T_0$: Number of pre-treatment time periods \end{itemize} \end{itemize} \subsection{Estimands and Estimators} \paragraph{Imputation for outcomes and propensity scores} \begin{itemize}[nosep] • $\hat{\mu}_{i,t}^{(0)}$ and $\hat{\mu}_{i,t}^{(1)}$: Imputation for counterfactual outcomes • $\hat{p}_{i,t}$: Estimator for propensity score \begin{itemize}[nosep] • $\hat{p}_t$: Estimator for probability of getting treated, $\hat{p}_t = \frac{1}{N}\sum\limits_{i=1}^N \frac{w_{i,t}}{N}$ \end{itemize} \end{itemize} \paragraph{Doubly robust exptected outcome estimation} \begin{itemize}[nosep] • Conditional expected outcome: \begin{align*} \mathbb{E}[Y_{i,t}(0) | w_{i,t} = 1] = \mathbb{E} \left\{\frac{w_{i,t} \cdot \mu_{i,t}^{(0)}}{p_t} + \frac{p_{i,t}}{1-p_{i,t}} \frac{(1-w_{i,t}) \cdot [Y_{i,t} - \mu_{i,t}^{(0)}]}{p_t} \right\} \end{align*} • Doubly robust imputation estimator: \begin{align*} \hat{\mathbb{E}}[Y_{i,t}(0) | w_{i,t} = 1] = \frac{1}{N} \sum\limits_{i=1}^{N} \left\{\frac{w_{i,t} \cdot \hat{\mu}_{i,t}^{(0)}}{\hat{p}_t} + \frac{\hat{p}_{i,t}}{1-\hat{p}_{i,t}} \frac{(1-w_{i,t}) \cdot [Y_{i,t} - \hat{\mu}_{i,t}^{(0)}]}{\hat{p}_t} \right\} \end{align*} \end{itemize} \paragraph{Doubly robust ATE estimation} \begin{itemize}[nosep] • $\mbox{ATE}_t$ and $\mbox{ATE}$: \begin{align*} ATE_t &= \mathbb{E}[ Y_{i,t}(1) - Y_{i,t}(0) ] = \mathbb{E} \left\{ \mu_{i,t}^{(1)} - \mu_{i,t}^{(0)} + w_{i,t} \frac{Y_{i,t}-\mu_{i,t}^{(1)}}{p_{i,t}} - (1-w_{i,t}) \frac{Y_{i,t}-\mu_{i,t}^{(0)}}{1-p_{i,t}} \right\} \\ ATE &= \mathbb{E}[ATE_t] \end{align*} • Doubly robust $\mbox{ATE}_t$ and $\mbox{ATE}$ estimators: \begin{align*} \hat{ATE}_t &= \frac{1}{N}\sum_{i=1}^{N} \left\{ \hat{\mu}_{i,t}^{(1)} - \hat{\mu}_{i,t}^{(0)} + w_{i,t} \frac{Y_{i,t}-\hat{\mu}_{i,t}^{(1)}}{\hat{p}_{i,t}} - (1-w_{i,t}) \frac{Y_{i,t}-\hat{\mu}_{i,t}^{(0)}}{1-\hat{p}_{i,t}} \right\} \\ \hat{ATE} &= \frac{1}{T-T_0}\sum_{i=T_0}^{T} \hat{ATE}_t \end{align*} \end{itemize}

Proposed Method

Our proposed estimation method uses pre-treatment outcomes to find untreated individuals whose latent factors $\alpha_i$ are similar to those of treated individuals. To motivate our approach, let us first suppose $\alpha_i$ were observed. Under Assumptions (ref) and (ref), we have

equation*[equation* omitted — 166 chars of source]

The objects on the right-hand sides of each equation above are regression functions. If $\alpha_{i,t}$ were observable, we could apply Nadaraya-Watson to non-parametrically estimate the functions $\mu^{(0)}(\cdot,\lambda_t)$ and $p_t(\cdot)$ and thus $p_{i,t}$ and $\mu^{(0)}_{i,t}$. To be precise, we could obtain the following estimates.

align*[align* omitted — 296 chars of source]

The estimates above are infeasible because in practice, we do not directly observe $\alpha_i$ for any individual $i$. In order to obtain feasible estimates, we replace the infeasible distance $\|\alpha_i-\alpha_j\|$ in the expressions above, with a feasible pseudo-distance. To define this pseudo-distance, suppose that no individuals in the population are treated prior to some period $T_0+1$ and let $\langle\cdot,\cdot \rangle$ be the Euclidean inner-product. We use the pseudo-distance defined below:

align*[align* omitted — 141 chars of source]

A pseudo-distance of the form above is employed in zhang2017estimating. To motivate the use of the pseudo-distance, suppose Assumption (ref) holds. Then conditional on the individual latent factors, the inner-product in the pseudo-distance is an unbiased estimate of the object on the right-hand side below. \[ \mathbb{E}[\langle Y_{k,1:T_0}, Y_{i,1:T_0} - Y_{j,1:T_0} \rangle|\alpha_i,\alpha_j,\alpha_k]=\int \mu^{(0)}(\alpha, \lambda) \big( \mu^{(0)}(\alpha_i, \lambda) - \mu^{(0)}(\alpha_j, \lambda) \big) \dd \pi(\lambda). \] In the above, $\pi$ is the stationary distribution of $\lambda_t$. As such, we can understand the pseudo-distance as the sample analogue of the population pseudo-distance below, where $\mathcal{A}$ is the support of $\alpha_i$.

align*[align* omitted — 313 chars of source]

The inequality relates the size of the pseudo-distance $d_{i,j}$ to a squared $L_2$ distance between the functions $\mu^{(0)}(\alpha_i, \cdot)$ and $\mu^{(0)}(\alpha_j, \cdot)$. Thus $d_{i,j}$ measures similarity of the latent factors to the extent that they impact outcomes. We provide sufficient conditions for the consistency of the sample pseudo-metric to this quantity.

It is worth contrasting the sample pseudo-distance above with the Euclidean distance between the history of pre-treatment outcomes, which is defined as $\|Y_{i,1:T_0} - Y_{j,1:T_0}\|^2$. The mean of the squared Euclidean distance conditional on the individual latent factors is given below, where we again assume that untreated potential outcomes of different individuals are independent conditional on the latent factors.

align*[align* omitted — 275 chars of source]

Thus the Euclidean distance is increasing in the conditional residual variance $E[u_{j,t}^2|\alpha_j,\lambda_t=\lambda]$. Thus in the presence of conditional heteroskedasticity, a small Euclidean distance may reflect that individual $j$'s outcomes have a low residual variance and not that individual $j$ and $i$ have similar latent factors. It is the need to be robust to conditional heteroskedasticity that motivates our use of the pseudo-metric.

Given the pseudo-distance, we may form feasible estimates of $p_{i,t}$ and $\mu^{(0)}_{i,t}$ as follows.

align*[align* omitted — 264 chars of source]

One could plug the estimates above into the the formula for the doubly-robust $\mathrm{ATT}_t$ estimate ((ref)). However, we instead employ a cross-fitting scheme to further de-bias our estimates. This is in-line with the extensive literature on double machine learning, which demonstrates the utility of cross-fitting for reducing bias and obtaining valid inference (e.g., chernozhukov2018double, abadie2024doubly). The full algorithm with cross-fitting is detailed below along with a variance estimate and confidence interval.

algorithm[algorithm omitted — 3,085 chars of source]

comment\begin{algorithm}[H] \caption{Doubly-Robust Estimation and Inference with Cross-Fitting} Inputs: Number of folds $J$, bandwidth $h$, confidence level $\alpha$. Returns: $\mathrm{ATT}_t$ estimate $\hat{\mathrm{ATT}}_t$, variance estimate $\hat{V}_t$, and level $1-\alpha$ confidence interval. \begin{algorithmic}[1] \STATE Randomly partition $[N] = \{1, 2, ..., N\}$ into $J$ folds $\{\mathcal{I}_j\}_{j=1}^J$ of size $\approx N/J$. Let $\mathcal{I}_{-j}=[N]\setminus\mathcal{I}_j$. \FORALL{ $j\in[J]$ and $i \in \mathcal{I}_{j}$} \STATE Calculate $\hat{\mu}_{i,t}^{(0)}$ and $\hat{p}_{i,t}$ as follows. \begin{align*} \hat{\mu}_{i,t}^{(0)} &\leftarrow \frac{\sum\limits_{j\in\mathcal{I}_{-j}; w_{j,t} = 0} K(\hat{d}_{i,j}/h)Y_{j,t}}{\sum\limits_{j\in\mathcal{I}_{-j}; w_{j,t} = 0}K(\hat{d}_{i,j}/h)}, \hat{p}_{i,t} \leftarrow \frac{\sum\limits_{j\in\mathcal{I}_{-j}} K(\hat{d}_{i,j}/h) w_{j,t}}{\sum\limits_{j\in\mathcal{I}_{-j}}K(\hat{d}_{i,j}/h)} \end{align*} \ENDFOR \STATE Construct $\hat{\mathrm{ATT}}_t$ using the formula below. \begin{align*} \hat{\mathrm{ATT}}_t&\leftarrow\frac{1}{N_{1,t}}\sum_{i=1}^N\bigg(Y_{i,t}w_{i,t}-\frac{(1-w_{i,t})Y_{i,t}\hat{p}_{i,t}+(w_{i,t}-\hat{p}_{i,t})\hat{\mu}^{(0)}_{i,t}}{1-\hat{p}_{i,t}}\bigg)\\ \hat{V}_t&\leftarrow\frac{N}{N_{1,t}^2}\sum_{i=1}^N\bigg(Y_{i,t}w_{i,t}-\frac{(1-w_{i,t})Y_{i,t}\hat{p}_{i,t}+(w_{i,t}-\hat{p}_{i,t})\hat{\mu}^{(0)}_{i,t}}{1-\hat{p}_{i,t}}-\frac{N_{1,t}}{N}\hat{\mathrm{ATT}}_t\bigg)^2 \end{align*} \STATE Form $1-\alpha$ level confidence interval $\big[\hat{\mathrm{ATT}}_t \pm Z_{(1-\alpha)/2}\sqrt{\hat{V}/N}\big]$. \end{algorithmic} \end{algorithm}

commentConsider the simplifed setting with only one post-treatment time period $T$. To tackle the nonidentifiable $\ell_2$ ditance $\Vert\alpha_i - \alpha_j\Vert_2$ with unobserved latent factor $\alpha_i$ and $\alpha_j$, we utilize the following pseudo distance: \begin{align*} d^2(i, j) = \sup_{\alpha_k \in \operatorname*{supp}(\alpha)} \left\vert \int_{\lambda \in \operatorname*{supp}(\lambda)} \mu^{(0)}(\alpha_k, \lambda) \left[ \mu^{(0)}(\alpha_i, \lambda) - \mu^{(0)}(\alpha_j, \lambda) \right] \dd \pi(\lambda) \right\vert \end{align*} where $\pi$ is the probability measure, and such pseudo distance can be estimated by \begin{align*} \hat{d}^2(i, j) = \frac{1}{T_0} \max\limits_{k \notin \{i, j\}} \vert \langle Y_{k,\cdot}, Y_{i,\cdot} - Y_{j,\cdot} \rangle \vert \end{align*} With the estimated pseudo distance, for each individual $i$ at post-treatment time period $T$, we impute the counterfactual outcomes and propensity scores \begin{align*} \hat{\mu}_{i,T}^{(w)} &= \frac{\dfrac{1}{Nh}\sum\limits_{j; w_{j,T} = w} K\left(\dfrac{\hat{d}_{i,j}}{h}\right)Y_{j,T}(w)}{\dfrac{1}{Nh}\sum\limits_{j; w_{j,T} = w}K\left(\dfrac{\hat{d}_{i,j}}{h}\right)} \hat{p}_{i,T} = \frac{\dfrac{1}{Nh}\sum\limits_{j} K\left(\dfrac{\hat{d}_{i,j}}{h}\right) w_{j,T}}{\dfrac{1}{Nh}\sum\limits_{j}K\left(\dfrac{\hat{d}_{i,j}}{h}\right)} \end{align*} where $h$ is the bandwidth for the kernel method. Lastly, we plug the imputations into the doubly robust AIPW score function and get the estimator: \begin{align*} \hat{ATE} = \frac{1}{N} \sum_{i=1}^{N} \left[ \hat{\mu}_{i,T}^{(1)} - \hat{\mu}_{i,T}^{(0)} + w_{i,T} \frac{Y_{i,T} - \hat{\mu}_{i,T}^{(1)}}{\hat{p}_{i,T}} - (1-w_{i,T}) \frac{Y_{i,T} - \hat{\mu}_{i,T}^{(0)}}{1-\hat{p}_{i,T}} \right] \end{align*}

comment\subsection{2-Fold Cross-Fitting} \subsubsection*{Algorithm (following previous simplified setting with one post-treatment period)} Input: $Y_{i,t}$, $W_i$ where \begin{small}$\begin{array}{l}i = 1,..., N \\ t = 1, ..., T_0, T\end{array}$\end{small} Output: Estimator and C.I. for $\mbox{ATE}$\\ \begin{itemize}[nosep] • Randomly partition $[N] = \{1, 2, ..., N\}$ data into two folds: $I = \{1, 2, ..., K\}$ and $I^c = \{K+1, K+2, ..., N\}$, where both partitions have sufficient amount of treated and control observations. • For every $i \in I$, calculate the distance to $j \in I^c$ by \begin{align*} \hat{d}^2(i, j) &= \frac{1}{T_0} \max\limits_{k \in I_2 \;|\; k \neq i, j} \vert \langle Y_{k,\cdot}, Y_{i,\cdot} - Y_{j,\cdot} \rangle \vert \end{align*} and use them to aquire $\hat{\mu}_{i,T}^{(0)}, \hat{\mu}_{i,T}^{(1)}, \hat{p}_{i,T}$ for all $i \in I$ by the proposed method. • Swap the roll of $I$ and $I^c$ and perform step 2 again to aquire $\hat{\mu}_{i,t}^{(0)}, \hat{\mu}_{i,t}^{(1)}, \hat{p}_{i,t}$ for all $i \in I^c$. • Construct the estimators for $\mbox{ATE}$ and $V$ \begin{align*} \hat{ATE} &= \frac{1}{N} \sum_{i=1}^{N} \left[ \hat{\mu}_{i,t}^{(1)} - \hat{\mu}_{i,t}^{(0)} + w_{i,T} \frac{Y_{i,T} - \hat{\mu}_{i,t}^{(1)}}{\hat{p}_{i,T}} - (1-w_{i,T}) \frac{Y_{i,T} - \hat{\mu}_{i,t}^{(0)}}{1-\hat{p}_{i,T}} \right] \\ \hat{V} &= \frac{1}{N-1} \sum_{i=1}^{N} \left\{ \left[ \hat{\mu}_{i,t}^{(1)} - \hat{\mu}_{i,t}^{(0)} + w_{i,T} \frac{Y_{i,T} - \hat{\mu}_{i,t}^{(1)}}{\hat{p}_{i,T}} - (1-w_{i,T}) \frac{Y_{i,T} - \hat{\mu}_{i,t}^{(0)}}{1-\hat{p}_{i,T}} \right] - \hat{ATE} \right\}^2 \end{align*} • The 95% confidence interval: $\left(\hat{\mbox{ATE}} \pm Z_{0.05/2} \frac{1}{\sqrt{N}} \sqrt{\hat{V}} \right)$ \end{itemize}

Large Sample Theory

The estimator $\hat{\mathrm{ATT}}_t$ is doubly robust and employs cross-fitting. An extensive literature (e.g., robins2008higher, chernozhukov2018double, abadie2024doubly) provides sufficient conditions for $\sqrt{N}$-consistency of such estimates and asymptotically correct coverage of the corresponding confidence intervals. A key condition is that the first stage nuisance-parameter estimates (in our case $\hat{\mu}^{(0)}_{i,t}$ and $\hat{p}_{i,t}$) converge sufficiently quickly. For this condition, the following rates suffice: \[ \max_{i \in [N]} \left\vert \hat{\mu}_{i,t}^{(0)} - \mu_{i,t}^{(0)} \right\vert = o_p(N^{-1/4})\hspace{30pt}\text{and}\hspace{30pt} \max_{i \in [N]} \left\vert \hat{p}_{i,t}- p_{i,t} \right\vert =o_p(N^{-1/4}) \]

We establish convergence rates for the first-stage estimates. These rates may be of interest per se because $\mu_{i,t}^{(0)}$ is the optimal prediction of individual $i$'s untreated potential outcome given the latent factors. We show that under certain conditions, the estimates achieve the stone1980optimal optimal rate for non-parametric regression on $\alpha_i$ under Lipschitz continuity. That is, we can achieve the same optimal rate attainable for the infeasible estimates that take $\alpha_i$ as known. However, in order to achieve this rate, we require that the number of pre-treatment periods $T_0$ grows sufficiently quickly with $N$.

comment\begin{align*} \left\Vert\hat{\mu}^{(0)} - \mu^{(0)}\right\Vert = \mathcal{O}_p\left(R_{\mu}\right) and \left\Vert\hat{p} - p\right\Vert = \mathcal{O}_p\left(R_p\right) \end{align*} stone1980optimal. The centric difficult inference procedure of nonparametric estimators is due to failing to achieve the root-$N$ rate of convergence, unlike other common parametric estimators. The augmented inverse propensity weighted (AIPW) estimator offers a valid approach for inference-making with nonparametric estimators, primarily due to its rate double robustness property related to convergence rates. Although the convergence rates for both the outcome model and the propensity score model are slower than the root-$N$ rate, the inference process remains accessible under certain regularity conditions if the product of their rates is sufficiently fast and attains the root-$N$ rate. That is, suppose the rates of convergence for the outcome and the propensity score model are denoted as \begin{align*} \left\Vert\hat{\mu}^{(0)} - \mu^{(0)}\right\Vert = \mathcal{O}_p\left(R_{\mu}\right) and \left\Vert\hat{p} - p\right\Vert = \mathcal{O}_p\left(R_p\right) \end{align*} Although both $R_{\mu}$ and $R_p$ are slower than root-$N$ when $\mu^{(w)}$ and $p$ are being nonparametrically estimated, the doubly robust estimator for the target causal estimand can still achieve asymptotic normality \begin{align*} \sqrt{N}\left(\rm{DR\;Estimator} - \rm{Causal\;Estimand}\right) \rightarrow N \left(0, \mathrm{AsyVar}\right) \end{align*} whenever $R_{\mu} = o(1)$, $R_p = o(1)$, and $R_{\mu} R_p = o(N^{-1/2})$. This concept, often referred to as rate double robustness in the literature, is discussed in works such as smucler2019unifying and utilized in studies like robins2008higher, chernozhukov2018double, abadie2024doubly. In Section (ref), we will formally illustrate the assumptions and the rate of convergence for the proposed method, which has been shown to achieve double robustness. Then, in Section (ref), we will utilize this rate of convergence to establish the asymptotic inference result, along with the corresponding doubly robust estimator.

Rate of Convergence for Outcome and Propensity Score

In order to derive convergence rates for our first-stage estimates of $\mu^{(0)}$ and $p_t$ we impose the following additional assumptions.

assu{4}[Model and Latent Factors] \\ \begin{enumerate}[label=(\roman*), ref=(\roman*)] • There is some $\underline{c}>0$ so that $p_{i,t} \geq \underline{c}$ almost surely for all $i$ and $t>T_0$. • $\mathbb{E}[u_{i,t} \vert \alpha, \lambda] = 0$, and for some $\delta > 0$, $\mathbb{E}[\exp(\eta u_{i,t}) | \alpha, \lambda] \leq \exp(\delta \eta^2)$ for all $\eta \in \mathbb{R}$ almost surely and likewise for $\epsilon_{i,t}$. These error terms are jointly independent across time and across individuals conditional on the latent factors. • $\operatorname*{supp}(\alpha) \subseteq \mathcal{A}$ where $\mathcal{A}$ is a compact subset of $\mathbb{R}^{d_\alpha}$. $\alpha_i$ and $\lambda_t$ are jointly i.i.d. across $i$ and $t$, respectively. There are constants $0<\underline{c}<\bar{c}<\infty $ so that for any fixed $\alpha\in\operatorname*{supp}(\alpha)$ and any $\eta$, $\underline{c}\eta^{d_\alpha}\leq\mathbb{P}(\|\alpha-\alpha_{i}\|\leq\eta)\leq\bar{c}\eta^{d_\alpha}$. • The functions $\mu^{(0)}, \mu^{(1)}, p: \mathcal{A} \times \mathcal{L} \rightarrow \mathbb{R}$ are uniformly bounded and satisfy, for some $ L_0, L_p <\infty$, \begin{align*} \vert \mu^{(0)}(\alpha_1, \lambda) - \mu^{(0)}(\alpha_2, \lambda) \vert &\leq L_0 \Vert \alpha_1 - \alpha_2 \Vert \\ \vert p_t(\alpha_1) - p_t(\alpha_2) \vert &\leq L_p \Vert \alpha_i - \alpha_j \Vert \end{align*} for all $\alpha_i, \alpha_j \in \mathcal{A}$, $\lambda \in \mathcal{L}$, and $t\in[T]$. • $K: \mathbb{R}_{+} \rightarrow \mathbb{R}$ is weakly positive and strictly positive at zero, supported on compact set and bounded by $\bar{K} < \infty$. $K$ satisfies $\vert K(z) - K(z') \vert \leq \bar{K}' \vert z - z' \vert$ for all $z, z' \in \mathbb{R}_{+}$ for some $\bar{K}' > 0$. • There exist constants $c_1,c_2>0$ so that for all $\alpha_1, \alpha_2 \in \mathcal{A}$ \begin{align} \int_{\lambda \in \operatorname*{supp}(\lambda)} \big( \mu^{(0)}(\alpha_1, \lambda) - \mu^{(0)}(\alpha_2, \lambda) \big)^2\dd \pi(\lambda) \geq c_1\|\alpha_1-\alpha_2\|^2 \end{align} and \begin{align} & \sup_{\alpha\in\operatorname*{supp}(\alpha)}\int_{\lambda \in \operatorname*{supp}(\lambda)} \mu^{(0)}(\alpha,\lambda)\big( \mu^{(0)}(\alpha_1, \lambda) - \mu^{(0)}(\alpha_2, \lambda) \big)\dd \pi(\lambda) \nonumber\\ \geq& c_2\sqrt{\int_{\lambda \in \operatorname*{supp}(\lambda)} \big( \mu^{(0)}(\alpha_1, \lambda) - \mu^{(0)}(\alpha_2, \lambda) \big)^2\dd \pi(\lambda)} \end{align} \end{enumerate}

Assumption (ref) (ref) is a standard overlap condition and would be required for regular estimation of the ATT even if $\alpha_i$ were observed. Note that because we are interested in the average effect of treatment on the treated, we only require that the conditional probability of treatment is bounded below away from zero and not above away from one. Assumption (ref) (ref) restricts that the tail behavior of the error terms, requiring them to be sub-Gaussian. The assumption also imposes that the residuals are independent over time given the time-specific factors. This condition allows us to apply particular concentration inequalities in order to obtain fast rates of convergence.

Assumption (ref) (ref) imposes that the latent factors are independent and identically distributed and that the individual-specific factors have compact support. Compact support of the individual-specific latent factors helps to ensure that with high probability, for each treated individual $i$, there exist untreated individuals in the sample whose latent factors are close to those of $i$. Independence and identical distribution of the individual-specific latent factors follows if we understand individuals to be drawn identically and independently from the underlying population. Independence of $\lambda_t$ over time may be plausible if time periods are sufficiently far apart. This restriction on $\lambda_t$ allows us to apply concentration inequalities and ensure fast convergence of the pseudo-distance.

Assumption (ref) (ref) imposes that the functions $\mu^{(0)}$ and $p_t$ are bounded and vary smoothly with the individual-specific factors. This smoothness assumption ensures that individuals with similar latent factors also have similar conditional-mean potential outcomes and treatments. Assumption (ref) (ref) stipulates properties of the kernel $K$. Kernels that satisfy this condition include are common in the literature, with the Epanechnikov kernel a particularly prevalent choices.

Assumption (ref) (ref) relates the pseudo-distance to the distance between individual latent factors. The assumption imposes first, that the mean-squared distance between $\mu^{(0)}(\alpha_1,\cdot)$ and $\mu^{(0)}(\alpha_2,\cdot)$ is at least proportional to the Euclidean distance between $\alpha_1$ and $\alpha_2$. In addition, the assumption states that the population pseudo-distance is at least proportional to this mean-squared distance. Consider the special case in which $\mu^{(0)}(\alpha,\lambda)=\alpha'M\lambda$ for some fixed matrix $M$. In this case, the condition ((ref)) holds so long as the matrix $M\int_{\lambda\in\operatorname*{supp}(\lambda)}\lambda\lambda'\dd\pi(\lambda) M'$ is strictly positive definite. If this is the case, then the second condition ((ref)) holds if there is a $c>0$ so that for any $\alpha_1,\alpha_2\in\mathcal{A}$ we have $\frac{c}{\|\alpha_{1}-\alpha_{2}\|}(\alpha_{1}-\alpha_{2})\in \mathcal{A}$. This is true, for example, if $\mathcal{A}$ is a Euclidean ball centered at zero.

comment\begin{assu}{5}[Distance Measurement and Kernel] Assume the distance measurement and the kernel satisfy the following conditions: \begin{enumerate}[label=(\roman*), ref=(\roman*)] • For each $\delta > 0$, $\exists C_{\delta} > 0$ such that \begin{align*} d^2(i, j) = \sup_{\alpha_k \in \operatorname*{supp}(\alpha)} \left\vert \int_{\lambda \in \operatorname*{supp}(\lambda)} \mu^{(0)}(\alpha_k, \lambda) \left[ \mu^{(0)}(\alpha_i, \lambda) - \mu^{(0)}(\alpha_j, \lambda) \right] \dd \pi(\lambda) \right\vert > C_{\delta} \end{align*} almost surely for $\alpha_i$ and $\alpha_j$ satisfying $\Vert \alpha_i-\alpha_j \Vert_2 \geq \delta$. • $K: \mathbb{R}_{+} \rightarrow \mathbb{R}$ is weakly positive and strictly positive at zero, supported on compact set and bounded by $\bar{K} < \infty$. $K$ sastisfies $\gamma_K^1 := \int K(\vert x \vert) \dd x > 0$, $\gamma_K^2 := \int K^2(\vert x \vert) \dd x > 0$, and $\vert K(z) - K(z') \vert \leq \bar{K}' \vert z - z' \vert$ for all $z, z' \in \mathbb{R}_{+}$ for some $\bar{K}' > 0$. • $h \rightarrow 0$ as $n \rightarrow \infty$, $\dfrac{1}{Nh} \rightarrow 0$, $\dfrac{\log(N)}{T} \rightarrow 0$, and $\dfrac{1}{h} \sqrt{\dfrac{\log(N)}{T}} \rightarrow 0$. \end{enumerate} \end{assu}

thm{1}[Rate of Convergence for the Outcome Model] Suppose Assumptions (ref) to (ref) hold and $\frac{1}{h}\sqrt{\frac{log(N)}{T_{0}}} \to 0$ and $\frac{\log(N)}{h^{2d_\alpha} N} \to 0$, then \begin{align*} \max_{i \in [N]} \left\vert \hat{\mu}_{i,t}^{(0)} - \mu_{i,t}^{(0)} \right\vert = O_{p}\left(h+\sqrt{\frac{log(N)}{h^{d_\alpha}N}}\right)\\ \max_{i \in [N]} \left\vert \hat{p}_{i,t}^{(0)} - p_{i,t}^{(0)} \right\vert = O_{p}\left(h+\sqrt{\frac{log(N)}{h^{d_\alpha}N}}\right) \end{align*}

The convergence rate in Theorem (ref) is identical to that of the Nadaraya-Watson estimator with observed covariates up to a log term. Setting $h\propto N^{-1/(d_\alpha+2)}$ we get convergence rate $N^{-1/(d_\alpha+2)}\sqrt{log(N)}$ which, up to a log term, is the optimal rate in stone1980optimal under Lipschitz continuity. However, the rate at which $h$ can converge to zero is restricted by the condition that $\frac{1}{h}\sqrt{\frac{log(N)}{T_{0}}} \to 0$. Thus the theorem above allows us to achieve the rate in stone1980optimal, up to log terms, if and only if $T_{0}$ grows sufficiently quickly that $N^{2/(d_\alpha+2)} / T_{0} \to 0$. In order to obtain $\sqrt{N}$-consistency and $\sqrt{N}$ centered asymptotic normality of the doubly-robust ATT estimator, we require $\max\limits_{i\in[N]} \left\vert \hat{\mu}_{i,t}^{(0)} - \mu_{i,t}^{(0)} \right\vert = o_{p}(N^{-1/4})$ and similarly for $\hat{p}_{i,t}$. In the case of a one-dimensional individual latent factor ($d_\alpha=1$), Theorem (ref) ensures that this holds for an appropriate choice of $h$ if and only if $T_{0}$ grows quickly enough that $\frac{N^{1/2}log(N)}{T_{0}} \to 0$.

comment\begin{thm}{1}[Rate of Convergence for the Outcome Model] Assume Assumptions (ref)-(ref) hold. Then we have \begin{align*} \max_{i} \left\vert \hat{\mu}_{i,t}^{(w)} - \mu_{i,t}^{(w)} \right\vert = \mathcal{O}_p\left( h + \sqrt{\frac{1}{Nh}} + \frac{1}{h}\sqrt{\frac{\log(N)}{T_0}} \right) \end{align*} where $w \in \{0,1\}$ \end{thm} \begin{coro}{1}[Rate of Convergence for the Propensity Score Model] Assume Assumptions (ref)-(ref) hold. Then we have \begin{align*} \max_{i} \left\vert \hat{p}_{i,t} - p_{i,t} \right\vert = \mathcal{O}_p\left( h + \sqrt{\frac{1}{Nh}} + \frac{1}{h}\sqrt{\frac{\log(N)}{T_0}} \right) \end{align*} \end{coro} \begin{coro}{2}[Requirement for Pre-treatment History] To achieve a fast enough rate, i.e., faster than $N^{-1/4}$, it is required that $\dfrac{T_0}{N \log(N)} \rightarrow \infty$. \end{coro} \begin{remk}{2} Depending on the smoothness that is assumed, the rate of convergencce will also change correspondingly. For example, instead of the regualr Lipschitz assumption, we now have the Holder class: \begin{align*} \Sigma(\beta, L) = \left\{ \mu^{(0)} : \vert \nabla^s \mu^{(0)}(\alpha_i, \lambda) - \nabla^s \mu^{(0)}(\alpha_j, \lambda) \vert \leq L \Vert \alpha_i - \alpha_j \Vert , for all $s$ such that $|s| = \beta - 1$ \right\} \end{align*} and the bi-Lipschitz part still need to be satisfied: \begin{align*} L_0' \Vert \alpha_i - \alpha_j \Vert \leq &\vert \mu^{(0)}(\alpha_i, \lambda) - \mu^{(0)}(\alpha_j, \lambda) \vert \end{align*} In this specific assumption, the convergence rate will become \begin{align*} \max_{i} \left\vert \hat{\mu}_{i,t}^{(w)} - \mu_{i,t}^{(w)} \right\vert = \mathcal{O}_p\left( h^\beta + \sqrt{\frac{1}{Nh}} + \frac{1}{h}\sqrt{\frac{\log(N)}{T_0}} \right) \end{align*} and the requirement for the number of the pre-treatment time periods $T_0$ will reduce. Furthermore, if we allow multi-dimensional latent factor $\alpha$, i.e., $\operatorname*{supp}(\alpha) \subseteq \mathcal{A}$ where $\mathcal{A}$ is a compact subset of $\mathbb{R}^{d_{\alpha}}$ with $d_{\alpha} > 1$, then the convergence rate will become \begin{align*} \max_{i} \left\vert \hat{\mu}_{i,t}^{(w)} - \mu_{i,t}^{(w)} \right\vert = \mathcal{O}_p\left( h^\beta + \sqrt{\frac{1}{Nh^{d_{\alpha}}}} + \frac{1}{h}\sqrt{\frac{\log(N)}{T_0}} \right) \end{align*} by adjusting the kernel based estimators with the assumption of $\dfrac{1}{Nh^{d_{\alpha}}} \rightarrow 0$. \end{remk}

Inference Framework with Rate Double Robustness

Theorem (ref) below applies well-established ideas from the literature on double machine learning to this setting. The result applies to the estimator with cross-fitting as specified in Algorithm (ref).

thm{2}[Asymptotic Normality] Suppose Assumptions (ref)-(ref) all hold. In addition, suppose $\max_{i\in[N]} \left\vert \hat{p}_{i,t}-p_{i,t} \right\vert = o_{p}(1)$, $\max_{i\in[N]} \left\vert \hat{\mu}_{i,t}-\mu_{i,t} \right\vert = o_{p}(1)$, $\max_{i\in[N]} \left\vert (p_{i,t}-\hat{p}_{i,t})(\mu_{i,t}^{(0)}-\hat{\mu}_{i,t}^{(0)}) \right\vert = o_{p}\left(N^{-1/2}\right)$. Finally, suppose that for each $k$, the size of fold $\mathcal{I}_{k}$ grows at the same rate as the same size, that is, $\frac{|\mathcal{I}_{k}|}{N}\to c_{k}$ for some finite $c_{k}>0$. It follows that \begin{align*} \sqrt{N}\hat{\mathrm{ATT}}_{t}=\frac{\sqrt{N}}{N_{1,t}}\sum_{i=1}^{N}\left(Y_{i,t}w_{i,t}-\frac{(1-w_{i,t})Y_{i,t}p_{i,t}+(w_{i,t}-p_{i,t})\mu_{i,t}^{(0)}}{1-p_{i,t}}\right) + o_{p}(1) \end{align*} and thus under the stated conditions, $\sqrt{N} \left(\hat{\mathrm{ATT}}_{t}-\mathrm{ATT}_{t}\right) = N\left(0,V_t\right) + o_{p}(1)$.

The first result in Theorem (ref) relates the estimate $\hat{\mathrm{ATT}}_t$ to an infeasible estimate in which the first stage estimates of $\mu_{i,t}^{(0)}$ and $p_{i,t}$ are replaced by their true values. In particular, the Theorem states that under certain conditions, the difference between $\hat{\mathrm{ATT}}_t$ and the infeasible estimate disappears strictly faster than $N^{-1/2}$. It then follows from standard results that $\hat{\mathrm{ATT}}_t$ is root-$N$ asymptotically normal and centered at $ATT_t$.

The proof of Theorem (ref) proceeds by similar steps to the results in chernozhukov2018double for DML2 estimators. A complicating factor is that, unlike in standard DML2 estimation of the ATT, in our setting the estimates $\hat{\mu}^{(0)}_{i,t}$ and $\hat{p}_{i,t}$ are constructed using the pseudo-distance which includes both outcome and treatment data for individual $i$ in the pre-treatment period. That is, unlike in standard DML2 ATT estimation, the first-stage estimates $\hat{\mu}^{(0)}_{i,t}$ and $\hat{p}_{i,t}$ depend on some of the outcome data from the fold that contains individual $i$. However, using the conditional serial independence of the errors in Assumption (ref), it is straight-forward to accommodate this dependence.

Theorem (ref) requires not only that the first stage estimates are uniformly consistent over the sample, but that the product of the estimation error in $\hat{\mu}^{(0)}_{i,t}$ and $\hat{p}_{i,t}$ goes to zero uniformly strictly faster that root-$N$. A sufficient condition, mentioned earlier in this section, is that each of these estimates converges strictly faster than $N^{-1/4}$. The theorem further requires that each fold grows at the same rate as the sample data. This is satisfied if the number of folds is fixed and all of the folds are of equal size. In recent work, velez2024asymptoticpropertiesdebiasedmachine provides conditions under which DML2 estimates are centered root-$N$ asymptotically normal even if leave-one-out cross-fitting is applied, and so it may be possible to weaken this condition on the fold size in Theorem (ref).

commentTo achieve rate double robustness, different sets of regularity conditions or methodological modifications are required. To avoid over-fitting bias while still attaining the $\sqrt{N}$ rate of the doubly robust estimator, two approaches can be considered. The first approach involves assuming strong Donsker conditions, i.e., strong smoothness conditions, for both the nuisance functions of the outcome model and the propensity score (e.g., robins2008higher). The second approach utilizes sample splitting techniques to retrieve independence and eliminate over-fitting bias (e.g., chernozhukov2018double, abadie2024doubly). Given that one advantage of the proposed method is to not rely on structural assumptions about the underlying data generating process and to impose minimal smoothness conditions, the first approach of imposing Donsker conditions defeats this purpose. Therefore, we adopt the sample splitting approach. Additionally, employing the sample splitting method with a cross-fitting algorithm as suggested in the aforementioned literature and in Algorithm (ref) can help us recover any efficiency loss in the estimation procedure within each fold. However, achieving rate double robustness is particularly challenging in the presence of latent confounders. abadie2024doubly provides an inference procedure for doubly robust estimators specifically designed for such situations. We will illustrate the assumptions and theorems from abadie2024doubly adapted to our notation. Also, it should be noted that these assumptions are consistent with the assumptions of our proposed method. \begin{assu}{AADS.1}[Positivity; Assumption 1 in abadie2024doubly] \protected@write \@auxout {\string \newlabel {asm:AADS.1}{{AADS.1}{\thepage}{AADS.1}{asm:AADS.1}}} \hypertarget{asm:AADS.1} The unknown propensity score matrix $P$ is such that \begin{align*} \delta \leq p_{i,t} \leq 1-\delta, \end{align*} for all $i \in \{1, 2, ..., N\}$ and $t \in \{1, 2, ..., T\}$, where $0 < \delta \leq 1/2$ is a constant. \end{assu} \begin{assu}{AADS.2}[Zero-mean, independent, and sub-Gaussian noise; Assumption 2 in abadie2024doubly] \protected@write \@auxout {\string \newlabel {asm:AADS.2}{{AADS.2}{\thepage}{AADS.2}{asm:AADS.2}}} \hypertarget{asm:AADS.2} Fix any $t \in \{1, 2, ..., T\}$. Then, \begin{enumerate}[label=(\roman*), ref=(\roman*)] • $\{ (u_{i,t}, v_{i,t}, \epsilon_{i,t}): i \in \{1, 2, ... ,N\} \}$ are mean zero and independent (across $i$), • $(u_{i,t}, v_{i,t}) \perp \!\!\! \perp \epsilon_{i,t}$ for every $i \in \{1, 2, ... ,N\}$, and • $u_{i,t}$ and $v_{i,t}$ have sub-Gaussian norm bounded by a constant $\bar{\sigma}$ for every $i \in \{1, 2, ... ,N\}$. \end{enumerate} \end{assu} \begin{assu}{AADS.3}[Assumption 3 in abadie2024doubly] \protected@write \@auxout {\string \newlabel {asm:AADS.3}{{AADS.3}{\thepage}{AADS.3}{asm:AADS.3}}} \hypertarget{asm:AADS.3} The estimated propensity score matrix $\hat{P}$ is such that \begin{align*} \bar{\lambda} \leq \hat{p}_{i,t} \leq 1-\bar{\delta}, \end{align*} for all $i \in \{1, 2, ..., N\}$ and $t \in \{1, 2, ..., T\}$, where $0 < \bar{\delta} \leq \delta$ is a constant. \end{assu} \begin{assu}{AADS.4}[Assumption 4 in abadie2024doubly] \protected@write \@auxout {\string \newlabel {asm:AADS.4}{{AADS.4}{\thepage}{AADS.4}{asm:AADS.4}}} \hypertarget{asm:AADS.4} Fix any $t \in \{1, 2, ..., T\}$. There exists a non-empty partition $(\mathcal{I}_1, \mathcal{I}_2)$ of the units $\{1, 2, ..., N\}$ such that \begin{align*} \left\{\left(\hat{p}_{i,t}, \hat{\mu}_{i,t}^{(w)}\right)\right\}_{i \in \mathcal{I}_j} \perp \!\!\! \perp \left\{\epsilon_{i,t}\right\}_{i \in \mathcal{I}_j} \end{align*} and \begin{align*} \left\{\hat{p}_{i,t}\right\}_{i \in \mathcal{I}_j} \perp \!\!\! \perp \left\{\left(\epsilon_{i,t}, \hat{\mu}_{i,t}^{(w)}\right)\right\}_{i \in \mathcal{I}_j} \end{align*} for every $w \in \{0,1\}$ and $j \in \{1, 2\}$ \end{assu} Assumptions (ref) and (ref) define the overlapping conditions between the population and its sample estimates for propensity scores. Assumption (ref) ensures the independence of error terms across individuals in both the outcome and propensity score models. Additionally, it restricts the tail behavior of the error term distribution, which is crucial for determining the bounds on the rate. Assumption (ref) specifies the allowable sample splitting algorithm that retrieves independence for eliminating over-fitting bias. This assumption is satisfied by the proposed $K$-fold sample splitting with a cross-fitting algorithm. With these assumptions, the asymptotic normality for the doubly robust estimator is guaranteed by the following theorem. \begin{thm}{AADS.1}[Asymptotic Normality for Doubly Robust Estimator; Theorem 2 in abadie2024doubly] \protected@write \@auxout {\string \newlabel {thm:AADS.1}{{AADS.1}{\thepage}{AADS.1}{thm:AADS.1}}} \hypertarget{thm:AADS.1} Denote \begin{align*} \mathcal{E}\left(\hat{P}\right) := \frac{\max\limits_{t \in \{1, ..., T\}} \left\{ \sum\limits_{i \in \{1, ..., N\}} \left( \hat{p}_{i,t} - p_{i,t} \right)^2 \right\}^{1/2}}{\sqrt{N}} \end{align*} and \begin{align*} \mathcal{E}\left(\hat{\mu}\right) := \sum_{w \in \{0, 1\}} \mathcal{E}\left(\hat{\mu}^{(w)}\right) := \sum_{w \in \{0, 1\}} \frac{\max\limits_{t \in \{1, ..., T\}} \left\{ \sum\limits_{i \in \{1, ..., N\}} \left( \hat{\mu}_{i,t}^{(w)} - \mu_{i,t}^{(w)} \right)^2 \right\}^{1/2}}{\sqrt{N}} \end{align*} Suppose Assumptions (ref)-(ref) and the following conditions hold, \begin{enumerate}[label=(\roman*), ref=(\roman*)] • $\mathcal{E}(\hat{P}) = o_p(1)$ and $\mathcal{E}(\hat{\mu}) = o_p(1)$ • $\mathcal{E}(\hat{P}) \mathcal{E}(\hat{\mu}) = o_p(N^{-1/2})$ • For every $i \in \{1, ..., N\}$ and $t \in \{1, ..., T\}$, let $\sigma_{i,t}^{(0)}$ and $\sigma_{i,t}^{(1)}$ be the standard deviation of $u_{i,t}$ and $v_{i,t}$, respectively. The sequence \begin{align*} \bar{\sigma}_t^2 := \frac{1}{N} \sum_{i \in \{1, ..., N\}} \frac{\left(\sigma_{i,t}^{(1)}\right)^2}{p_{i,t}} + \frac{1}{N} \sum_{i \in \{1, ..., N\}} \frac{\left(\sigma_{i,t}^{(0)}\right)^2}{1 - p_{i,t}}, \end{align*} is bounded away from zero as $N$ increases. \end{enumerate} Then, for all $t \in \{1, ..., T\}$, \begin{align*} \sqrt{N} \left(\hat{\mathrm{ATE}}_{\cdot,t}^{\mathrm{DR}} - \mathrm{ATE}_{\cdot,t}\right) / \bar{\sigma}_t \xrightarrow{d} \mathcal{N}(0,1), \end{align*} as $N \rightarrow \infty$. \end{thm} While the theorem is established with the average treatment effect being the target estimand, as mentioned in abadie2024doubly, other causal estimands, such as our main parameter of interest $\mathrm{ATT}_t$, can also achieve asymptotic normality for inference. This is attainable when having proper doubly robust score functions and the appropriate form of variance estimator, as the result is immediate. Employing Theorem (ref) along with the rate of convergence result in Theorem (ref), we establish the asymptotic inference for the proposed estimator: \begin{thm}{2}[Asymptotic Inference for the Proposed Estimator] Assume Assumptions (ref)-(ref) hold, $\frac{N^{1/2}log(N)}{T_{0}}\to0$ and $h$ is chosen so that $\frac{1}{h} \sqrt{\frac{\log(N)}{T_0}} \to 0$, it follows that the product rate of the outcome model and the propensity score model becomes faster than $N^{-1/4}$. With the results from Theorem (ref), the aforementioned doubly robust estimator for $\mathrm{ATT}_t$, outlined in Algorithm (ref), achieves asymptotic normality. \end{thm}

Simulation Study

In this section, we illustrate the finite sample properties of the proposed method in a number of numerical experiments. We also compare it with (1) the workhorse TWFE approach, and (2) the method proposed by feng2024causal.

Comparison with the TWFE approach

The TWFE approach is widely used in applied research and is designed to adjust for unobserved additive individual and time fixed effects. In this section, we demonstrate that the proposed method performs comparably well in settings in which the TWFE approach is valid. We also show that the proposed method maintains good finite sample properties in a more complicated setting in which the TWFE approach fails.

Consider a large panel data setting with individuals labeled $i = 1, ..., N$ over time periods $t = 1, ..., T$, where $N \in \{50, 250\}$ and $T_0 \in \{50, 250\}$ with only one post-treatment time period in this numerical study. For the outcome model, we consider the following data generating processes:

align*[align* omitted — 259 chars of source]

where the individual latent factor is distributed as $\alpha_i \sim \rm{Uniform}(-1, 1)$, the time latent factor as $\lambda_t \sim \rm{Uniform}(-1, 1)$, and $u_{i,t} \sim N(0, 0.5^2)$. Here, the target estimand, the average treatment effect on the treated, is $\mathrm{ATT}_t = 0.5$. Models 1 and 2 are standard additive and interactive fixed effects models commonly used in panel data.

The treatment assignment mechanism depends on unobserved latent characteristics for each $i$. Specifically, the propensity score$p_{i,t}$ follows the model

align*[align* omitted — 65 chars of source]

With the configuration of factors located on compact supports, it results in roughly $p_{i,T} \in [0.2689, 0.7311]$, which satisfies the overlapping condition required for the proposed method.

We simulate $500$ replications from these models. In these simulations, we utilize the Epanechnikov kernel, that satisfies the requirements of Theorem 1 and is defined as $K(x) = \frac{3}{4} (1 - x^2) \cdot \mathbf{1}\{\vert x \vert \leq 1 \}$. We use this kernel for both the outcome and propensity score imputations in our proposed method. Using this kernel, we compared the performance of different methods: (i) TWFE; (ii) the local PCA method proposed by feng2024causal with leave-one-out cross validation for the choice of nearest neighbors\footnote{The implementation of the method proposed by feng2024causal follows the replication files available on the author's website: \url{https://github.com/yingjieum/replication-Feng_2024}.}; (iii) an infeasible doubly robust estimates in which the pseudo-distance is replaced by the oracle $\ell_2$ distance between the true $\alpha$s and in which the estimated propensity score is replaced with the true propensity score; (iv) an infeasible doubly robust estimates in which the pseudo-distance is replaced by the $\ell_2$ distance between the true $\alpha$s but in which the propensity score is unknown and estimated using the oracle distance; (iv) our proposed method.

In the numerical experiments, the bandwidth is selected between $0.05$ and $5$. velez2024asymptoticpropertiesdebiasedmachine suggests that leave-one-out estimation is the optimal cross-fitting procedure for DML2 estimators in terms of both bias and the second-order asymptotic mean squared error under certain conditions. Hence, in addition to the results from 2-fold cross-fitting and cross-validation bandwidth selection, we also implemented the proposed method using the leave-one-out procedure for this numerical exercise. In Table (ref) and (ref), we report the 95% confidence interval coverage and other statistics across $500$ replications in our simulation for each outcome model and each method. Additionally, we note that due to the non-existence of moments, some mean versions of the statistics can be distorted by outliers in the replications. This distortion is more pronounced when the number of individual $N$ is small. For this reason, we report the median absolute deviation and median confidence interval length rather than the means.

table[table omitted — 6,106 chars of source]

For model 1, the additive fixed effects model, the proposed method shows performance similar to the two-way fixed effects method in terms of 95% confidence interval coverage, which is under an ideal data generating process for TWFE. Additionally, the proposed method, which incorporates the pseudo distance, effectively captures the distance between unobserved latent characteristics, which is evinced by the performance of our method relative to the oracle estimators which use the true propensity score and/or the true distance in the $\alpha$s. However, there is some loss of efficiency when the sample size $N$ is small. This efficiency loss can be mitigated by increasing the number of folds or by using the leave-one-out procedure for bandwidth selection.

table[table omitted — 6,109 chars of source]

For model 2 with the interactive fixed effects, the TWFE method struggles to achieve approximately correct 95% confidence interval coverage. This problem worsens as $N$ increases and there are more individual latent characteristics to capture. In contrast, the proposed method demonstrates relatively good coverage for all different choices of $N$ and $T$.

Comparison with Method Proposed by feng2024causal

We compare the performance of our approach with the method proposed by feng2024causal, which, to the best of our knowledge, is one of the few available inferential methods for average causal effects in nonlinear factor models. In this subsection, we use precisely the same data generating processes employed in feng2024causal but we consider different values of $N$ and $T$. The implementation is based on the replication code provided by the author. These numerical exercises are designed to investigate how the considered methods perform under relatively small sample sizes which may be more typical of microeconomic applications.

Again, consider a large panel data setting with individuals labeled $i = 1, ..., N$ over time periods $t = 1, ..., T$, where $N \in \{50, 250\}$ and $T_0 \in \{50, 250, 1000\}$ with only one post-treatment time period. Our target estimand is the counterfactual mean $\mathbb{E}[Y_{i,T}(0) | w_{i,T} = 1]$ which is the object of interest in the simulations of feng2024causal. Note that this object is readily calculated from our $\hat{\mathrm{ATT}}_t$ estimate because $\mathbb{E}[Y_{i,T} | w_{i,T} = 1]$ is trivially identified. We consider the following DGPs for the untreated potential outcomes:

align*[align* omitted — 397 chars of source]

for $t \in \{1, ..., T_0\}$. In the post-treatment time period, all three models follow

align*[align* omitted — 128 chars of source]

where the individual latent factor $\alpha_i \sim \rm{Uniform}(-0, 1)$ and time latent factor $\lambda_t \sim \rm{Uniform}(0, 1)$, $u_{i,t} \sim N(0, 0.5^2)$, and $\epsilon_{i,t} \sim N(0, 1)$. The propensity score is given by

align*[align* omitted — 152 chars of source]

For this DGP, we roughly have $p_{i,T} \in [0.4378, 0.6792]$, which satisfies the required overlapping condition.

Model 3 is a nonlinear factor model. Models 4 and 5 are the Gaussian kernel and exponential kernel. A similar simulation design is also used fernandez2021low.

The results are provided in Tables (ref), (ref), and (ref) below. We consider the same methods and report the same statistics as for the previous numerical experiments. The number of replications is $500$.

table[table omitted — 7,403 chars of source]

table[table omitted — 7,395 chars of source]

table[table omitted — 7,398 chars of source]

For models 3 to 5, which involve different nonlinear data generating processes, the performance of both the proposed method and the local PCA method proposed by feng2024causal differs across various combinations of $N$ and $T_0$. In models 3 and 4 with the nonlinear and Gaussian kernel DGPs, the proposed method provides good coverage across for different sample sizes. Specifically, we find that the cross-fitting procedure substantially improves the confidence interval coverage, and the leave-one-out procedure improves the accuracy of point estimation. Here, the proposed method either performs similarly or better than the local PCA method. In contrast, the local PCA method tends to underperform when the pre-treatment history is relatively short compared to the number of individuals. This issue gets amplified in model 5 with the exponential kernels DGP. In configurations with $N = 250$, the local PCA method requires a larger pre-treatment history of $T_0$ being $1000$ to achieve satisfactory confidence interval coverage. On the other hand, the proposed method with the leave-one-out procedure achieves effective coverage with a shorter $T_0$ requirement.

In summary, the proposed method provides a simple and transparent estimation and inference procedure for different types of causal estimands with their associated doubly robust score function. With the aforementioned algorithm, we can effectively collect individuals who are similar to each other in terms of their latent characteristics in a more transparent manner. This simulation study showcases that the proposed method can match or even outperform the workhorse two-way fixed effects model and the method proposed by feng2024causal across different data generating processes, requiring fewer amount of pre-treatment history in certain settings.

\addcontentsline{toc}{section}{References}