2.1 The Average Treatment Effect Localization

We consider a panel data setting with $N+1$ units observed over T periods. The treated unit (the first unit) is untreated for $T_0$ period and remains treated thereafter. The N control units remain untreated for all periods. Let $Y^I_{it}$ and $Y^{N}_{it}$ be the potential outcomes of unit i in period t with and without treatment, respectively. The individual treatment effect for the ith unit at time t is defined as

$$ \begin{align*} \alpha_{it} = Y^{I}_{it} - Y^{N}_{it}. \end{align*} $$

The fundamental challenge for causal inference is that we cannot observe both outcomes $Y^{I}_{it}$ and $Y^{N}_{it}$ simultaneously. Thus, the observed outcomes $Y_{it}$ are related to potential outcomes as $Y_{it} = D_{it}Y^{I}_{it} + (1 - D_{it})Y^{N}_{it}$ , where $D_{it} =1$ if $i=1$ and $t>T_0$ and $D_{it}=0$ otherwise.

In the absence of treatment, outcomes are generated by the following nonparametric model:

(3) $$ \begin{align} Y_{it}^N = h_t(\eta_{it}, X_{it}) + u_{it}, \text{ for } i=1,\ldots,N+1,\text{ } t=1,\ldots,T, \end{align} $$

where the covariates $(\eta _{it}, X_{it})$ are only partially observed. Specifically, $\eta _{it}$ is a vector of unobserved variables, and $X_{it}$ is a vector of observed covariates. This model captures the nonlinearity and time-varying effects from partially observed covariates $(\eta _{it}, X_{it})$ through the time-varying unknown function $h_t(\cdot ,\cdot )$ . The term $u_{it}$ represents the idiosyncratic error, which is independent of $(\eta _{it},X_{it})$ .

After the treated unit receives a policy intervention or treatment, the treated unit outcome for the post-treatment periods is

(4) $$ \begin{align} Y_{1t} = Y^I_{1t} = h_t(\eta_{1t}, X_{1t}) + \alpha_{1t} + u_{1t},\text{ for } t = T_0+1,\ldots,T. \end{align} $$

Given that the treated unit receives the treatment after $T_0$ , the conventional causal analysis often focuses on the long-term average treatment effect on the treated. It is formally defined as

$$ \begin{align*} \alpha_{\text{long-term}} = \frac{1}{T_1}\sum_{t=T_0+1}^T \mathbb{E}(Y^I_{1t} - Y^N_{1t}), \end{align*} $$

where $T_1 = T - T_0$ is the number of post-treatment periods.

In contrast to the long-term average treatment effect, this article focuses on the policy intervention that is localized at the time of treatment $T_0$ , to facilitate the analysis of short-term policy impacts. The parameter of interest is the ATEL, defined as

(5) $$ \begin{align} \alpha = \frac{1}{T_1}\sum_{t=T_0+1}^T \mathbb{E}(Y^I_{1t} - Y^N_{1t})K_h \left(\frac{t-T_0}{T_1}\right), \end{align} $$

where $K_h(\cdot )$ is a kernel function with bandwidth parameter h. The kernel function assigns weights to post-treatment periods, localizing the treatment effect around $T_0$ , so $\alpha $ represents the short-term effect soon after $T_0$ . In particular, if $K_h(\cdot )$ follows a uniform kernel, $\alpha $ simplifies to the average treatment effect over the $\lfloor T_1 h \rfloor $ post-treatment periods. Additional details about the kernel function are provided in (19).

To motivate why short-term policy evaluation is important, our empirical example examines how RTC laws affect violent crime shortly after they are implemented. These laws usually allow people to carry concealed firearms in public with a permit. Supporters believe that RTC laws can help citizens protect themselves and deter crime, whereas critics argue that more guns in public increase the risk of violence and accidental shootings. Although many empirical studies have investigated the long-term average treatment effects of RTC laws, the findings remain inconclusive (e.g., Lott Jr. and Mustard, Reference Lott and Mustard1997; Aneja et al., Reference Aneja, Donohue and Zhang2011; Donohue et al., Reference Donohue, Aneja and Weber2019). Given that the United States has one of the highest rates of firearm violence among developed countries, it is crucial to evaluate RTC laws promptly.

The short-term policy impact especially matters to policymakers, such as governors, legislators, and mayors, who often serve short terms of 2–4 years. If an RTC law is linked to a sudden rise in violent crime soon after it goes into effect, it may cause a political backlash, no matter what the long-term data eventually show. Public perception also plays a large role in decisions, sometimes more so than academic findings. A short-term spike in gun violence could lead to the repeal or tightening of RTC laws. Therefore, policymakers need to pay close attention to the policy’s short-term impact when deciding whether it is effective.

In practice, $\alpha $ represents the short-term policy effect right after $T_0$ , which is useful for policymakers and business leaders who must make prompt decisions. In contrast, researchers often focus on the long-term average treatment effect on the treated. While $\alpha _{\text {long-term}}$ considers the treatment effect across all post-treatment periods, it can be diminished over time as treated units adapt or respond strategically, making it harder to measure the policy’s direct impact. For example, in February 2025, the Trump administration reduced the USAID budget, forcing multiple labs to shut down and causing an immediate drop in academic output. Over time, academic institutions responded with various mitigation strategies, such as fundraising and layoffs, making it difficult to capture the policy’s original impact.

As reported by The New York Times, these cuts at Johns Hopkins University, funded largely by USAID grants for public health and agricultural research, led to the elimination of research projects on tuberculosis, AIDS, and cervical cancer. In response, the university laid off more than 2,000 employees to cope with the reduced budget (Saul, Reference Saul2025). For researchers examining how USAID budget reductions influence academic output, these institutional measures and other external adjustments can distort the long-term average treatment effect, obscuring the policy’s direct influence in a long run. By contrast, our main objective, $\alpha $ , measures the policy intervention shortly after $T_0$ , which remains largely unaffected by these subsequent actions and countermeasures.

2.2 Low-Rank Time-Varying Factor Structure: Approximation and Estimation

In this section, we discuss how the nonparametric model (3) can be approximated by a low-rank time-varying factor model with a slowly growing rank, and we outline its estimation methods.

We motivate the sieve approximation for the nonlinear time-varying function $h_t(\eta _{it},X_{it})$ with two considerations. First, standard nonparametric regression is infeasible due to the presence of unobserved covariates $\eta _{it}$ . Unobserved heterogeneity is typically handled via linear transformations (e.g., first-differencing) in linear panel data models to remove unit- or time-invariant components. However, because $\eta _{it}$ varies across both units and time and $h_t(\cdot )$ imposes a nonlinear structure, there is no general rule of transformation to remove $\eta _{it}$ . Second, the sieve expansion with truncation allows us to approximately represent the matrix formed by the elements $h_t(\eta _{it},X_{it})$ using a low-rank matrix. This low-rank structure is fundamental to our procedure, as it enables us to recover the relevant low-dimensional components while avoiding fully nonparametric estimation of $h_t(\cdot )$ .

Formally, let $\Theta $ be the $(N+1) \times T$ matrix where the $(i,t)$ element of $\Theta $ is $h_t(\eta _{it},X_{it})$ . Suppose the unknown function $h_t(\eta _{it},X_{it})$ has the following sieve approximation:

(6) $$ \begin{align} h_t(\eta_{it},X_{it}) = \sum_{j=1}^J \phi_j(\eta_{it},X_{it}) f_{jt} + r_{it} = \beta^{\prime}_{it} F_t + r_{it}, \end{align} $$

where $\beta _{it} = (\phi _1(\eta _{it},X_{it}),\ldots ,\phi _{J}(\eta _{it},X_{it}) )' \in \mathbb {R}^J$ is a set of sieve transformation of $\eta _{it},X_{it}$ using $\phi _j$ as the basis function (e.g., B-spline, Fourier series, wavelets, and polynomial series), $F_t = (f_{1t},\ldots ,f_{Jt})'$ is the vector of sieve coefficients, and $r_{it}$ is the sieve approximation error. J denotes the number of sieve terms which grow slowly as $N,T \rightarrow \infty $ . Let $\Theta _0$ be the $(N+1) \times T$ matrix with elements $\beta _{it}' F_t$ , and R be the $(N+1) \times T$ matrix of $r_{it}$ . In matrix notation, we represent (6) as

$$ \begin{align*} \Theta= \Theta_0 + R. \end{align*} $$

For sufficiently smooth functions $\{h_t(\cdot )\}$ , the sieve approximation error $r_{it} \xrightarrow {} 0$ as $J \rightarrow \infty $ , so the matrix $\Theta $ can be well approximated by the low-rank matrix $\Theta _0$ . For example, suppose $h_t$ belongs to a Hölder class with parameter $b,C,d> 0$ :

$$ \begin{align*} \{ h: \| D^b h(x) - D^b h(y)\| \leq C \|x - y\|^{d}, \text{for all } x,y\}. \end{align*} $$

Then, for common bases such as B-splines, trigonometric polynomials, or polynomial series, the approximation error satisfies

$$ \begin{align*} \max_{it} |r_{it}| < C J^{-a},\quad a = (b+d) / (dim(\eta_{it}) + dim(X_{it})), \end{align*} $$

which can be made arbitrarily small for sufficiently smooth functions even if J grow slowly. The smoothness condition justifies the low-rank approximation of $h_t(\eta _{it},X_{it})$ because it guarantees that the approximation error vanishes as J increases, allowing the truncation of the sieve expansion to the first J terms. In many empirical applications, the function $h_t(\cdot )$ represents the average outcome conditional on the characteristics $\eta _{it}$ and $X_{it}$ . It is plausible to assume that such function is sufficiently smooth so that it belongs to the Hölder class. Then it is well-known that such functions can be well-approximated by sieve functions with a small number of sieve bases, yielding a low-rank structure.

Combining (3) and (6), the untreated potential outcomes follow an approximate time-varying factor structure:

(7) $$ \begin{align} Y^N_{it} = \beta^{\prime}_{it} F_t + r_{it} + u_{it}. \end{align} $$

Several studies have considered low-rank factor structure to represent nonparametric function using sieve approximation (see, e.g., Fernández-Val et al., Reference Fernández-Val, Freeman and Weidner2021; Chernozhukov et al., Reference Chernozhukov, Hansen, Liao and Zhu2023; Choi et al., Reference Choi, Kwon and Liao2024). However, instead of a time-varying factor structure, these analyses approximate functions as pure factor structure, which restrict the unobserved covariates either be unit-invariant or time-invariant, while reasonable in certain applications, may be unrealistic in many practical settings.

For simplicity, we assume that the rank of the low-rank factor structure, J, is known. Since the time-varying function is approximated through a low-rank factor structure via a sieve representation, the rank equals the sieve dimension. When the rank is unknown in an empirical study, consistent estimators of rank J can be obtained using a conditioning argument (see, e.g., Bai and Ng, Reference Bai and Ng2002; Ahn and Horenstein, Reference Ahn and Horenstein2013; Su and Wang, Reference Su and Wang2017).

Building on the factor representation, we separately discuss the estimation of factors and factor loadings. To estimate factors, we adopt the DP approach proposed by Fan and Liao (Reference Fan and Liao2022), which reduces dimensionality by projecting data onto low-dimensional subspaces using a diversified weight matrix W. Unlike conventional PCA, which requires a strong signal-to-noise ratio to identify factor estimators, the DP method estimates factors through weighted averages of observations rather than relying on eigenvectors. This projection is conducted across individuals, making the estimator robust to serial conditions. The factor estimator is given by

(8) $$ \begin{align} \hat{F}_t = \frac{1}{N} \sum_{i=2}^{N+1} W_{it} Y^N_{it}, \text{ for } t = 1,\ldots,T, \end{align} $$

where $W_{it}=(w_{it1},\ldots ,w_{itJ})' \in \mathbb {R}^J$ that should satisfy:

1. (i) $W_{it}$ is independent of $u_{it}$ .

2. (ii) $W_{it}$ is correlated with $\beta _{it}$ .

Conditions (i) and (ii) are diversified weight’s analogous to the exclusion restriction and relevance restriction in instrument variables. They are necessary to consistently estimate $F_t$ up to affine transformation. Substituting (3) and (6) into $\hat {F}_t$ gives

(9) $$ \begin{align} \begin{aligned} \hat{F}_t &= \left(\frac{1}{N} \sum_{i=2}^{N+1} W_{it} \beta^{\prime}_{it}\right) F_t + \frac{1}{N}\sum_{i=2}^{N+1} W_{it} r_{it} + \frac{1}{N} \sum_{i=2}^{N+1} W_{it} u_{it} \\ & = H_t F_t + e_t, \end{aligned} \end{align} $$

where $H_t$ is the $J \times J$ time-varying transformation matrix, and $e_t$ is the estimation error that includes the sieve approximation error $r_{it}$ and the idiosyncratic error $u_{it}$ . Under the mild conditions of W, $e_t \xrightarrow {p}0$ as $N \rightarrow \infty $ .

A natural question is how to acquire the diversified weight matrix that satisfies these two conditions. A compelling choice is through observed characteristics $X_{it}$ , which have been extensively explored in factor model. For instance, Connor et al. (Reference Connor, Hagmann and Linton2012) introduced a characteristic-based factor model to analyze stock returns. This model incorporates factor loadings influenced by firm-specific attributes, such as market capitalization and book-to-price ratio. Fan et al. (Reference Fan, Liao and Wang2016) propose a more general form in the factor loading that can be explained by both observed covariates and unobserved components.

We can then define the diversified weight $W_{it}$ as the transformation of the observed characteristics $X_{it}$ with a set of functions $\varphi _j(\cdot )$ :

$$ \begin{align*} w_{itj} = \varphi_j(X_{it}), \text{ for } i\leq N+1,\text{ }t\leq T, \text{ } j\leq J, \end{align*} $$

where $\varphi _j(\cdot )$ is a sieve basis function, such as B-spline, trigonometric series, and polynomial series. Given that the factor loading is expressed as the sieve transformation of covariates (i.e., $\beta _{it} = \phi (\eta _{it},X_{it})$ ) and $u_{it}$ is independent of $X_{it}$ , both conditions required for $W_{it}$ are satisfied. However, certain choices of $W_{it}$ may violate these conditions. For instance, Fan and Liao (Reference Fan and Liao2022) suggest constructing diversified weights from the initial outcome $Y_{i1}$ , which has the advantage of not requiring additional observed characteristics. Since $Y_{i1}$ contains the information of $u_{i1}$ , if $u_{it}$ is serially correlated, the resulting weights can be correlated with $u_{it}$ , and the exogeneity condition fails.

The DP approach used here differs from that in Fan and Liao (Reference Fan and Liao2022) in several ways. First, their DP method is robust to over-estimating the number of factors, it does not require exact knowledge of this number for consistent factor estimation, as long as the rank of the diversified weight, $W_{it}$ , is at least as large as the true number of factors. Intuitively, since DP constructs factors as cross-sectional weighted averages rather than eigenvectors, over-estimating the number of factors typically adds only redundant projections onto the same factor space. This differs from PCA, where additional eigenvectors correspond to non-spiked eigenvalues and are often driven by idiosyncratic noise. In contrast, our framework assumes that the number of factors, J, is known, allowing us to construct a diversified weight matrix with rank equal to J. Consequently, our DP method can be viewed as an exactly identified IV approach. It is interesting to extend this framework to overidentified IV settings. We examine robustness to over-estimating the number of factors in the simulation and empirical sections and leave the associated theoretical development for future work.

Second, Fan and Liao’s (Reference Fan and Liao2022) diversified weights are restricted to being time-invariant, with the recommended choice for $W_i$ based on individual-specific characteristics. Similar restrictions can also be found in some semi-parametric factor models (e.g., Connor et al., Reference Connor, Hagmann and Linton2012; Fan et al., Reference Fan, Liao and Wang2016). However, such restrictions may be unrealistic in practice. As discussed in Connor et al. (Reference Connor, Hagmann and Linton2012), financial applications often require certain asset characteristics to vary over time. Their work demonstrates the time-series evolution of cross-sectional correlations among asset characteristics, indicating the need for time variation in asset attributes. Similarly, in macro-level data, characteristics often change over time, as they are aggregated from individual or household-level data. For instance, characteristics such as gender, race/ethnicity, or educational attainment are time-invariant at the individual level but vary over time after aggregation. Thus, the DP approach we consider is more adaptable to practical applications.

For practical implementation of $W_{it}$ , since the number of factors, J, is equal to the dimension of the diversified weight, we recommend a data-driven rank estimation procedure for choosing J. Specifically, we adopt the information criteria suggested by Su and Wang (Reference Su and Wang2017), which are tailored for time-varying factor structures. For the diversified weight functions $\varphi _j(\cdot )$ , we recommend using B-spline bases.Footnote ¹ Although standard sieve bases such as polynomial series, trigonometric series, and B-splines can all generate smooth functions, B-splines offer better numerical stability and more suitable boundary behavior. High-order polynomial series tend to induce severe multicollinearity as the sieve dimension J increases. Trigonometric bases impose a periodic boundary condition on the covariate support, which is typically unnatural for economic covariates. In contrast, B-spline basis functions only affect a local region of the covariate space, which provides stable boundary behavior and well-conditioned weight matrices. For an excellent review of the sieve basis function, we refer the readers to Chen (Reference Chen, Heckman and Leamer2007).

Table 1

Simulated coverage probability of 95% confidence interval for the standardized ATEL estimates and MSEs for different sieve basis functions $\varphi _j(\cdot )$ (B-splines, trigonometric series, and polynomial series) under DGP1–DGP3.

Table 1 Long description

The table consists of three main columns representing D G P 1, D G P 2, and D G P 3, each subdivided into Coverage and M S E.

* B-spline row: D G P 1 shows Coverage 0.9350 and M S E 0.2023. D G P 2 shows Coverage 0.9323 and M S E 0.0675. D G P 3 shows Coverage 0.9400 and M S E 0.0653.

* Trigonometric row: D G P 1 shows Coverage 0.9330 and M S E 0.2269. D G P 2 shows Coverage 0.9300 and M S E 0.0709. D G P 3 shows Coverage 0.9353 and M S E 0.0691.

* Polynomial row: D G P 1 shows Coverage 0.9327 and M S E 0.2055. D G P 2 shows Coverage 0.9373 and M S E 0.0681. D G P 3 shows Coverage 0.9407 and M S E 0.0655.

A note below the table specifies a sample size N equals 100, T sub 0 equals 50, T sub 1 equals 50, and the number of factors J equals 2, based on 3,000 repetitions.

Note: The sample size is $N=100$ , $T_0 = 50$ , $T_1= 50$ , and the number of factors $J=2$ . Based on simulations with 3,000 repetitions.

Table 1 examines the sensitivity of the ATEL estimator to the choice of sieve basis functions, $\varphi _j(\cdot )$ , used in the construction of $W_{it}$ . We evaluate this across different data-generating processes (DGPs) with the number of factors $J=2$ . Specifically, DGP1 specifies $h_t(\cdot )$ as a linear function, DGP2 uses a polynomial-series specification, and DGP3 uses a trigonometric-series specification. More details on these DGPs are provided in Section 4. Table 1 reports the simulation coverage probabilities of the 95% confidence intervals for the standardized ATEL estimates and the mean squared errors (MSEs). Overall, the simulation results are satisfactory. Across all DGPs, empirical coverage is close to the nominal level, and MSEs remain uniformly small. These findings suggest that our estimator is insensitive to the specific choice of sieve basis functions. Additionally, we present further simulations in Section 4 across different number of factors J for the construction of $W_{it}$ . Those results indicate that our ATEL estimates remain robust to the choice of J as well.

After estimating the factors, the next step is to estimate the time-varying factor loadings. The conventional PCA may fail under time-varying factor loadings. To see this, let $\tilde {F}_t$ be the PCA estimator of $F_t$ . Under regular conditions, the PCA estimator for the factor loading is

$$ \begin{align*} \tilde{\beta}_i = \frac{1}{T}\sum_{t=1}^T \tilde{F}_t Y_{it} = \frac{1}{T} \sum_{t=1}^T \tilde{F}_t F_t' \beta_{it} +o_P(1). \end{align*} $$

This representation shows that $\tilde {\beta }_i$ is a weighted average (over time) of $\beta _{it}$ . As a result, the PCA estimator captures only a global average factor loading rather than the spontaneous, time-varying parameter $\beta _{it}$ , and thereby failing to capture the variation in factor loadings over time.

To estimate time-varying factor loadings, we use local linear regression, a widely used nonparametric method that mitigates boundary bias (see, e.g., Fan, Reference Fan1996; Li and Racine, Reference Li and Racine2007). By assuming factor loading is a smooth function of $t/T$ , we have

$$ \begin{align*} \beta_{it} = \beta_i\left(\frac{t}{T}\right). \end{align*} $$

Since $\beta _i:[0,1]\rightarrow \mathbb {R}$ is a smooth function of $t/T$ , we can approximate $\beta _i(\frac {t}{T})$ by $\beta _i(\frac {r}{T})$ if $r/T$ is close to $t/T$ enough for $r\in \{1,\ldots ,T\}$ . Specifically, through Taylor expansion, we have

$$ \begin{align*} \beta_i(\frac{t}{T}) = \beta_i\left(\frac{r}{T}\right) + \beta_i^{(1)}\left(\frac{r}{T}\right)\left(\frac{t-r}{T}\right) + d_{i}(t,r), \end{align*} $$

where $\beta ^{(1)}_{i}(\cdot )$ is the first-order derivative of $\beta _i(\cdot )$ and $d_i(t,r)$ is the Taylor approximation error.

Based on this local linear approximation, the estimator for the time-varying factor loading is obtained by minimizing the locally weighted sum of squared residuals:

(10) $$ \begin{align} \min_{\beta_{1r},\beta^{(1)}_{1r}} \frac{1}{T_0} \sum_{t=1}^{T_0} \left[Y^N_{1t} - \beta_{1r}'\hat{F}_t - \beta^{(1)'}_{1r}\hat{F}_t \left(\frac{t-r}{T_0}\right) \right]^2 K_h\left(\frac{t-r}{T_0}\right), \text{ for }r = 1,\ldots,T_0, \end{align} $$

where $K_h(\cdot )$ is a kernel function and $h=h_{T_0}$ is the bandwidth parameter. To handle boundary issues and ensure uniform results, we adopt the kernel function as defined in Su and Wang (Reference Su and Wang2017).

Let $q_{tr}=\Big (\hat {F}^{\prime }_t,(\frac {t-r}{T_0})\hat {F}^{\prime }_t\Big )'$ be a $2J\times 1$ vector. Stacking these over $T_0$ periods, we define $Q_r = (q_{1r},\ldots ,q_{T_0r})'$ as a $T_0\times 2J$ matrix. Additionally, let $M_r = diag\Big (K_h(\frac {1-r}{T_0}),\ldots ,K_h(\frac {T_0-r}{T_0}) \Big )$ be a $T_0\times T_0$ matrix, and $Y_1 = (Y^N_{11},\ldots ,Y^N_{1T_0})'$ be a $T_0\times 1$ vector. Then the minimizer of (10) is given by

(11) $$ \begin{align} \begin{pmatrix} \hat{\beta}_{1r}\\ \hat{\beta}^{(1)}_{1r} \end{pmatrix} = (Q^{\prime}_rM_rQ_r)^{-1} Q_r'M_rY_1. \end{align} $$

The local linear estimator at the boundary point where $r=T_0$ is used for post-treatment periods. To better illustrate the idea, the Taylor expansion of $\beta _{1t}$ is

$$ \begin{align*} \beta_{1}\left(\frac{t}{T_1}\right) = \beta_{1}\left(\frac{T_0}{T_1}\right) + \beta_1^{(1)}\left(\frac{T_0}{T_1}\right)\left(\frac{t-T_0}{T_1}\right) + d_1(t,T_0). \end{align*} $$

This motivates our estimator of time-varying factor loading in the post-treatment periods is

(12) $$ \begin{align} \hat{\beta}_{1t} = \hat{\beta}_{1T_0} + \hat{\beta}_{1T_0}^{(1)}\left(\frac{t-T_0}{T_1}\right), \text{ for } t = T_0+1,\ldots,T. \end{align} $$

2.3 Comparison with Synthetic Control Method

The SCM, proposed by Abadie and Gardeazabal (Reference Abadie and Gardeazabal2003) and Abadie et al. (Reference Abadie, Diamond and Hainmueller2010), is a widely used approach for estimating the treatment effect. Given the structural similarities in how SCM and our DP approach construct estimators, it is important to compare these methods. SCM also conducts a cross-sectional weighted average of untreated outcomes to construct the estimator:

$$ \begin{align*} \hat{Y}^{SC}_{1t} =W^{SC'} Y_{t}, \end{align*} $$

where $W^{SC} = (W^{SC}_{2},\ldots ,W^{SC}_{N+1})$ and $Y_{t} = (Y^N_{2t},\ldots ,Y^N_{N+1,t})$ are $N\times 1$ vectors. The weights are selected by best fitting the outcome of the treated unit using pre-treatment data, and the weights are nonnegative and sum to one. Specifically, we focus on the case where all pre-treatment outcomes are included as predictors:

(13) $$ \begin{align} \min_{W \in \Delta } \sum_{t=1}^{T_0}\left( Y^N_{1t} - W' Y^{N}_{t} \right)^2, \end{align} $$

where $\Delta =\{W \in \mathbb {R}^N : W_i \geq 0 \ \forall i, \sum _{i=2}^{N+1} W_i = 1\}$ .

There are two fundamental differences between our DP approach and SCM. First, SCM uses only one weighted vector while we use the J weighting vectors. This difference arises because we use DP to estimate multiple factors, whereas SCM uses a weight to directly estimate the counterfactual outcomes. Second, SCM is primarily motivated by a time-invariant factor structure. To make the comparison transparent, we focus on the case where the untreated outcome follows a factor model with time-invariant factor loading with just one factor:

$$ \begin{align*} Y_{it}^N = \beta_i F_t+ u_{it},\quad i=1,\ldots,N+1. \end{align*} $$

Suppose we just use the SCM weight $W_i^{SC}$ to form our DP:

$$ \begin{align*} \hat F_t = \sum_{i=2}^{N+1}W_iY_{it}^N. \end{align*} $$

SCM uses $\hat F_t $ directly as the estimated counterfactual outcome:

$$ \begin{align*} \hat F_t = \hat{Y}^{SC}_{1t}, \end{align*} $$

while our DP approach uses $\hat F_t $ as the estimated factor.

Consequently, our approach requires an additional step to estimate the factor loadings of the treated unit by regressing the treated unit’s outcome on the estimated factors over the pre-treatment periods:

$$ \begin{align*} \hat\beta_1= \left(\sum_t \hat F_t^2\right)^{-1}\sum_t \hat F_tY_{1t}^N. \end{align*} $$

Abadie et al. (Reference Abadie, Diamond and Hainmueller2010) provide conditions under which the existence of $W^{SC}$ requires the factor loadings of the treated unit to be perfectly reconstructed via a convex combination of the factor loadings in the control group; that is,

(14) $$ \begin{align} \beta_1 = \sum_{i=2}^{N+1} W^{SC}_i \beta_i. \end{align} $$

In this case,

$$ \begin{align*} \begin{aligned} \hat F_t &=\sum_{i=2}^{N+1}W_i \beta_i F_t+ \sum_{i=2}^{N+1}W_i u_{it}\\ &=\beta_1 F_t+ o_P(1). \end{aligned} \end{align*} $$

This shows that, under the SCM weight condition, the DP estimated factors $\hat F_t$ indeed can be used as the counterfactual estimator, which is consistent for $\beta _1 F_t.$

Meanwhile, we use $\hat F_t$ as the factor estimator and apply one further step to estimate $\beta _1$ . Under the SCM condition (14), our approach is also consistent. To see this, note that

$$ \begin{align*} \begin{aligned} \hat\beta_1 &=\left(\sum_t \hat F_t^2\right)^{-1}\sum_t \hat F_t (F_t\beta_1+u_{1t}) \\&=\left(\sum_t \hat F_t^2\right)^{-1}\sum_t \sum_{i=2}^{N+1}\hat F_tF_t\beta_i W_i +\left(\sum_t \hat F_t^2\right)^{-1}\sum_t \hat F_t u_{1t} \\&=\left(\sum_t \hat F_t^2\right)^{-1}\sum_t \sum_{i=2}^{N+1}\hat F_t(Y_{it}^N-u_{it})W_i +\left(\sum_t \hat F_t^2\right)^{-1}\sum_t \hat F_t u_{1t} \\&=1+\zeta=1+o_P(1), \end{aligned} \end{align*} $$

where $\zeta = (\sum _t \hat F_t^2)^{-1}\sum _t \hat F_t u_{1t}-(\sum _t \hat F_t^2)^{-1}\sum _t \sum _{i=2}^{N+1}\hat F_tu_{it}W_i.$

As a result,

$$ \begin{align*} \hat\beta_1 \hat F_t= \hat F_t+o_P(1) = \beta_1 F_t+ o_P(1). \end{align*} $$

Therefore, if the treated unit’s factor loading lies within the convex hull of the control group’s factor loadings, both methods are consistent. However, when this condition fails, SCM is not consistent. SCM requires the weights to be well-chosen to satisfy the exact matching condition of factor loadings. In contrast, because we flexibly estimate $\beta _1$ , our framework does not require the weights to be precisely chosen. Indeed, our requirement for W is weaker than SCM.

Finally, we acknowledge that SCM is simpler to implement which only requires one set of weights. The DP, however, requires choosing the number of factors and choosing multiple weights to satisfy our conditions.

2.4 Estimation Procedure

With low-rank approximations, untreated potential outcomes can be modeled as a factor model. Following the factor-based approach proposed by Gobillon and Magnac (Reference Gobillon and Magnac2016) and Li and Sonnier (Reference Li and Sonnier2023), we separately estimate the factors and factor loadings (the treated unit). The detailed estimation steps are outlined below:

1. Step 1: (Initialization) Pre-determine the diversified weights $W_{it}$ as described in Section 2.2.

2. Step 2: (Factors $F_t$ ) Using the DP method, the factor $F_t$ is estimated by projecting the untreated outcomes $Y^N_{it}$ of the control group onto the space defined by the diversified weights $W_{it}$ :

   (15) $$ \begin{align} \hat{F}_t=\frac{1}{N}W_t'Y_t, \text{ for } t = 1,\ldots,T, \end{align} $$
   where $W_t = (W_{2t},\ldots ,W_{N+1,t})'$ is the $N\times J$ matrix.

3. Step 3: (Factor Loading $\beta _{1t}$ ) Run local linear estimator $Y^N_{1t}$ on $[\hat {F}_t, \hat {F}_t (\frac {t-T_0}{T_0})]$ with pre-treatment periods, obtain $\hat {\beta }_{1T_0}$ and $\hat {\beta }^{(1)}_{1T_0}$ , then

   $$ \begin{align*} \hat{\beta}_{1t} = \hat{\beta}_{1T_0} + \hat{\beta}^{(1)}_{1T_0}\left(\frac{t-T_0}{T_1}\right), \text{ for }t=T_0+1,\ldots,T. \end{align*} $$

4. Step 4: (Counterfactual Outcome) With the estimated factors and factor loading, the counterfactual outcomes for the treated unit $i=1$ during the post-treatment period $t=T_0+1,\ldots ,T$ are constructed as

   (16) $$ \begin{align} \hat{Y}^N_{1t} = \hat{\beta}^{\prime}_{1t} \hat{F}_t, \text{ for } t = T_0+1,\ldots,T. \end{align} $$

5. Step 5: (ATEL Estimator) Using the estimated counterfactual outcomes, the treatment effect for each post-treatment period $(t=T_0+1,\ldots ,T)$ is given by

   (17) $$ \begin{align} \hat{\alpha}_{1t} = Y_{1t} - \hat{Y}^N_{1t}. \end{align} $$

   The ATEL estimator is then computed as

   (18) $$ \begin{align} \hat{\alpha} = \frac{1}{T_1} \hat{\alpha}_{1}' \mathbf{K}_{T_1}, \end{align} $$
   where $\hat {\alpha }_{1} = (\hat {\alpha }_{1,T_0+1},\ldots ,\hat {\alpha }_{1T})'$ and $\mathbf {K}_{T_1} = \Big (K_h(\frac {(T_0+1)-T_0}{T_1}),\ldots ,K_h(\frac {T-T_0}{T_1})\Big )'$ are $T_1\times 1$ vector.

Our estimation procedure is computationally efficient and scales well to large panels. It consists of two main steps: (i) factor estimation via DP, a cross-sectional weighted average of untreated outcomes, which can be computationally less costly than PCA because it avoids eigen-decomposition for large matrices and (ii) factor-loading estimation via local linear smoothing (weighted linear regression), which is a standard procedure and remains easy to compute even if $T_0$ is large, since it uses only the $\lfloor T_0h \rfloor $ observations in a neighborhood of the boundary point $T_0$ . In practice, both steps can be implemented efficiently via parallelization over t.Footnote ²

3.1 Inferential Theory

This section presents the inferential theory. We provide the asymptotic normality of the ATEL. Remind the following notation:

$$ \begin{align*} h_t(\eta_{it}, X_{it}) = \sum_{j=1}^J \phi_j(\eta_{it},X_{it}) f_{jt}+ r_{it} = \beta_{it}'F_t + r_{it}. \end{align*} $$

Let $W_t = (W_{2t},\ldots , W_{N+1,t})'$ be the $N \times J$ matrix. Each of its rows $W_{it}=(w_{it1},\ldots ,w_{itJ})'$ is a $J\times 1$ vector for $i=2,\ldots ,N+1$ and $t=1,\ldots ,T$ . $\beta _{t} = (\beta _{2t},\ldots ,\beta _{N+1,t})'$ is the $N \times J$ matrix. $F_t = (f_{t1},\ldots ,f_{tJ})'$ is a $J\times 1$ vector. $u_t = (u_{1t},\ldots ,u_{N+1,t})'$ is $(N+1) \times 1$ vectors for $t = 1,\ldots ,T$ .

We now state the conditions and asymptotic results of estimators.

Assumption 1 presents some conditions for sieve approximation. Assumption 1(i) restricts $\{h_t(\cdot )\}_{t\leq T}$ to a Hilbert space with a uniform $L_2$ bound such that $\max _{t\leq T} \sum _{j=1}^\infty f^2_{tj}< \infty $ . Moreover, $F_t$ is nonrandom since it’s sieve coefficients. Then we illustrate that, under some sensible conditions on the functional space and the sieve bases, Assumption 1(ii) is well satisfied. Suppose $h_t$ belongs to a Hölder class: For some $b,C,d> 0$ , $\{ h: \| D^b h(x) - D^b h(y)\| \leq C \|x - y\|^{d}, \text {for all } x,y\}$ . Moreover, suppose the sieve bases follow some common basis functions, such as B-spline, trigonometric polynomials, or polynomial series. Then, with the sufficiently smooth functions $\{h_t(\cdot ) \}$ , sieve approximation error will satisfy $\max _{it} |r_{it}| < C J^{-a}$ , $a = (b+d) / (dim(\eta _{it}) + dim(X_{it}))$ . Assumption 1(iii) can be satisfied if the basis functions are bounded such as trigonometric basis or the support of unobserved and observed covariates $\eta _{it}$ , $X_{it}$ is bounded. It implies that $\| \beta _{it}\|^2 = O_P(J)$ since factor loading is a vector of sieve basis function and the number of sieve terms grows slowly as $N,T \rightarrow \infty .$

Assumption 2 ensures that the diversified idiosyncratic error, $\frac {1}{N}\sum _{i=2}^{N+1} W_{it}u_{it}$ , satisfies the conditions for a cross-sectional central limit theorem when the idiosyncratic errors ( $u_{2t},\ldots ,u_{N+1,t}$ ) are cross-sectionally weakly dependence. Furthermore, according to Assumption 1(ii), the diversified sieve approximation error $\frac {1}{N}\sum _{i=2}^{N+1}W_{it} r_{it}$ is negligible than the diversified idiosyncratic error such that

$$ \begin{align*} \hat{F}_t - H_tF_t = \frac{1}{N}\sum_{i=2}^{N+1} W_{it} u_{it} + o_p(1). \end{align*} $$

Therefore, the $\sqrt {N} (\hat {F}_t - H_tF_t)$ is asymptotically normal when the idiosyncratic errors ( $u_{2t},\ldots ,u_{N+1,t}$ ) are cross-sectionally weakly dependence for each $t \leq T$ . In particular, this asymptotic normality result still holds even in the finite T, which differs from conventional PCA requiring both $N, T \rightarrow \infty $ .

Assumption 3 imposes conditions on factors and factor loadings, similar to those in Su and Wang (Reference Su and Wang2017). The key difference is that they assume random factors and nonrandom factor loadings, whereas we assume nonrandom factors and random factor loadings, derived from the sieve approximation. Assumption 3(i) ensures the nonsingularity of $\frac {1}{T}\sum _{t=1}^T F_t F_t'$ for sufficiently large T. Assumption 3(ii) imposes pervasive condition on factor loadings.

Assumption 4 imposes the common smoothness restriction on time-varying factor loadings typically assumed in the existing literature. The main difference is that, while many studies directly assume that factor loadings are functions of $t/T$ , we approximate the unknown function $h_t(\eta _{it},X_{it})$ via a sieve approach, where the factor loadings are represented by sieve basis functions of the covariates $\eta _{it}$ and $X_{it}$ . The assumption is plausible because aggregate variables often evolve slowly. In macroeconomic applications, variables such as GDP, inflation, and interest rates are commonly assumed to have trends or low-frequency components that evolve smoothly (see, e.g., Nelson and Plosser, Reference Nelson and Plosser1982; Hodrick and Prescott, Reference Hodrick and Prescott1997).

When $\eta _{it}$ and $X_{it}$ follow smooth trajectories over time, that is, $\eta _{it} = \eta _i(\frac {t}{T})$ and $X_{it} = X_i(\frac {t}{T})$ , we have

$$ \begin{align*} \beta_{it} = \phi(\eta_{it},X_{it}) = \phi\left(\eta_i\left(\frac{t}{T}\right),X_i\left(\frac{t}{T}\right) \right) = \beta_i\left(\frac{t}{T}\right). \end{align*} $$

Likewise, the diversified weight $W_{it}$ is also a sieve basis function of $X_{it}$ , it follows that $W_{it}$ is a smooth function of time. Accordingly, Assumption 4 implies that $H_t = \frac {1}{N}\sum _{i=2}^{N+1}W_{it}\beta _{it}$ is continuous with respect to t.

We impose the $\alpha $ -mixing condition. Let $\mathcal {F}^0_{-\infty }$ and $\mathcal {F}^{\infty }_T$ denote the $\sigma $ -algebra generated by $\{ u_t:t\leq 0\}$ and $\{u_t: t \geq T\},$ respectively. Define the mixing coefficient

$$ \begin{align*}\alpha(T) = sup_{A\in \mathcal{F}^0_{\infty}, B\in\mathcal{F}^{\infty}_{T} } |P(A)P(B) -P(AB) |.\end{align*} $$

Assumption 5(i) and (ii) is a standard requirement for moments and the mixing coefficient for $\alpha $ -mixing time series. In particular, since diversified weight $W_{it}$ is chosen based on $X_{it}$ , the independence condition in Assumption 5(i) implies Assumption 2(iii). Assumption 5(iii) allows the idiosyncratic error $u_{it}$ to be weak dependence on both time series and cross section, which are in the same spirit of approximate factor analysis (see, e.g., Bai and Ng, Reference Bai and Ng2002; Bai, Reference Bai2003).

Assumption 6(i) is employed to bound the approximation error of the sieve approximation, requiring $\{h_t(\cdot )\}_{t\leq T}$ to be sufficiently smooth to satisfy this condition.

Local linear regression is applied to estimate $\hat {\beta }_{it}, \hat {\beta }_{it}^{(1)}$ . Standard practice involves imposing mild conditions on the kernel functions and the convergence properties of the bandwidth, h.

Assumption 7 is the standard assumption for local linear regression. Several types of kernel function satisfy Assumption 7(i), such as Epanechnikov, uniform, and quartic kernels. To handle boundary issues, we define the kernel function $K(\cdot )$ based on $K^*(\cdot )$ following Su and Wang (Reference Su and Wang2017):

(19) $$ \begin{align} \begin{aligned} K_h(\frac{t-r}{T})= \frac{1}{h}K(\frac{t-r}{Th}) = \begin{cases} \frac{1}{h} K^*(\frac{t-r}{Th}) / \int^1_{-r/(Th)} K^*(u)du, \quad\text{ if }r\in [1, \lfloor Th \rfloor )\\ \frac{1}{h} K^*(\frac{t-r}{Th}), \quad\text{ if }r\in [\lfloor Th \rfloor , T-\lfloor Th \rfloor ]\\ \frac{1}{h} K^*(\frac{t-r}{Th}) / \int^{(1-r/T)/h}_{-1} K^*(u)du, \quad\text{ if }r\in (T-\lfloor Th \rfloor, T ].\\ \end{cases} \end{aligned} \end{align} $$

To establish an asymptotic result, define $\mu _j=\int u^j K(u)du$ and $\nu _j = \int u^jK^2(u)du$ . Specifically, $\mu _0=1$ uniformly.

Formally, the following theorem shows that $\hat {\alpha }$ is asymptotic normal.

Theorem 1. With Assumptions 1–7, as $N,T_0,T_1\rightarrow \infty $ ,

$$ \begin{align*} (\frac{1}{T_0h} \Sigma_1 + \frac{1}{T_1h} \Sigma_2)^{\frac{-1}{2}}(\hat{\alpha} - \alpha) \xrightarrow{d} N(0, 1), \end{align*} $$

where $\Sigma _1 {\kern-1pt}={\kern-1pt} \nu _0 \lambda '\xi ^{-1}\Omega \xi ^{-1}\lambda $ , $\lambda {\kern-1pt}={\kern-1pt} \lim _{T_1 \rightarrow \infty }\frac {1}{T_1}{\kern-1pt}\sum _{t=T_0+1}^T F_t$ , $\xi {\kern-1pt}={\kern-1pt} \lim _{T_0 \rightarrow \infty } \frac {1}{T_0}{\kern-1pt}\sum _{t=1}^{T_0}{\kern-1pt} F_tF_t' $ , $\Omega = \lim _{T_0 \rightarrow \infty } \left [ \frac {1}{T_0}\sum _{t=1}^{T_0} F_tF^{\prime }_t \mathbb {E}(u^2_{1t}) + \frac {2}{T_0}\sum _{t=1}^{T_0}\sum _{s=t+1}^{T_0}F_tF^{\prime }_s \mathbb {E}(u_{1t} u_{1s})\right ]$ , $\Sigma _2 = \nu _0 \sigma ^2$ , $\sigma ^2 =\lim _{T_1\rightarrow \infty } \left [ \frac {1}{T_1} \sum _{t=T_0+1}^T \mathbb {E}(u_{1t}^2) + \frac {2}{T_1}\sum _{t=T_0+1}^T \sum _{s=t+1}^T \mathbb {E}(u_{1t} u_{1s})\right ]$ .

The asymptotic result holds when N, $T_0$ , and $T_1$ are large, corresponding to the factor-based estimation of the average treatment effect (e.g., Bai and Ng, Reference Bai and Ng2021; Li and Sonnier, Reference Li and Sonnier2023). The asymptotic variance comprises two terms. Because the two variance terms vanish at different rates, the rate of convergence of $\hat {\alpha }$ is $1/\delta _{T_0,T_1,h}$ , where $\delta _{T_0,T_1,h} = \min \{\sqrt {T_0h},\sqrt {T_1h} \}$ . It is important to note that the ATEL estimator is always asymptotically normal regardless of the limit relationship between $T_0$ and $T_1$ . To see this, the asymptotic distribution can be rewritten as

$$ \begin{align*} \begin{aligned} \left(\frac{\delta^2_{T_0,T_1,h}}{T_0h}\Sigma_1 + \frac{\delta^2_{T_0,T_1,h}}{T_1h} \Sigma_2 \right)^{-\frac{1}{2}} (\delta_{T_0,T_1,h} (\hat{\alpha} - \alpha) ) \xrightarrow{d} N(0, 1). \end{aligned} \end{align*} $$

In particular, if $T_0/T_1\rightarrow 0$ , $\sqrt {T_0h}(\hat {\alpha } - \alpha ) \xrightarrow {d} N(0, \Sigma _1)$ ; if $T_1/T_0 \rightarrow 0$ , $ \sqrt {T_1h}(\hat {\alpha } - \alpha ) \xrightarrow {d} N(0, \Sigma _2)$ .

The validity of the asymptotic inference does not require the consistent estimation of $H_t$ because, similar to the standard time-invariant factor model, the rotation matrix is not identifiable. This issue is analogous to the rotation problem in the usual PCA setting, where factors are estimated up to an unknown rotation H that itself cannot be consistently estimated. However, this lack of identification does not affect the validity of asymptotic inference, because the asymptotic variance would depend on H only through either $H F_t$ or $H^{-1} \beta _{i}$ . In either case, we can replace these terms with their consistent estimators, $\hat {F}_t$ and $\hat {\beta }_{i}$ , to consistently estimate the variance. This analogy carries over to our setting where the rotation matrix $H_t$ varies over time. In other words, the asymptotic variance is rotation-free and can be consistently estimated, even though the rotation matrix $H_t$ itself is not.

We now turn to the consistent estimation of the asymptotic covariance in Theorem 1. Given the weak dependence of the time series on the idiosyncratic error, an estimator of $\Sigma _1$ is given by

$$ \begin{align*} \hat{\Sigma}_1 = \hat{\lambda}'\hat{\xi}^{-1}\hat{\Omega} \hat{\xi}^{-1}\hat{\lambda,} \end{align*} $$

where $\hat {\lambda } = \frac {1}{T_1\mu _0}\sum _{t=T_0+1}^{T} \hat {F}_t K_h(\frac {t-T_0}{T_1})$ , $\hat {\xi }= \frac {1}{T_0 \mu _0 } \sum _{t=1}^{T_0} \hat {F_t}\hat {F_t}' K_h(\frac {t-T_0}{T_0}) $ , and

$$ \begin{align*} \hat{\Omega} = \frac{h}{T_0} \sum_{t=1}^{T_0} \hat{u}^2_{1t} \hat{F_t} \hat{F_t}' K^2_h\left(\frac{t-T_0}{T_0}\right) + \frac{2h}{T_0} \sum_{j=1}^{T_0-1} \sum_{t = j+1}^{T_0} \hat{u}_{1t} \hat{u}_{1j} \hat{F}_t \hat{F}^{\prime}_{j} K_h\left(\frac{t-T_0}{T_0}\right)K_h\left(\frac{j-T_0}{T_0}\right), \end{align*} $$

where $\hat {u}_{1t} = Y_{1t} - \hat {\beta }^{\prime }_{1t} \hat {F}_t$ .

An estimator of $\Sigma _2$ is given by

$$ \begin{align*} \hat{\Sigma}_2 = \frac{h}{T_0} \sum_{t=1}^{T_0} \hat{u}^2_{1t} K^2_h\left(\frac{t-T_0}{T_0}\right) + \frac{2h}{T_0} \sum_{j=1}^{T_0-1} \sum_{t=j+1}^{T_0} \hat{u}_{1t} \hat{u}_{1j}K_h\left(\frac{t-T_0}{T_0}\right)K_h\left(\frac{j-T_0}{T_0}\right). \end{align*} $$

If $u_{1t}$ is serially uncorrelated, $\hat {\Sigma }_1$ and $\hat {\Sigma }_2$ can be simplified:

$$ \begin{align*} \hat{\Sigma}_1 = \ \hat{\lambda}'\hat{\xi}^{-1}\hat{V} \hat{\xi}^{-1} \hat{\lambda},\quad \hat{\Sigma}_2 = \frac{h}{T_0} \sum_{t=1}^{T_0} \hat{u}^2_{1t} K^2_h(\frac{t-T_0}{T_0}), \end{align*} $$

where $\hat {V} = \frac {h}{T_0} \sum _{t=1}^{T_0} \hat {u}^2_{1t} \hat {F_t} \hat {F_t}' K^2_h(\frac {t-T_0}{T_0})$

In practice, a data-driven approach is recommended for bandwidth selection. One effective method is leave-one-out cross-validation, which determines the optimal bandwidth, $h_{cv}$ , by minimizing the following objective function:

(20) $$ \begin{align} h_{cv} = \arg\min_{h} \sum_{t=1}^{T_0}( y_{1t} - \hat{\beta}^{(-t)}_{1t} \hat{F}_t )^2, \end{align} $$

where $\hat {\beta }^{(-t)}_{1t}$ is the leave-one-out estimator. Specifically, $\hat {\beta }^{(-t)}_{1t}$ is the counterpart of $\hat {\beta }_{1t} = \hat {\beta }_{1T_0} + \hat {\beta }^{(1)}_{1T_0}(\frac {t-T_0}{T_0})$ , but excludes the t-th observation in the local linear estimation process.

While our main focus is inference for the ATEL, practitioners may also be interested in inference on the estimated counterfactual at each post-treatment period by constructing pointwise confidence intervals. This requires further assumption on $\beta _{it}$ . Consider the estimated counterfactual outcomes at $t>T_0$ :

$$ \begin{align*} \begin{aligned} \hat{Y}^N_{1t} &= -\underbrace{\left(\hat{\beta}_{1T_0}-H^{-1}_{T_0}\beta_{1T_0}\right)'H_tF_t}_{O_P\left(\frac{1}{\sqrt{T_0h}}\right)} - \underbrace{\left(\hat{\beta}^{(1)}_{1T_0}-H_{T_0}^{-1}\beta_{1T_0}^{(1)}\right)'H_tF_t\left(\frac{t-T_0}{T_1}\right)}_{O_P\left(\frac{1}{\sqrt{T_0h}}\right)}+\underbrace{\beta^{\prime}_{1t}F_t}_{O_P(1)}+o_P(1). \end{aligned} \end{align*} $$

The first two terms, estimation error in the factor loading for the pre-treatment periods, are the source of $\Sigma _1$ (before averaging over post-treatment periods), while the third term, $\beta _{1t}'F_t$ , is the common components of the factor model. Since we are not considering the treatment effect, the common component $\beta _{1t}'F_t$ is dominant here. If we further assume that $\beta _{1t}$ are i.i.d. and follow a normal distribution, the statistical inference can still be implemented. For example, suppose $\beta _{1t} \sim N(\bar {\beta }, \Sigma _{\beta } )$ with $\Sigma _{\beta } = diag\Big ( \sigma ^2_{1},\ldots ,\sigma ^2_{J} \Big )$ . Then, for each post-treatment period, the variance of the estimated counterfactual outcome is

$$ \begin{align*} \begin{aligned} \tilde{\Sigma}_t = \tilde{\Sigma}_{1,t} + \tilde{\Sigma}_{2,t}, \end{aligned} \end{align*} $$

where $\tilde {\Sigma }_{1,t}=\frac {1}{T_0h}\nu _0 F_t' \xi ^{-1}\Omega \xi ^{-1}F_t$ and $\tilde {\Sigma }_{2,t} = F_t' \Sigma _{\beta } F_t$ .

The estimators for $\tilde {\Sigma }_{1,t}$ and $\tilde {\Sigma }_{2,t}$ are given by

$$ \begin{align*} \hat{\tilde{\Sigma}}_{1,t} = \frac{1}{T_0h} \hat{F}_t \hat{\xi}^{-1}\hat{\Omega}\hat{\xi}^{-1}\hat{F}_t, \quad \hat{\tilde{\Sigma}}_{2,t} = \hat{F}^{\prime}_t \hat{\Sigma}_{\beta}\hat{F}_t, \end{align*} $$

where $\hat {\Sigma }_{\beta } = \frac {1}{T_1}\sum _{t=T_0+1}^{T} ( \hat {\beta }_{1t}-\hat {\bar {\beta }}) ( \hat {\beta }_{1t}-\hat {\bar {\beta }})'K_h(\frac {t-T_0}{T_1})$ and $\hat {\bar {\beta }} = \frac {1}{T_1}\sum _{t=T_0+1}^{T} \hat {\beta }_{1t} K_h(\frac {t-T_0}{T_1})$ . Accordingly, the pointwise confidence intervals for the estimated counterfactual outcome can be constructed as

(21) $$ \begin{align} \begin{aligned} \hat{Y}^N_{1t}\pm c_{\alpha} \hat{\tilde{\Sigma}}^{1/2}_t, \text{ for }t=T_0+1,\ldots,T, \end{aligned} \end{align} $$

where $c_{\alpha }$ is the critical value of a standard normal or t-distribution for the given confidence interval, and $\hat {\tilde {\Sigma }}_t$ is the estimator for $\tilde {\Sigma }_t$ .

3.2 Extensions

We conclude this section with a few extensions of our estimator that highlight the flexibility of our approach.

3.2.1 Dynamic ATEL

The original ATEL parameter focuses on the short-term effect localized around the intervention date. However, researchers are often also interested in the dynamic treatment effect, that is, how the policy impact evolves after the initial intervention (e.g., Callaway and Sant’Anna, Reference Callaway and Sant’Anna2021; Sun and Abraham, Reference Sun and Abraham2021). In particular, one may ask: How does the effect of treatment vary with the length of exposure? For example, does ATEL increase or decrease as more time elapses after treatment begins?

The ATEL framework can be naturally generalized to a dynamic version that localizes the treatment effect at later post-treatment periods. Specifically, for a given post-treatment horizon m, define

$$ \begin{align*} \alpha_m = \frac{1}{T_1}\sum_{t=T_0+1}^T \mathbb{E}(Y^I_{1t}-Y^N_{1t}) K_h\left(\frac{t-(T_0+m)}{T_1}\right), \end{align*} $$

which denotes the dynamic ATEL, capturing the short-term effect localized around the post-treatment $T_0+m$ . The corresponding estimator is constructed analogously as

$$ \begin{align*} \hat{\alpha}_m = \frac{1}{T_1}\sum_{t=T_0+1}^T \hat{\alpha}_{1t} K_h\left(\frac{t-(T_0+m)}{T_1}\right). \end{align*} $$

The asymptotic normality result in Theorem 1 remains valid for $\hat {\alpha }_m$ , provided that m is fixed or $m/T_1h \rightarrow 0$ .

3.2.2 ATEL with Carryover Effects

In many applications, treatment effects may carry over in time (Imai et al., Reference Imai, Kim and Wang2023). Treatment received in earlier periods can continue to influence outcomes in later periods, even if the current treatment status is the same. In other words, the outcome at time t may depend not only on whether the unit is treated at t, but also on how long it has been treated and on its past treatment history.

We extend the ATEL framework to accommodate such carryover effects. Let the binary treatment path up to time t be denoted by $D_{1:t}\in \{0,1\}^t$ . In the potential outcomes framework, the potential outcome is written as $Y_{1t}(D_{1:t})$ , representing the outcome that the treated unit would experience at time t under the treatment path $D_{1:t}$

Following Imai et al. (Reference Imai, Kim and Wang2023), we impose a limited (finite) carryover assumption that the realization of the potential outcome depends only on the most recent L lags of the treatment path.

Assumption 8 (Finite Carryover)

For all treatment paths $D_{1:t}\in \{0,1\}^t$ ,

$$ \begin{align*} Y_{1t}(D_{1:t}) = Y_{1t}\left(D_{t-L+1},\ldots,D_{t}\right). \end{align*} $$

This assumption explicitly restricts carryover effects to at most L periods. Since ATEL measures a local treatment effect using a narrow window in the post-treatment period, it can be generalized to this finite carryover setting. The estimand of interest is now defined as

$$ \begin{align*} \alpha(L) = \frac{1}{T_1}\sum_{t=T_0+L+1}^T \mathbb{E}(Y_{1t}(\boldsymbol{1})-Y_{1t}(\boldsymbol{0})) K_h\left(\frac{t-(T_0+L)}{T_1}\right), \end{align*} $$

where $\boldsymbol {1}$ and $\boldsymbol {0}$ denote L-dimensional vectors of ones and zeros, respectively. Thus, $Y_{1t}(\boldsymbol {1})$ and $Y_{1t}(\boldsymbol {0})$ are the potential outcomes for the treated unit that has been continuously treated or untreated over the last L periods prior to t.

The key point is that we exclude the first L post-treatment periods, $t=T_0+1,\ldots ,T_0+L$ , because the corresponding potential outcomes $Y_{1t}(\boldsymbol {1})$ involve a counterfactual history in which the unit would have already been treated for L periods, which is incompatible with the observed treatment path immediately after $T_0$ . We therefore construct the ATEL estimator with carryover effects as

$$ \begin{align*} \hat{\alpha}(L) = \frac{1}{T_1}\sum_{t=T_0+L+1}^T \hat{\alpha}_{1t} K_h(\frac{t-(T_0+L)}{T_1}). \end{align*} $$

This estimator is a kernel-weighted average of the post-treatment effects after the policy has been implemented for at least L periods. The asymptotic normality result in Theorem 1 remains valid for $\hat {\alpha }(L)$ , since trimming a fixed number L of post-treatment periods is asymptotically negligible.