2.1 The Model

The scalar form of our model is given by

(2.1) $$ \begin{align} y_{it}=(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0})^{\prime }x_{it}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it},\quad i\in \lbrack N],~t\in \lbrack T], \end{align} $$

where $x_{it}$ , $f_{t}^{0}$ , $\lambda _{i}^{0}$ , and $u_{it}$ are defined as above. Note that $\beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}$ are not separately identified. In other words, equation (2.1) cannot be directly employed to construct a least-squares estimator. To achieve identification, we assume $E(\beta _{i\cdot }^{0})=0$ and leave $ E(\beta _{\cdot t}^{0})$ unrestricted. Define $\dot {\beta }_{i\cdot }=\beta _{i\cdot }^{0}-\frac {1}{N}\sum _{j=1}^{N}\beta _{j\cdot }^{0}$ and $\dot {\beta }_{\cdot t}=\beta _{\cdot t}^{0}+\frac {1}{N}\sum _{j=1}^{N}\beta _{j\cdot }^{0}$ . By construction, $\sum _{j=1}^{N}\dot {\beta }_{j\cdot }=0$ . We can then rewrite the model as

(2.2) $$ \begin{align} y_{it} &=(\dot{\beta}_{i\cdot }+\dot{\beta}_{\cdot t})^{\prime }x_{it}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it},\quad i\in \lbrack N-1],~t\in \lbrack T], \end{align} $$

(2.3) $$ \begin{align} y_{Nt} &=\Big(-\sum_{j=1}^{N-1}\dot{\beta}_{j\cdot }+\dot{\beta}_{\cdot t} \Big)^{\prime }x_{Nt}+\lambda _{N}^{0\prime }f_{t}^{0}+u_{Nt},\quad t\in \lbrack T]. \end{align} $$

Let $\dot {B}^{(1)}\equiv \lbrack \dot {\beta }_{1\cdot }^{\prime },\ldots , \dot {\beta }_{N-1\cdot }^{\prime }]^{\prime }\in \mathbb {R}^{(N-1)K}$ and $ \dot {B}^{(2)}\equiv \lbrack \dot {\beta }_{\cdot 1}^{\prime },\ldots ,\dot { \beta }_{\cdot T}^{\prime }]^{\prime }\in \mathbb {R}^{TK}$ . Equations (2.2) and (2.3) can be employed to construct a least squares estimator.

We establish two alternative representations of models (2.2) and (2.3). The vector representation facilitates the analysis of the least-squares estimators for elements in $\dot {B}\equiv \lbrack \dot {B} ^{(1)\prime },\dot {B}^{(2)\prime }]^{\prime }$ in the following discussion, while the matrix representation is particularly useful for studying the estimator associated with the factor structure.

To proceed, we introduce some notation. Let $\mathbf {y}_{i\cdot }\equiv \lbrack y_{i1},y_{i2},\ldots ,y_{iT}]^{\prime } \in \mathbb {R}^{T}$ , $ \mathbf {y}_{\cdot t}\equiv \lbrack y_{1t},y_{2t},\ldots ,y_{Nt}]^{\prime } \in \mathbb {R}^{N}$ , and $Y\equiv \lbrack \mathbf {y}_{1\cdot }^{\prime }, \mathbf {y}_{2\cdot }^{\prime },\ldots ,\mathbf {y}_{N\cdot }^{\prime }]^{\prime }\in \mathbb {R}^{NT}$ . Define $\mathbf {u}_{i\cdot },\mathbf {u} _{\cdot t},$ and U analogously. Let $\mathbf {x}_{i\cdot }\equiv \lbrack x_{i1},x_{i2},\ldots ,x_{iT}]^{\prime }\in \mathbb {R}^{T\times K}$ and $ \mathbf {x}_{\cdot t}\equiv \lbrack x_{1t},x_{2t},\ldots ,x_{Nt}]^{\prime }\in \mathbb {R}^{N\times K}$ be the covariate matrices that may appear in equations for individual i and time t, respectively. Define

$$ \begin{align*} X^{(1)}\equiv \underbrace{ \begin{bmatrix} \mathbf{x}_{1\cdot } & & \\ & \ddots & \\ & & \mathbf{x}_{N-1\cdot } \\ -\mathbf{x}_{N\cdot } & \cdots & -\mathbf{x}_{N\cdot } \end{bmatrix} }_{NT\times (N-1)K},\text{ }\mathbf{x}_{i\cdot }^{(2)}\equiv \underbrace{ \left[ \begin{array}{ccc} x_{i1}^{\prime } & & \\ & \ddots & \\ & & x_{iT}^{\prime } \end{array} \right] }_{T\times TK}\text{, and }X^{(2)}\equiv \underbrace{\left[ \begin{array}{c} \mathbf{x}_{1\cdot }^{(2)} \\ \mathbf{x}_{2\cdot }^{(2)} \\ \vdots \\ \mathbf{x}_{N\cdot }^{(2)} \end{array} \right] }_{NT\times TK}. \end{align*} $$

Let $X^{(1,2)}\equiv \lbrack X^{(1)},X^{(2)}]$ . Given the above notations, the model in (2.2) and (2.3) can be rewritten in an $NT$ -vector form:

(2.4) $$ \begin{align} \mathbf{Vector}\text{ }\mathbf{representation}:\quad Y=X^{(1)}\dot{B} ^{(1)}+X^{(2)}\dot{B}^{(2)}+\theta ^{0}+U=X^{(1,2)}\dot{B}+\theta ^{0}+U, \end{align} $$

where $\theta ^{0}\equiv $ vec $(\Theta ^{0})$ . The vector $\dot { B}$ stacks all free unit- and time-specific slope coefficients into a single $(N+T-1)K\times 1$ vector. Accordingly, for each pair $(i,t)$ , the corresponding row of $X^{(1,2)}$ consists of an expanded regressor vector, constructed by interacting $x_{it}$ with the full set of unit and time dummies, while imposing the normalization constraint $\sum _{i=1}^{N}\dot { \beta }_{i\cdot }=0$ . Consequently, the product $X^{(1,2)}\dot {B}$ is equivalent to a regression specification that includes all unit–time interactive effects involving $x_{it}$ . Next, we introduce the matrix representation. Denote the kth entries of $x_{it}, \dot {\beta }_{i\cdot }$ , and $\dot {\beta }_{\cdot t}$ as $x_{k,it}$ , $\dot {\beta }_{k,i\cdot }$ , and $ \dot {\beta }_{k,\cdot t},$ respectively. Let $\mathbf {X}_{k}$ be a $T\times N$ matrix with $(t,i)$ th entry $x_{k,it}$ , for $k\in \lbrack K]$ . Similarly, we define $\mathbf {Y}$ and $\mathbf {U}$ . Then we have the matrix representation:

(2.5) $$ \begin{align} \mathbf{Matrix}\text{ }\mathbf{representation}:\quad \mathbf{Y} =\sum_{k=1}^{K}(\dot{\mathbf{B}}_{k,T}\mathbf{X}_{k}+\mathbf{X}_{k}\dot{ \mathbf{B}}_{k,N})+\Theta ^{0}+\mathbf{U}, \end{align} $$

where $\dot {\mathbf {B}}_{k,T}\equiv \mathrm {diag}(\dot {\beta }_{k,\cdot 1},\dot { \beta }_{k,\cdot 2},\ldots ,\dot {\beta }_{k,\cdot T})$ and $\dot {\mathbf {B}} _{k,N}\equiv \mathrm {diag}(\dot {\beta }_{k,1\cdot },\dot {\beta }_{k,2\cdot },\ldots ,\dot {\beta }_{k,N-1\cdot },-\sum _{i=1}^{N-1}\dot {\beta }_{k,i\cdot }) $ are $T\times T$ and $N\times N$ diagonal matrices, respectively.

We now outline the main steps for estimating model (2.1):

1. 1. Solve the nuclear-norm-regularized problem defined in (2.7) below. This step yields a consistent estimator of the matrix $\Theta ^{0}$ in terms of Frobenius norm.

2. 2. Determine the number of factors using (2.8) or other methods discussed below and obtain the initial estimators of the factors and factor loadings as in (2.9). This step provides consistent estimators of the number of factors, factors, and loadings.

3. 3. Apply Algorithm 1, introduced in Section 2.3, to obtain the final estimators of the slope coefficients and latent factors and loadings. This step refines the estimators through iterations and produces the final estimators used for statistical inference.

2.2 Initial Estimation: Nuclear Norm Regularized Estimation

In this section, we produce initial estimators of the parameters in the model. Define the objective function as

(2.6) $$ \begin{align} \mathcal{L}_{NT}(B,\Theta )\equiv \frac{1}{2NT}||Y-X^{(1,2)}B-\theta ||^{2}+\gamma _{NT}\frac{\left\vert \left\vert \Theta \right\vert \right\vert _{\text{*}}}{\sqrt{NT}}, \end{align} $$

where $B\in \mathbb {R}^{(N-1+T)K}$ , $\theta \equiv $ vec $(\Theta )$ , $\Theta \in \mathbb {R}^{T\times N}$ , and $\gamma _{NT}$ is a regularization parameter such that $\gamma _{NT}\rightarrow 0$ as $\left (N,T\right ) \rightarrow \infty $ . Then the initial estimator is given by

(2.7) $$ \begin{align} (\hat{B}_{o},\hat{\Theta}_{o})\equiv \arg \min_{(B,\Theta )}\mathcal{L} _{NT}(B,\Theta ). \end{align} $$

At any fixed $\Theta $ , minimizing $\mathcal {L}_{NT}(B,\Theta )$ is a least-squares problem. Formally, we can write

$$ \begin{align*} \mathcal{\tilde{L}}_{NT}(\Theta )=\min_{B}\frac{1}{2NT}||Y-X^{(1,2)}B-\theta ||^{2}=\frac{1}{2NT}(Y-\theta )^{\prime }\mathbb{M}_{x}(Y-\theta ), \end{align*} $$

where $\mathbb {M}_{x}\equiv \mathbb {M}_{X^{(1,2)}}$ is the $NT\times NT$ projection matrix that projects out the columns of $X^{(1,2)}.$ Hence, minimizing $\mathcal {L}_{NT}(B,\Theta )$ can be achieved by minimizing $ \mathcal {\tilde {L}}_{NT}(\Theta )+\gamma _{NT}\frac {\left \vert \left \vert \Theta \right \vert \right \vert _{\text {*}}}{\sqrt {NT}}.$ That is,

$$\begin{align*}\hat{\Theta}_{o} &=\arg \min_{\Theta }\mathcal{\tilde{L}}_{NT}(\Theta )+\gamma _{NT}\frac{\left\vert \left\vert \Theta \right\vert \right\vert _{\text{*}}}{ \sqrt{NT}} \\ &=\arg \min_{\Theta }\frac{1}{2NT}(Y-\theta )^{\prime }\mathbb{M} _{x}(Y-\theta )+\gamma _{NT}\frac{\left\vert \left\vert \Theta \right\vert \right\vert _{\text{*}}}{\sqrt{NT}}. \end{align*}$$

Given the initial estimators $\hat {B}_{o}$ and $\hat {\Theta } _{o} $ , one can use various methods to determine the number of factors $ R^{0}.$ For example, we can estimate $R^{0}$ via the singular value thresholding as follows:

(2.8) $$ \begin{align} \hat{R}=\sum_{i=1}^{N\wedge T}\mathbf{1}\left\{ \mu _{i}(\hat{\Theta} _{o})\geq \left(\gamma _{NT}\sqrt{NT}\left\Vert \hat{\Theta}_{o}\right\Vert _{\text{sp}}\right) ^{1/2}\right\} , \end{align} $$

where $\mu _{i}(\hat {\Theta }_{o})$ denotes the ith largest singular value of $\hat {\Theta }_{o}.$ Alternatively, we can estimate $R^{0}$ by some other methods, including the information criteria and panel criteria methods of Bai and Ng (Reference Bai and Ng2002), the edge distribution estimator of Onatski (Reference Onatski2010), the eigenvalue ratio and growth ratio method of Ahn and Horenstein (Reference Ahn and Horenstein2013), and the cross-validation method of Jin, Miao, and Su (Reference Jin, Miao and Su2021).

Given $\hat {R}$ , we can obtain preliminary estimators of the factors and factor loadings. Consider the singular value decomposition (SVD) of $\hat {\Theta }_{o}: \hat {\Theta }_{o}=\hat {L}\hat {D}\hat {V}^{\prime },$ where $\hat {D}=$ diag $(\mu _{1}(\hat {\Theta }_{o}),\ldots ,\mu _{N\wedge T}(\hat { \Theta }_{o})),$ and $\hat {L}$ and $\hat {V}$ stack the left and right singular vectors of $\hat {\Theta }_{o}$ , respectively. We estimate $F^{0}$ and $\Lambda ^{0}$ , respectively, by

(2.9) $$ \begin{align} \hat{F}_{o}=\sqrt{T}\hat{L}_{\ast ,[\hat{R}]}\text{ and }\hat{\Lambda}_{o}= \hat{V}_{\ast ,[\hat{R}]}\hat{D}_{[\hat{R}],[\hat{R}]}/\sqrt{T}, \end{align} $$

where $\hat {L}_{\ast ,[\hat {R}]}$ denotes a submatrix of $\hat {L}$ that contains its first $\hat {R}$ columns and $\hat {D}_{[\hat {R}],[\hat {R} ]}\equiv $ diag $(\mu _{1}(\hat {\Theta }_{o}),\dots , \mu _{\hat {R}}(\hat {\Theta } _{o})).$

2.3 Post Nuclear Norm Regularized Estimation

In this section we propose an iterative procedure to obtain updated estimators of the regression coefficients, the factors, and the factor loadings. The algorithm goes as follows:

Note that the iterations in Algorithm 1 alternate between updating $\hat {B} _{s}$ and $(\hat {F}^{(s)},\hat {\Lambda }^{(s)})$ for each iteration index s. However, the algorithm does not require numerical convergence. Instead, it can be terminated after any finite number of iterations $s^{\ast }\geq 2$ . The advantages of this approach are twofold:

1. 1. Theoretical Perspective: We do not require the solution to correspond to the global minimum of an optimization problem. As long as the initial estimators of the factors and factor loadings are consistent, we can show that a finite number of iterations yields the desired asymptotic properties. Moreover, our theoretical analysis indicates that stopping at $s^{\ast }=2$ is sufficient, as additional iterations do not improve the asymptotic efficiency of the estimators in theory.

2. 2. Practical Perspective: In implementation, it is unnecessary to verify numerical convergence, nor is there a concern that the algorithm may become trapped in a local minimum.

3.1 Asymptotic Properties of the Initial Estimators

To study the asymptotic properties of the initial estimators, we make the following assumptions.

Assumption 1 imposes some regularity conditions on $\mathbb {M} _{x}U.$ It is close to the condition imposed on the tuning parameter in Lemma 2 of Moon and Weidner (Reference Moon and Weidner2026), which studies the model with homogeneous slope coefficients. As one can argue, it suffices to choose $\gamma _{NT}=C_{\gamma }(N^{-1/2}+T^{-1/2})$ for some sufficiently large positive constant $C_{\gamma }$ to ensure Assumption 1 to be satisfied under some primitive conditions on $\{x_{it}, u_{it}\}$ . Such a rate ensures that the NNR estimator of $\Theta $ exhibits the optimal rate of convergence in terms of Frobenius norm (see Theorem 3.1 below). In our simulations and application, we follow the lead of Miao et al. (Reference Miao, Phillips and Su2023) and set $\gamma _{NT}=0.5\max (N^{-1/2},T^{-1/2})$ for simplicity, and leave the optimal choice of $\gamma _{NT}$ for future research.

Assumption 2 imposes the celebrated restricted strong convexity (RSC) condition that was first introduced by Negahban and Wainwright (Reference Negahban and Wainwright2012) and Negahban et al. (Reference Negahban, Ravikumar, Wainwright and Yu2012). We conjecture that this condition can be verified under more primitive assumptions by adapting the approach of Chernozhukov et al. (Reference Chernozhukov, Hansen, Liao and Zhu2026), accounting for the additional complexity introduced by the high-dimensional structure of our slope coefficients (see also Hong, Su, and Jiang, Reference Hong, Su and Jiang2023; Moon and Weidner, Reference Moon and Weidner2026).

The following theorem studies the consistency of the nuclear-norm regularized estimator $\hat {\Theta }_{o}$ in terms of the Frobenius norm.

Theorem 3.1 indicates that we can estimate the large-dimensional matrix $\Theta ^{0}$ consistently in terms of the Frobenius norm. Note that $ \Theta ^{0}$ is a high-dimensional matrix with $NT$ entries and $\Vert \Theta ^{0}\Vert =O_{P}(\sqrt {NT})$ in general. This explains why we need to normalize $\Vert \hat {\Theta }_{o}-\Theta ^{0}\Vert $ by $\sqrt {NT}$ in Theorem 3.1.

Let $\max _{i}=\max _{1\leq i\leq N}$ and $\max _{t}=\max _{1\leq t\leq T}$ . Define $\max _{i,j}$ , $\max _{i,t}$ , $\max _{t,s}$ , $\min _{i}$ , and $\min _{t}$ analogously. Next, we impose an assumption on the common factors and factor loadings.

Assumption 3(i) and (ii) posits a strong factor structure, a standard condition widely adopted in the literature (see, e.g., Bai and Ng, Reference Bai and Ng2002; Bai, Reference Bai2003, Reference Bai2009; Su et al., Reference Su, Shi and Phillips2016; Moon and Weidner, Reference Moon and Weidner2017; Miao et al., Reference Miao, Phillips and Su2023, among others). Armstrong et al. (Reference Armstrong, Weidner and Zeleneev2026) propose a robust method for the estimation and inference of the slope coefficients in homogeneous panels involving weak factors. Assumption 3(iii) further imposes a uniform bound on the latent factors and factor loadings, which serves to facilitate the theoretical analysis in the subsequent sections. It holds when both $\lambda _{i}^{0}$ and $f_{t}^{0}$ are sub-exponentially distributed.

The following theorem reports the consistency of $\hat {R}, \hat {F}_{o},$ and $\hat {\Lambda }_{o}.$

Theorem 3.2(i) indicates that $\hat {R}$ is a consistent estimator of $R^{0}.$ Theorem 3.2(ii) indicates that the columns of $\hat {F}_{o}$ and $F^{0}$ span the same space asymptotically, so do the columns of $\hat {\Lambda }_{o}$ and $\Lambda ^{0}.$

3.2 Asymptotic Properties of the Infeasible Estimator

Despite the consistency of the initial estimators $\hat {F}^{(0)}$ and $\hat { \Lambda }^{(0)},$ it is challenging to study the asymptotic properties of $ \hat {\beta }_{i\cdot }$ and $\hat {\beta }_{\cdot t}$ directly. Below we first study the asymptotic properties of their infeasible version as given in (3.1) and then study those of the feasible one. Our results show that the difference between the feasible and infeasible versions is asymptotically negligible. As an intermediate result, we first study the infeasible estimator in (3.1) in this section. Now we formally introduce the infeasible estimator. Given $F^{0}$ and $\Lambda ^{0}$ are observed, we can estimate $\dot {B}$ by

(3.1) $$ \begin{align} \check{B}\equiv \left[ X^{(1,2)\prime }(\mathbb{M}_{\Lambda ^{0}}\otimes \mathbb{M}_{F^{0}})X^{(1,2)}\right] ^{-1}X^{(1,2)\prime }(\mathbb{M} _{\Lambda ^{0}}\otimes \mathbb{M}_{F^{0}})Y. \end{align} $$

Similar to the definition $\dot {B}=[\dot {\beta }_{1\cdot }^{\prime },\ldots ,\dot { \beta }_{N-1\cdot }^{\prime },\dot {\beta }_{\cdot 1}^{\prime }, \ldots ,\dot { \beta }_{\cdot T}^{\prime }]^{\prime }$ , we set $\check {B}\equiv \lbrack \check {\beta }_{1\cdot }^{\prime }, \ldots ,\check {\beta }_{N-1\cdot }^{\prime }, \check {\beta }_{\cdot 1}^{\prime },\ldots ,\check {\beta }_{\cdot T}^{\prime }]^{\prime }$ . The infeasible estimator for $\dot {\beta }_{N\cdot }$ is given by $\check {\beta }_{N\cdot }=-\sum _{i=1}^{N-1}\check {\beta }_{i\cdot }$ . In Section S.1 of the Supplementary Material, we show that $X^{(1,2)\prime }(\mathbb {M}_{\Lambda ^{0}}\otimes \mathbb {M}_{F^{0}})X^{(1,2)}$ is invertible asymptotically. Therefore, $\check {B}$ is well defined in large samples.

Remark 3.3. The infeasible estimator in equation (3.1) can be applied to estimate a heterogeneous panel data model that includes both observed time-invariant regressors $\lambda _{i}^{0}$ and factors $f_{t}^{0}$ :

$$ \begin{align*} y_{it}=(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0})^{\prime }x_{it}+\lambda _{i}^{0\prime }g_{t}+\gamma _{i}^{\prime }f_{t}^{0}+u_{it}, \end{align*} $$

where $g_{t}$ and $\gamma _{i}$ are unknown coefficients.Footnote ³ In particular, the heterogeneous panel data model with two-way fixed effects (TWFE) model, as in Lu and Su (Reference Lu and Su2023), arises as a special case. When $\lambda _{i}^{0}=f_{t}^{0}=1,$ the IFE structure simplifies to $\lambda _{i}^{0\prime }g_{t}+\gamma _{i}^{\prime }f_{t}^{0}=g_{t}+\gamma _{i}$ , which represents the standard TWFE specification. An important distinction lies in the assumptions regarding the regressors $x_{it}.$ Lu and Su (Reference Lu and Su2023) impose the condition that $x_{it}$ ’s behave as idiosyncratic errors after projecting out the individual and time effects. Our assumptions on $x_{it}$ are more general and less restrictive. Consequently, the theoretical results developed in this section also apply directly to inference on the slope coefficients $ \beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}$ within the specification described above.

To study the asymptotic properties of the infeasible estimator, we first introduce some notation and impose two regularity assumptions. Let $\dot { \mathbf {x}}_{i\cdot }\equiv \mathbb {M}_{F^{0}}\mathbf {x}_{i\cdot }$ , and let $\dot {x}_{it}$ denote the transpose of the t’th row of $\dot {\mathbf {x}} _{i\cdot }$ . Similarly, define $\ddot {\mathbf {x}}_{\cdot t}\equiv \mathbb {M} _{\Lambda ^{0}}\mathbf {x}_{\cdot t}$ and $\ddot {x}_{it}$ . Define

(3.2) $$ \begin{align} \dot{\mathbf{s}}_{i\cdot } &\equiv \frac{1}{T}\sum_{t=1}^{T}\dot{x}_{it} \dot{x}_{it}^{\prime }=\frac{1}{T}\mathbf{x}_{i\cdot }^{\prime }\mathbb{M} _{F^{0}}\mathbf{x}_{i\cdot },\text{ }\ddot{\mathbf{s}}_{\cdot t}\equiv \frac{ 1}{N}\sum_{i=1}^{N}\ddot{x}_{it}\ddot{x}_{it}^{\prime }=\frac{1}{N}\mathbf{x} _{\cdot t}^{\prime }\mathbb{M}_{\Lambda ^{0}}\mathbf{x}_{\cdot t}, \notag \\ \dot{q}_{i} &\equiv \mathbf{x}_{i\cdot }^{\prime }\mathbb{M}_{F^{0}}\mathbf{ u}_{i\cdot },\text{ and }\ddot{q}_{t}\equiv \mathbf{x}_{\cdot t}^{\prime } \mathbb{M}_{\Lambda ^{0}}\mathbf{u}_{\cdot t}. \end{align} $$

Let $\mathbf {Q}_{12}$ be an $(N-1)\times T$ block-partitioned matrix whose $ (i,t)$ -th block (of size $K\times K$ ) is given by $\dot {x}_{it}\ddot {x} _{it}^{\prime }-\frac {1}{N-1}\sum _{j=1}^{N-1}\dot {x}_{jt}\ddot {x} _{jt}^{\prime },\text { for }i\in \lbrack N-1], t\in \lbrack T].$ Let $ \delta _{NT}=\min (\sqrt {N},\sqrt {T})$ .

Assumption 4 imposes regularity conditions on $\dot {\mathbf {s}} _{i\cdot }$ , $\ddot {\mathbf {s}}_{\cdot t}$ , $x_{it}$ , $f_{t}^{0}$ , and $ \lambda _{i}^{0}$ . Assumption 4(i) introduces specific regularity requirements on $\dot {\mathbf {s}}_{i\cdot }$ and $\ddot {\mathbf {s}}_{\cdot t} $ , which are constructed from the de-factored and de-factor-loading versions of $x_{it}$ , denoted by $\dot {x}_{it}$ and $\ddot {x}_{it}$ , respectively. These conditions can be verified under suitable primitive assumptions, although they do exclude certain cases. Assumption 4 (ii) holds under mild moment and weak dependence conditions on the sequence $ \{x_{it}\}$ . Assumption 4(iii) is central for ensuring the asymptotic non-singularity of the Gram matrix $X^{(1,2)\prime }(\mathbb {M} _{\Lambda }\otimes \mathbb {M}_{F})X^{(1,2)}$ . For discussions on low-rank covariates, see, for example, Moon and Weidner (Reference Moon and Weidner2017) and Moon and Weidner (Reference Moon and Weidner2026).

Remark 3.4. Assumption 4(iii) eliminates the possibility of low-rank regressors in $x_{it}$ as in Bai (Reference Bai2009) and Armstrong et al. (Reference Armstrong, Weidner and Zeleneev2026). In particular, rank-one regressors, including $x_{it}$ that vary only across i or only across t will certainly lead to an identification issue. We present a simple example here. Consider the model in Remark 3.3. Suppose a scalar $x_{it}$ has rank one. Then it can be expressed as $x_{it}=a_{i}e_{t}$ , where $a_{i}$ and $e_{t}$ are both scalars that take value 0 with zero probability. For this case, the model specification suffers from a fundamental identification problem. For any vectors $\delta _{1}$ and $\delta _{2}$ of conformable dimensions, we have

$$ \begin{align*} \dot{\beta}_{\cdot t}x_{it}+\gamma _{i}^{\prime }f_{t}^{0}=\dot{\beta} _{\cdot t}a_{i}e_{t}+\gamma _{i}^{\prime }f_{t}^{0}=\left(\dot{\beta} _{\cdot t}+\frac{\delta _{1}^{\prime }f_{t}^{0}}{e_{t}}\right) a_{i}e_{t}+(\gamma _{i}-a_{i}\delta _{1})^{\prime }f_{t}^{0}=\dot{\beta} _{\cdot t}^{\ast }x_{it}+\gamma _{i}^{\ast \prime }f_{t}^{0}, \end{align*} $$

and similarly,

$$ \begin{align*} \dot{\beta}_{i\cdot }x_{it}+\lambda _{i}^{0\prime }g_{t}=\dot{\beta}_{i\cdot }a_{i}e_{t}+\lambda _{i}^{0\prime }g_{t}=\left(\dot{\beta}_{i\cdot }+\frac{ \delta _{2}^{\prime }\lambda _{i}^{0}}{a_{i}}\right) a_{i}e_{t}+\lambda _{i}^{0\prime }(g_{t}-e_{t}\delta _{2})=\dot{\beta}_{i\cdot }^{\ast }x_{it}+\lambda _{i}^{0\prime }g_{t}^{\ast }, \end{align*} $$

where $\dot {\beta }_{\cdot t}^{\ast }\equiv \dot {\beta }_{\cdot t}+\frac { \delta _{1}^{\prime }f_{t}^{0}}{e_{t}}, \gamma _{i}^{\ast }\equiv \gamma _{i}-a_{i}\delta _{1}, \dot {\beta }_{i\cdot }^{\ast }\equiv \dot {\beta } _{i\cdot }+\frac {\delta _{2}^{\prime }\lambda _{i}^{0}}{a_{i}},\quad $ and $ g_{t}^{\ast }\equiv g_{t}-e_{t}\delta _{2}.$ This algebraic equivalence demonstrates that the parameters in the model described in Remark 3.3 are not identified in the presence of rank-one regressors. This issue directly results in singular Gram matrix $X^{(1,2)\prime }(\mathbb {M} _{\Lambda ^{0}}\otimes \mathbb {M}_{F^{0}})X^{(1,2)}$ . Further details are provided in Section S.6 of the Supplementary Material.

Assumption 5 imposes some high-level conditions on terms related to $u_{it}$ ’s. Assumption 5(i) is standard in the literature on factor models (see, e.g., Lu and Su, Reference Lu and Su2016; Moon and Weidner, Reference Moon and Weidner2017, Reference Moon and Weidner2026; Su and Wang, Reference Su and Wang2017; Chernozhukov et al., Reference Chernozhukov, Hansen, Liao and Zhu2026). Moon and Weidner (Reference Moon and Weidner2017) provide several examples to verify it. The rate conditions in Assumption 5(ii)–(vii) can be verified under some primitive conditions on $\left \{ x_{it},\text { } f_{t}^{0},\text { }u_{it}\right \} $ by calling up certain exponential inequalities. Assumption 5(ii) imposes pointwise asymptotic rate conditions on the relevant partial sums. In particular, $ \sum _{t=1}^{T}x_{it}u_{it}=O_{P}(\sqrt {T})$ is interpreted for each fixed i, and $\sum _{i=1}^{N}x_{it}u_{it}=O_{P}(\sqrt {N})$ is interpreted for each fixed t; neither condition is intended to hold uniformly over its free index. The aggregate rate conditions required in our analysis are imposed separately in Assumption 5(iii) and (iv), whereas Assumption 5(v)–(vii) provides the required uniform bounds through maxima over the relevant indices.

The following theorem establishes the asymptotic linear expansions of $ \check {\beta }_{i\cdot }$ and $\check {\beta }_{\cdot t}.$

Theorem 3.3 establishes the asymptotic leading term for the infeasible estimators. Under Assumptions 3–5, our theoretical analysis focuses on bounding the higher-order remainders in the expansions of $\check {\beta }_{i\cdot }-\dot {\beta }_{i\cdot }$ and $\check { \beta }_{\cdot t}-\dot {\beta }_{\cdot t}$ . For the relative growth of N and $ T$ , it suffices that $N^{-1/2}\log ^{2}T+T^{-1/2}\log ^{2}N\rightarrow 0$ , which contrasts with the stronger rate requirement in Theorem 3.4.

3.3 Asymptotic Properties of the Feasible Estimators

The next theorem reports the linear expansions of $\hat {\beta }_{i\cdot }- \dot {\beta }_{i\cdot }$ and $\hat {\beta }_{\cdot t}-\dot {\beta }_{\cdot t}.$

Theorem 3.4 presents the uniform asymptotic expansions of $\hat { \beta }_{i\cdot }-\dot {\beta }_{i\cdot }$ and $\hat {\beta }_{\cdot t}-\dot {\beta }_{\cdot t}$ under Assumptions 1–5, along with an additional rate condition on $(N,T)$ . The required rate condition is stronger than that in Theorem 3.3, reflecting the need to estimate the latent factors rather than treat them as observed. A standard rate condition in the literature on factor models is $(\sqrt {N}/T+\sqrt {T} /N)\rightarrow 0$ . Given the greater complexity of our model and estimation procedure, we impose the slightly stronger condition $(\sqrt {N}/T+\sqrt {T} /N)\log ^{2}(NT)\rightarrow 0$ to satisfy the technical demands of our asymptotic analysis.

We impose the following assumption to enable inference on the slope coefficients, factors, and factor loadings.

Assumption 6(i) imposes moment and conditional mean restrictions on $\beta _{i\cdot }^{0}$ and $u_{it}$ , while Assumption 6(ii) imposes uniform rate conditions on $\dot {\mathbf {s}}_{i\cdot }$ and $\ddot {\mathbf {s}}_{\cdot t}$ defined in (3.2). Assumption 6(iii) imposes central limit theorems for the aggregate slope heterogeneity and the score terms entering the estimators of the unit- and time-specific slope coefficients. Assumption 6(iv)–(vi) imposes probability limits and central limit theorems for terms entering the estimation of the factors and factor loadings. All convergence statements in Assumption 6(iii)–(vi) that contain a free index are pointwise. These conditions do not assert joint convergence over multiple units or time periods. Accordingly, Theorem 3.5 establishes pointwise asymptotic normality for each fixed unit-specific coefficient and each fixed time-specific coefficient, rather than joint or simultaneous inference over $ i$ or t.

The following theorem reports the asymptotic distributions of $\hat {\beta } _{i\cdot }$ and $\hat {\beta }_{\cdot t}.$

Theorem 3.5 provides pointwise inference only for each fixed i or t; it does not establish joint or simultaneous inference across multiple units or time periods. Such inference, including inference after data-dependent selection, requires additional conditions. For the cell-level slope $\beta _{it}^{0}=\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}$ , the natural rate of $\hat {\beta }_{it}-\beta _{it}^{0}$ is $\min \{ \sqrt {N},\sqrt {T}\}$ . If $N/T\rightarrow \infty $ (resp. $ N/T\rightarrow 0$ ), the unit-specific (resp. time-specific) error is the only leading term, and its marginal limit in Theorem 3.5 determines the limiting distribution. If $N/T\rightarrow \kappa \in (0,\infty )$ , both terms contribute. Because Theorem 3.5 does not guarantee their joint asymptotic normality, in this case, the expansions in Theorem 3.4 must be supplemented by the additional condition that the two leading terms are asymptotically uncorrelated.

Uniform confidence bands for the collection $\{\beta _{i\cdot }^{0}:1\leq i\leq N\}$ may also be of practical interest. However, the uniform asymptotic expansions in Theorem 3.4 and the marginal convergence results in Theorem 3.5 are not sufficient to justify such inference. Uniform inference requires additional conditions and theoretical results, including a uniform Gaussian approximation and the validity of an appropriate bootstrap procedure. Moreover, this problem is more challenging than that in Lu and Su (Reference Lu and Su2023), as the latent factors and factor loadings must be estimated, and their estimation errors need to be controlled uniformly alongside those of the heterogeneous slope coefficients. We therefore leave the development of uniform inference to future research.

Given the results in Theorem 3.5, it is standard to construct pointwise confidence intervals (CIs) for the elements of $\beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}$ . To this end, we need to estimate the parameters that appear in the variance matrices $V_{\beta ,i\cdot }$ and $ V_{\beta ,\cdot t}$ :

• • The quantities $s_{i\cdot }^{0}$ and $s_{\cdot t}^{0}$ can be estimated by $\hat {\mathbf {s}}_{i\cdot }\equiv x_{i\cdot }^{\prime }M_{\hat {F }}x_{i\cdot }/T$ and $\hat {\mathbf {s}}_{\cdot t}\equiv x_{\cdot t}^{\prime }M_{\hat {\Lambda }}x_{\cdot t}/N,$ respectively.

• • Under the additional assumption that the $\beta _{i\cdot }^{0}$ ’s are independent across i, $\Sigma _{1\beta }$ can be estimated by the sample variance of $\hat {\beta }_{i\cdot }$ .

• • Under the additional assumption that conditional on $(\Lambda ^{0},x_{\cdot t}), u_{it}$ ’s are uncorrelated across i, $\Phi _{\cdot t}$ equals the probability limit of $\frac {1}{N}\sum _{i=1}^{N}\ddot {x}_{it}\ddot { x}_{it}^{\prime }u_{it}^{2}.$ An estimator of $\Phi _{\cdot t}$ is given by $ \hat {\Phi }_{\cdot t}\equiv \frac {1}{N}\sum _{i=1}^{N}\widehat {\ddot {x}}_{it} \widehat {\ddot {x}}_{it}^{\prime }\hat {u}_{it}^{2},$ where $\widehat {\ddot {x}} _{it}\equiv x_{it}-x_{\cdot t}^{\prime }\hat {\Lambda }(\hat {\Lambda }^{\prime } \hat {\Lambda })^{-1}\hat {\lambda }_{i}$ .

• • Under the additional assumption that conditional on $(F^{0},x_{i\cdot }), u_{it}$ ’s are uncorrelated across t, $\Phi _{i\cdot }$ equals the probability limit of $\frac {1}{T}\sum _{t=1}^{T}\dot {x}_{it}\dot {x} _{it}^{\prime }u_{it}^{2}.$ An estimator of $\Phi _{i\cdot }$ is given by $ \hat {\Phi }_{i\cdot }\equiv \frac {1}{T}\sum _{t=1}^{T}\widehat {\dot {x}}_{it} \widehat {\dot {x}}_{it}^{\prime }\hat {u}_{it}^{2},$ where $\widehat {\dot {x}} _{it}\equiv x_{it}-x_{i\cdot }^{\prime }\hat {F}(\hat {F}^{\prime }\hat {F} )^{-1}\hat {f}_{t}$ .

  If serial correlation is a concern, $\Phi _{i\cdot }$ represents the long-run variance of $\frac {1}{\sqrt {T}}\sum _{t=1}^{T}\dot {x}_{it}u_{it}$ . In such cases, a heteroskedasticity-and-autocorrelation-consistent (HAC) estimator should be employed to estimate $\Phi _{i\cdot }$ .

Remark 3.6. We now compare the results in Theorem 3.5 with those in Chernozhukov et al. (Reference Chernozhukov, Hansen, Liao and Zhu2026). Since we assume additive structure, $\beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}$ can be separately identified under the location normalization $E(\beta _{i\cdot }^{0})=0$ . This makes it possible to obtain the linear expansions of $\hat { \beta }_{i\cdot }-\beta _{i\cdot }^{0}$ and $\hat {\beta }_{\cdot t}-\beta _{\cdot t}^{0}$ in Theorem 3.4 and to derive the asymptotic distributions of $\sqrt {T}(\hat {\beta }_{i\cdot }-\beta _{i\cdot }^{0})$ and $ \sqrt {N}(\hat {\beta }_{\cdot t}-\beta _{\cdot t}^{0})$ in Theorem 3.5. Given these results, we can also, under additional assumptions, study the average effects by averaging $\beta _{i\cdot }^{0}$ over a group of individuals or $\beta _{\cdot t}^{0}$ over a group of time periods. In contrast, Chernozhukov et al. (Reference Chernozhukov, Hansen, Liao and Zhu2026) assume a low-rank structure for $\beta _{it}^{0}$ and focus on the estimation of certain group average effects at time t. For example, one object of interest is $\bar {\beta }_{\mathcal {G},t}^{0}=\frac {1}{\left \vert \mathcal {G} \right \vert }\sum _{i\in \mathcal {G}}\beta _{it}^{0},$ where $\mathcal {G}$ is a subset of $\left [ N\right ] $ and $\left \vert \mathcal {G}\right \vert $ denotes the cardinality of the set $\mathcal {G}$ . Under the restrictive condition that $\left \vert \mathcal {G}\right \vert =o(N)$ along with some other restrictive conditions, they derive the asymptotic normality of

$$ \begin{align*} \left\vert \mathcal{G}\right\vert ^{1/2}\left(\frac{1}{\left\vert \mathcal{ \ G}\right\vert }\sum_{i\in \mathcal{G}}\hat{\beta}_{it}-\bar{\beta}_{ \mathcal{\ G},t}^{0}\right). \end{align*} $$

When $\beta _{it}^{0}=\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}, \bar {\beta }_{\mathcal {G},t}^{0}=\frac {1}{\left \vert \mathcal {G}\right \vert }\sum _{i\in \mathcal {G}}\beta _{it}^{0}=\frac {1}{\left \vert \mathcal {G}\right \vert } \sum _{i\in \mathcal {G}}\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}.$ In order to identify the time effect $\beta _{\cdot t}^{0}$ , one needs to specify $ \mathcal {G}$ to be $\left [ N\right ] ,$ which is unfortunately ruled out by the condition that $\left \vert \mathcal {G}\right \vert =o(N)$ in Chernozhukov et al. (Reference Chernozhukov, Hansen, Liao and Zhu2026). So the results in our article and Chernozhukov et al. (Reference Chernozhukov, Hansen, Liao and Zhu2026) are quite different from each other and they are derived under different sets of assumptions. See the next remark.

For the group average $\bar {\beta }_{\mathcal {G},t}^{0}$ , define $\widehat {\bar {\beta }}_{\mathcal {G},t}\equiv |\mathcal {G} |^{-1}\sum _{i\in \mathcal {G}}\hat {\beta }_{it}$ . Then

$$ \begin{align*} \widehat{\bar{\beta}}_{\mathcal{G},t}-\bar{\beta}_{\mathcal{G},t}^{0}=\frac{1 }{|\mathcal{G}|}\sum_{i\in \mathcal{G}}\left(\hat{\beta}_{i\cdot }-\beta _{i\cdot }^{0}\right) +\left(\hat{\beta}_{\cdot t}-\beta _{\cdot t}^{0}\right). \end{align*} $$

Conducting inference on this quantity thus requires an aggregate central limit theorem for the first term, as well as a characterization of its joint limiting behavior with the second. Under the conventional $\sqrt {T|\mathcal {G}|}$ rate for the first term, the relative magnitudes of $T|\mathcal {G}|$ and N determine which component dominates asymptotically. Furthermore, if the target of inference is a population—rather than realized—group mean, an additional central limit theorem is needed for the centered mean of the true effects $\beta _{i\cdot }^{0}$ ’s. These requirements are not implied by Theorems 3.4 and 3.5; rather, they play a role analogous to that of Assumption 8(vi).

To conduct inference on the factors and factor loadings, we first introduce some notation. Define $V=\mathrm {diag}(v_{1},\ldots ,v_{R^{0}})$ and $ Q=V^{1/2}\Upsilon ^{\prime }\Sigma _{\Lambda }^{-1/2}$ , where $ v_{1}>v_{2}>\cdots >v_{R^{0}}>0$ are the eigenvalues of $\Sigma _{\Lambda }^{1/2}\Sigma _{F}\Sigma _{\Lambda }^{1/2}$ , and $\Upsilon $ is the corresponding matrix of eigenvectors satisfying $\Upsilon ^{\prime }\Upsilon =\mathbb {I}_{R^{0}}$ . Let $\tilde {\mathbf {Y}}=\mathbf {Y}-\mathrm {mat}\!\big (X^{(1,2)}\hat {B}\big )$ , and let $V_{NT}$ denote the $R^{0}\times R^{0}$ diagonal matrix containing the largest $R^{0}$ eigenvalues of $(NT)^{-1} \tilde {\mathbf {Y}}\tilde {\mathbf {Y}}^{\prime }$ , arranged in decreasing order. Define $H=\big (N^{-1}\Lambda ^{0^{\prime }}\Lambda ^{0}\big )\big (T^{-1}F^{0^{\prime }}\hat {F}\big )V_{NT}^{-1}.$ The following theorem establishes the asymptotic distributions of $\hat {f}_{t}$ , $\hat {\lambda } _{i} $ , and the product $\hat {\lambda }_{i}^{\prime }\hat {f}_{t}$ .

Theorem 3.6 provides pointwise distribution results. It does not establish joint or simultaneous inference across multiple time periods or units. In particular, uniform inference on the factor trajectory $\{f_{t}^{0}\}_{t\leq T}$ under the adopted normalization (equivalently, on the rotated trajectory $\{H^{\prime }f_{t}^{0}\}_{t\leq T}$ ) requires additional conditions and a uniform asymptotic linear expansion for $\hat {f}_{t}-H^{\prime }f_{t}^{0}$ , analogous to the expansions in Theorem 3.4, together with a uniform Gaussian approximation and valid critical values for the corresponding maximum statistic. We leave the development of such uniform inference to future research.

The covariance matrices $\Gamma _{\cdot t}$ and $\Gamma _{\cdot i}$ , as in Assumption 6(v) and (vi), are different from that often appear in the factor literature. This point indicates the slow convergences of $\hat {\beta } _{\cdot t}$ and $\hat {\beta }_{i\cdot }$ enter the asymptotics of $\hat {f} _{t} $ and $\hat {\lambda }_{i},$ respectively, and both enter the asymptotic distribution of the estimator of the common component $\lambda _{i}^{0\prime }f_{t}^{0}.$

4.1 Basic Setup

One key application of our asymptotic results is to estimate various heterogeneous causal effects. To illustrate the main idea, we consider a simple setup where we have $N=N_{1}+N_{2}$ individuals where the first $ N_{1} $ individuals $\left (i\in \lbrack N_{1}]\right ) $ never receive any treatment, and the last $N_{2}$ individuals $(i\in \left [ N\right ] \backslash \left [ N_{1}\right ] \equiv \{N_{1}+1,\ldots ,N\})$ receive treatment at some potentially different time points. We focus on the case where once a unit adopts the treatment, it stays treated afterward. We denote by $T_{0,i}$ and $T_{1,i}$ the number of control and treated time periods for unit i such that $T_{0,i}+T_{1,i}=T.$ If individual i never receives treatment, then $T_{0,i}=T$ and $T_{1,i}=0$ . Following the literature of the potential outcome framework, we use $y_{it}^{1}$ and $ y_{it}^{0}$ to denote the potential outcome for individual i who receives and does not receive the treatment at time t, respectively. Let $y_{it}$ be the observed outcome: $y_{it}=y_{it}^{1}d_{it}+y_{it}^{0}\left (1-d_{it}\right ) ,$ where $d_{it}=\mathbf {1}\{$ individual i receives treatment at time $t\}$ is the treatment indicator. Let $\mathcal {D} ^{0}\equiv \{(i,t)\mid d_{it}=0,(i,t)\in \lbrack N]\times \lbrack T]\}$ and $ \mathcal {D}^{1}\equiv ([N]\times \lbrack T])\backslash \mathcal {D}^{0}$ . Thus, we observe $y_{it}^{0}$ ’s for $(i,t)\in \mathcal {D}^{0}$ and $ y_{it}^{1}$ ’s for $(i,t)\in \mathcal {D}^{1}$ , respectively.

The individual treatment effect measures the difference between the treated and control outcome:

$$ \begin{align*} \tau _{it}=y_{it}^{1}-y_{it}^{0}\text{. } \end{align*} $$

We do not impose any assumptions on $y_{it}^{1}$ and therefore $ \tau _{it}$ . In particular, $y_{it}^{1}$ does not need to follow a model like that for $y_{it}^{0}$ . $\tau _{it}$ can be heterogeneous across both i and t and it does not necessarily have a factor structure as assumed in Xiong and Pelger (Reference Xiong and Pelger2023). We are only interested in the treatment effect on the treated, that is, $\tau _{it}$ for $(i,t)\in \mathcal {D}^{1}$ . For them, we directly observe $y_{it}^{1},$ and therefore we only need to model $ y_{it}^{0}.$ We assume that $y_{it}^{0}$ is generated via the following data-generating process (DGP):

(4.1) $$ \begin{align} y_{it}^{0}=\left(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}\right) ^{\prime }x_{it}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it}\text{ for }i\in \left[ N \right] \text{ and }t\in \left[ T\right] , \end{align} $$

where $x_{it}$ , $\beta _{i\cdot }^{0},\beta _{\cdot t}^{0}, \lambda _{i}^{0}, f_{t}^{0},$ and $u_{it}$ are as defined in the last section. As we do not observe all $y_{it}^{0}$ ’s for $(i,t)\in \lbrack N]\times \lbrack T]$ , model (4.1) assists in imputing unobserved $y_{it}^{0}$ ’s.

Nevertheless, identification of the counterfactual untreated outcomes relies on model (4.1), which is assumed to hold for all $ (i,t)$ , including treated cells for which $y_{it}^{0}$ is counterfactual. In particular, conditional on $x_{it}$ , the time-specific slope $\beta _{\cdot t}^{0}$ and common factor $f_{t}^{0}$ are shared across all units. Thus, rather than dispensing with a parallel-trends-type restriction, our conditional factor model for untreated outcomes imposes a factor-augmented stability condition in the same spirit as the conditional parallel-trends assumption in the staggered difference-in-differences literature.

Let $\mathcal {G}_{t}\equiv \left \{ i\in \left [ N\right ] \backslash \left [ N_{1}\right ] ,\text { }d_{it}=1\right \} $ and $\bar {T}_{1}=\frac {1}{N_{2}} \sum _{i=N_{1}+1}^{N}T_{1,i}.$ Let $\left \vert \mathcal {G}_{t}\right \vert $ denote the cardinality of the set $\mathcal {G}_{t},$ that is, the number of individuals who receive treatment at time $t.$ In the following, we consider different levels of ATT:

1. 1. Average treatment effect on the treated for i: $\tau _{i\cdot }= \frac {1}{T_{1,i}}\sum _{t=T_{0,i}+1}^{T}\tau _{it}.$

2. 2. Average treatment effect on the treated at time $t:\tau _{\cdot t}= \frac {1}{\left \vert \mathcal {G}_{t}\right \vert }\sum _{i\in \mathcal {G} _{t}}\tau _{it}.$

3. 3. Aggregate average treatment effect on the treated: $\bar {\tau }=\frac {1 }{N_{2}\bar {T}_{1}}\sum _{i=N_{1}+1}^{N}\sum _{t=T_{0,i}+1}^{T}\tau _{it}.$

Since we do not make any assumption on $\tau _{it},$ it is impossible to eliminate the effect of $u_{it}$ on the estimator of $\tau _{it}$ for a single cross-sectional unit at a time period. This motivates the above consideration of ATT. Of course, one can also consider other forms of treatment effects, say, weighted treatment effects, or treatment effects for certain selected groups of individuals or time periods. For each of the above three types of treatment effects, we derive the asymptotic distributions based on which we can test for the null hypothesis of no causal effects.

Nevertheless, if one would like to assume that $u_{it}$ has the same distribution over either i or $t,$ or both, we argue below that it is still possible to make inference on $\tau _{it}$ by relying on some other method.

4.2 Estimation of Heterogeneous Treatment Effects

Due to the presence of the two-dimensional unobserved slope heterogeneity in (4.1), we assume that both $N_{1}$ and T pass to infinity. This enables us to apply the asymptotic results in Section 3 to derive the asymptotic properties of certain estimators based on the data $ \left \{ \left (y_{it},x_{it}\right ) \right \} _{i\in \left [ N_{1}\right ] ,t\in \left [ T\right ] }.$ Then we can propose an estimator for the counterfactual outcome $y_{it}^{0}$ for $i=N_{1}+1,\ldots ,N$ and $ t=T_{0,i}+1,\ldots ,T,$ based on model (4.2) below.

Let $\hat {\tau }_{it}=y_{it}-\hat {y}_{it}^{0}$ for $(i,t)\in \mathcal {D}^{1}.$ $\hat {\tau }_{it}$ is an estimator of $\tau _{it}$ that is asymptotically unbiased but generally inconsistent. Then we can estimate $\tau _{i\cdot },$ $\tau _{\cdot t}$ , and $\bar {\tau }$ as follows:

$$ \begin{align*} \hat{\tau}_{i\cdot }=\frac{1}{T_{1,i}}\sum_{t=T_{0,i}+1}^{T}\hat{\tau}_{it}, \text{ }\hat{\tau}_{\cdot t}=\frac{1}{\left\vert \mathcal{G}_{t}\right\vert } \sum_{i\in \mathcal{G}_{t}}\hat{\tau}_{it},\text{ and }\hat{\tau}=\frac{1}{ N_{2}\bar{T}_{1}}\sum_{i=N_{1}+1}^{N}\sum_{t=T_{0,i}+1}^{T}\hat{\tau}_{it}. \end{align*} $$

In the next section, we will study the asymptotic properties of $\hat {\tau } _{i\cdot }, \hat {\tau }_{\cdot t},$ and $\hat {\tau }$ under some regularity conditions, based on which one can make valid inferences on $\tau _{i\cdot }, \tau _{\cdot t}$ , and $\bar {\tau }.$

For notational simplicity, here we consider the staggered treatment timing. Our approach can actually handle any irregular treatment pattern as long as the untreated period is reasonably long for the treated units. In this case, we simply run a time-series regression using the data from the untreated period in Step 2 of Algorithm 2 and estimate the counterfactual outcome for the treated periods in Step 3. Neither step requires regularly spaced time-series data.

4.3 Asymptotic Distributions

In this section, we study the asymptotic distributions of $\hat {\tau } _{i\cdot }, \hat {\tau }_{\cdot t}$ , and $\hat {\tau }.$ We first state the potential outcome assumptions formally.

Assumption 7. (Conditions for Potential Outcomes)

(i) The untreated and treated potential outcomes, $y_{it}^{0}$ and $y_{it}^{1},$ are given by

$$ \begin{align*} \left\{ \begin{array}{l} y_{it}^{0}=\left(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}\right) ^{\prime }x_{it}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it} \\ y_{it}^{1}=y_{it}^{0}+\tau _{it}, \end{array} \right. \end{align*} $$

where $i\in \left [ N\right ] $ and $t\in \left [ T\right ] .$

(ii) $E\left (u_{it}|x_{it},\beta _{i\cdot }^{0},\beta _{\cdot t}^{0},F^{0},\Lambda ^{0},{{\{d_{is}\}_{s\leq t}}}\right ) =0.$

Assumption 7(i) is common in the synthetic control literature. For example, it is similar to Assumption 1 in Ferman and Pinto (Reference Ferman and Pinto2021) and Assumption 1 in Xu (Reference Xu2017) except that our model specification is more general with covariates and we allow for two-dimensional heterogeneity in the slope coefficients. Assumption 7(ii) strengthens the exogeneity Assumption 6(i) by requiring that $u_{it}$ be mean-independent of the treatment history $ \{d_{is}\}_{s\leq t}$ . This is similar to Assumption 3 in Ferman and Pinto (Reference Ferman and Pinto2021) and Assumption 2 in Xu (Reference Xu2017).

We now state some technical assumptions. With a little bit abuse of notation, we continue to use $\hat{\boldsymbol{B}}_{k,T}=$ diag( $\hat {\beta } _{k,\cdot 1}, \ldots ,\hat {\beta }_{k,\cdot T}$ ) and $\hat{\boldsymbol{B}} _{k,N_{1}}=$ diag( $\hat {\beta }_{k,1\cdot },\ldots ,\hat {\beta }_{k,N_{1}\cdot }$ ) to denote estimators of $\boldsymbol{\dot {B}}_{k,T}=$ diag $(\dot {\beta }_{k,\cdot 1}, \ldots ,\dot {\beta }_{k,\cdot T})$ and $\boldsymbol{\dot {B}}_{k,N_{1}}=$ diag( $ \dot {\beta }_{k,1\cdot },\ldots , \dot {\beta }_{k,N_{1}\cdot }$ ) based on the data $\left \{ \left (y_{it},x_{it}\right ) \right \} _{i\in \left [ N_{1}\right ] ,t\in \left [ T\right ] }.$ The corresponding matrices $\mathbf {Y}$ (or $ \mathbf {X}_{k})$ and B (or $\hat {B})$ are now of dimensions $ T\times N_{1}$ and $\left (N_{1}+T-1\right ) K\times 1,$ respectively. Let $ \boldsymbol{\tilde {Y}=} \mathbf {Y}-\sum _{k=1}^{K}(\hat{\boldsymbol{B}}_{k,T} \mathbf {X}_{k}+\mathbf {X}_{k}\hat{\boldsymbol{B}}_{k,N_{1}})=\mathbf {Y}-$ mat $ (X^{\left (1,2\right ) }\hat {B}).$ Let $V_{N_{1}T}$ be the $R^{0}\times R^{0}$ diagonal matrix of the first $R^{0}$ largest eigenvalues of $\left (N_{1}T\right ) ^{-1}\boldsymbol{\tilde {Y}\tilde {Y}}^{\prime }$ in decreasing order along its main diagonal line. Let $H =(N_{1}^{-1}\Lambda ^{0\prime }\Lambda ^{0})(T^{-1}F^{0\prime }\hat {F})V_{N_{1}T}^{-1},$ where $\Lambda ^{0}\equiv \Lambda _{N_{1}}^{0}\equiv \lbrack \lambda _{1}^{0},\ldots ,\lambda _{N_{1}}^{0}]^{\prime }.$ Let $z_{it}=\left (x_{it}^{\prime },f_{t}^{0\prime }H\right ) ^{\prime }$ and $\bar {z}_{it}=\left (x_{it}^{\prime },f_{t}^{0\prime }\right ) ^{\prime }.$ Similarly, $\mathbf {x}_{\cdot t}=[x_{1t},\ldots ,x_{N_{1}t}]^{\prime }$ is now an $N_{1}\times K$ matrix. Let $ \mathbb {S}_{1}\in \mathbb {R}^{K\times \left (K+R^{0}\right ) }$ and $\mathbb {S }_{2}\in \mathbb {R}^{R^{0}\times \left (K+R^{0}\right ) }$ denote selection matrices such that $\mathbb {S}_{1}A$ and $\mathbb {S}_{2}A$ select the first $ K$ rows and the last $R^{0}$ rows of $A,$ respectively. Define

$$ \begin{align*} \mu _{i\cdot ,x} &=E(x_{it}),\text{ }\mu _{f}=E\left(f_{t}^{0}\right) , \text{ }\mu _{\lambda }=E(\lambda _{i}^{0}),\text{ } \\ \sigma _{i\cdot ,u}^{2} &=\lim_{T_{1,i}\rightarrow \infty }\frac{1}{T_{1,i}} \sum_{t=T_{0,i}+1}^{T}\sum_{s=T_{0,i}+1}^{T}E\left(u_{it}u_{is}\right) , \\ \text{ }S_{i,T_{0,i}} &=\frac{1}{T_{0,i}}\sum_{t=1}^{T_{0,i}}z_{it}z_{it}^{ \prime },\text{ }\hat{\pi}_{t,x\lambda }=(\Lambda ^{0\prime }\Lambda ^{0})^{-1}\Lambda ^{0\prime }\mathbf{x}_{\cdot t}, \\ S_{i} &=\underset{T_{0,i}\rightarrow \infty }{\text{plim}}S_{i,T_{0,i}}, \\ \text{ }\pi _{t,x\lambda } &=\text{plim}_{N_{1}\rightarrow \infty }\hat{\pi} _{t,x\lambda },\text{ }\alpha _{i\cdot ,T}=\frac{T_{1,i}}{\sqrt{T_{0,i}\bar{T }_{1}}}(\mu _{i\cdot ,x}^{\prime }\mathbb{S}_{1}+\mu _{f}^{\prime }H\mathbb{S }_{2})S_{i,T_{0,i}}^{-1},\text{ } \\ \mu _{\cdot t,x\lambda } &=\frac{1}{N_{1}}\sum_{j=1}^{N_{1}}E(x_{jt}\lambda _{j}^{0\prime }),\text{ }\nu _{\cdot t,x\lambda }=\text{plim}_{\left\vert \mathcal{G}_{t}\right\vert \rightarrow \infty }\frac{1}{\left\vert \mathcal{G }_{t}\right\vert }\sum_{i\in \mathcal{G}_{t}}[x_{it}-\mathbf{x}_{\cdot t}^{\prime }\Lambda ^{0}\left(\Lambda ^{0\prime }\Lambda ^{0}\right) ^{-1}{ \lambda }_{i}^{0}], \\ \Omega _{i,zu} &=\lim_{T_{0,i}\rightarrow \infty }\frac{1}{T_{0,i}} \sum_{t=1}^{T_{0,i}}\sum_{s=1}^{T_{0,i}}\operatorname*{\mathrm{blkdiag}}\left(\mathbb{I} _{K},Q^{-1}\right) E\left[ \bar{z}_{it}\bar{z}_{is}^{\prime }u_{it}u_{is} \right] \operatorname*{\mathrm{blkdiag}}\left(\mathbb{I}_{K},(Q^{\prime })^{-1}\right) , \\ \Xi _{\cdot t} &=\lim_{N_{1}\rightarrow \infty }\frac{1}{N_{1}} \sum_{j=1}^{N_{1}}E\left\{ [\nu _{\cdot t,x\lambda }^{\prime }(\mathbf{s} _{\cdot t}^{0})^{-1}(x_{jt}-\pi _{t,x\lambda }^{\prime }\lambda _{j}^{0})-\mu _{\lambda }^{\prime }\Sigma _{\Lambda }^{-1}\lambda _{j}^{0}]u_{jt}\right\} ^{2}, \\ \Xi &=\underset{\left(N,T\right) \rightarrow \infty }{\text{plim}}\frac{1}{ N_{2}}\sum_{i=N_{1}+1}^{N}\alpha _{i\cdot ,T}\frac{1}{T_{0,i}} \sum_{t=1}^{T_{0,i}}\sum_{s=1}^{T_{0,i}}\operatorname*{\mathrm{blkdiag}}\left(\mathbb{I} _{K},Q^{-1}\right) \\ &\quad \times E\left[ \bar{z}_{it}\bar{z}_{is}^{\prime }u_{it}u_{is} \right] \operatorname*{\mathrm{blkdiag}}\left(\mathbb{I}_{K},(Q^{\prime })^{-1}\right) \alpha _{i\cdot ,T}^{\prime }, \end{align*} $$

where ${\mathrm {blkdiag}}\left (\cdot \right ) $ signifies a block diagonal matrix. Let $\underline {T}_{0}=\min _{i\in \left [ N\right ] \backslash \left [ N_{1}\right ] }T_{0,i}$ and $\bar {T}_{1}=\frac {1}{N_{2}}\sum _{i=N_{1}+1}^{N}T_{1,i}.$ We add the following assumption.

Assumption 8 imposes some conditions of $\left \{ x_{it},u_{it},\lambda _{i}^{0},f_{t}^{0}\right \} $ that can be verified under some primitive conditions. Assumption 8 contains three types of asymptotic conditions. First, parts (iv) and (v) impose pointwise convergence conditions. These conditions support parts (i) and (ii), respectively, of Theorem 4.1. Second, part (vi) imposes aggregate convergence conditions jointly over treated units and post-treatment periods and is needed only for the overall average estimator in Theorem 4.1(iii). Third, parts (i)–(iii) and (vii) impose uniform bounds over the relevant unit and time indices. In establishing Theorem 4.1(iii), the aggregate limits in part (vi), together with these uniform bounds, control the accumulation of the leading terms and estimation remainders when averaging over both dimensions.

The following theorem states the asymptotic distributions of $\hat {\tau } _{i\cdot }, \hat {\tau }_{\cdot t}$ , and $\hat {\tau }.$

Theorem 4.1. Suppose Assumptions 1–8 hold for $ \left \{ \left (y_{it},x_{it}\right ) \right \} _{i\in \left [ N_{1}\right ] ,t\in \left [ T\right ] }$ with obvious changes to take into account the change in the sample size. Suppose that $(\sqrt {N_{1}}/T+\sqrt {T}/N_{1})\log ^{2}(N_{1}T)\rightarrow 0$ and $\frac {T_{0,i}}{T}\rightarrow r_{i}\in (0,1]$ as $(N_{1},T)\rightarrow \infty $ . Then

(i) For each fixed treated unit i, as $T_{1,i}\rightarrow \infty ,$

$$ \begin{align*} \sqrt{T_{1,i}}\left(\hat{\tau}_{i\cdot }-\tau _{i\cdot }\right) \overset{d}{ \rightarrow }N\left(0,\sigma _{i\cdot ,u}^{2}+\frac{1-r_{i}}{r_{i}}\Xi _{i\cdot }\right) , \end{align*} $$

where $\Xi _{i\cdot }\equiv (\mu _{i\cdot ,x}^{\prime }\mathbb {S}_{1}+\mu _{f}^{\prime }Q^{-1}\mathbb {S}_{2})S_{i}^{-1}\Omega _{i,zu}S_{i}^{-1}(\mathbb {S}_{1}^{\prime }\mu _{i\cdot ,x}+\mathbb {S}_{2}^{\prime }(Q^{-1})^{\prime }\mu _{f}).$

(ii) For each fixed time period $t\in \left [ T\right ] \backslash \left [ \min _{i\in \left [ N\right ] }T_{0,i}\right ] $ , as $ \left \vert \mathcal {G}_{t}\right \vert \rightarrow \infty $ and $\frac { \left \vert \mathcal {G}_{t}\right \vert }{N_{1}}\rightarrow r_{\cdot t}\in \lbrack 0,\infty ),$

$$ \begin{align*} \sqrt{\left\vert \mathcal{G}_{t}\right\vert }\left(\hat{\tau}_{\cdot t}-\tau _{\cdot t}\right) \overset{d}{\rightarrow }N\left(0,\sigma _{\cdot t,u}^{2}+r_{\cdot t}\Xi _{\cdot t}\right). \end{align*} $$

(iii) As $\bar {T}_{1}\rightarrow \infty $ and $N_{2}$ is either fixed or diverges to infinity such that $N_{2}=o(\min (N_{1}, \underline {T} _{0}/\log N, \delta _{N_{1}T}^{4}/\{\bar {T}_{1}[\log (N_{1}T)]^{2}\})),$

$$ \begin{align*} \sqrt{N_{2}\bar{T}_{1}}\left(\hat{\tau}-\bar{\tau}\right) \overset{d}{ \rightarrow }N\left(0,\sigma _{u}^{2}+\Xi \right). \end{align*} $$

Theorem 4.1 studies the asymptotic distributions of $\hat {\tau } _{i\cdot }, \hat {\tau }_{\cdot t}$ , and $\hat {\tau }.$ To make inferences, we need to consistently estimate the asymptotic variances. Let $\hat {z} _{it}=(x_{it}^{\prime },\hat {f}_{t}^{\prime })^{\prime }$ for $(i,t)\in [N]\times [T]$ and $\hat {u}_{it}=y_{it}-x_{it}^{ \prime }(\hat {\beta } _{i\cdot }+\hat {\beta }_{\cdot t})-\hat {\lambda } _{i}^{\prime }\hat {f}_{t}$ for $(i,t)\in \mathcal {D}^0$ .

Case (1). Since $\sigma _{i\cdot ,u}^{2}$ and $\Omega _{i,zu}$ denote the long-run variance of $\left \{ u_{it}\right \} $ and $ \left \{ z_{it}u_{it}\right \} ,$ it is standard to use $\left \{ \hat {u} _{it}\right \} $ and $\left \{ \hat {z}_{it}\hat {u}_{it}\right \} $ to estimate them by HAC estimators. The main difficulty here lies in the fact that the above definition of $\hat {u}_{it}$ only works for the pre-treatment period when individual i receives the treatment in period $t.$ When $(i,t)\in \mathcal { D}^{1}$ , we cannot consistently estimate $u_{it}$ . Under the covariance stationarity of $\left \{ \left (u_{it},\bar {z}_{it}u_{it}\right ) ,t\geq 1\right \} ,$ we can estimate $\sigma _{i\cdot ,u}^{2}$ and $\Omega _{i,zu}$ , respectively, by

$$ \begin{align*} \hat{\sigma}_{i\cdot ,u}^{2} &=\frac{1}{T_{0,i}}\left[ \sum_{t=1}^{T_{0,i}} \hat{u}_{it}^{2}+2\sum_{r=1}^{T_{0,i}-1}w\left(r/m\right) \sum_{t=r+1}^{T_{0,i}}\hat{u}_{it}\hat{u}_{i,t-r}\right] \ \text{and} \\ \hat{\Omega}_{i,zu} &=\frac{1}{T_{0,i}}\sum_{t=1}^{T_{0,i}}\hat{z}_{it}\hat{ z}_{it}^{\prime }\hat{u}_{it}^{2}+\sum_{r=1}^{T_{0,i}-1}w\left(r/m\right) \sum_{t=r+1}^{T_{0,i}}\left(\hat{z}_{it}\hat{z}_{i,t-r}^{\prime }+\hat{z} _{i,t-r}\hat{z}_{it}^{\prime }\right) \hat{u}_{it}\hat{u}_{i,t-r}, \end{align*} $$

where $w\left (\cdot \right ) $ is a kernel function and $m=m_{T}$ is a truncation parameter. Of course, if $\left \{ (u_{it},\bar {z} _{it}u_{it})\right \} $ share the same variance–covariance matrix across $ i\in \left [ N\right ] \backslash \left [ N_{1}\right ] ,$ one can pool over the treated cross-sectional units to estimate the above variances more accurately in finite samples. Then, we can estimate the asymptotic variance of $\sqrt {T_{1,i}}\left (\hat {\tau }_{i\cdot }-\tau _{i\cdot }\right ) $ by

$$ \begin{align*} \widehat{\text{avar}}(\sqrt{T_{1,i}}\hat{\tau}_{i\cdot })=\hat{\sigma} _{i\cdot ,u}^{2}+\frac{T_{1,i}}{T_{0,i}}(\hat{\mu}_{i\cdot ,x}^{\prime } \mathbb{S}_{1}+\hat{\mu}_{f}^{\prime }\mathbb{S}_{2})\hat{S}_{i,T_{0,i}}^{-1} \hat{\Omega}_{i,zu}\hat{S}_{i,T_{0,i}}^{-1}(\mathbb{S}_{1}^{\prime }\hat{\mu} _{i\cdot ,x}+\mathbb{S}_{2}^{\prime }\hat{\mu}_{f}), \end{align*} $$

where $\hat {\mu }_{i\cdot ,x}=\frac {1}{T_{1,i}}\sum _{t=T_{0,i}+1}^{T}x_{it},$ $\hat {\mu }_{f}=\frac {1}{T}\sum _{t=1}^{T}\hat {f}_{t},$ and $\hat {S} _{i,T_{0,i}}=\frac {1}{T_{0,i}}\sum _{t=1}^{T_{0,i}}\hat {z}_{it}\hat {z} _{it}^{\prime }.$ Alternatively, under the covariance stationarity of $ \left \{ \bar {z}_{it},t\geq 1\right \} ,$ we can also use the following estimators $\tilde {\mu }_{i\cdot ,x}=\frac {1}{T} \sum _{t=1}^{T}x_{it}$ and $\tilde {S}_{i}=\frac {1}{T}\sum _{t=1}^{T}\hat {z} _{it}\hat {z}_{it}^{\prime }$ to replace $\hat {\mu }_{i\cdot ,x}$ and $\hat {S} _{i,T_{0,i}},$ respectively.

Case (2). Under the condition that $\left \{ u_{it},x_{it},\lambda _{i}^{0}\right \} $ are independent across $i,$ we can estimate $\sigma _{\cdot t,u}^{2}$ and $\Xi _{\cdot t}$ , respectively, by

$$ \begin{align*} \hat{\sigma}_{\cdot t,u}^{2} &=\frac{1}{\left\vert [N]\backslash \mathcal{G} _{t}\right\vert }\sum_{i\in \lbrack N]\backslash \mathcal{G}_{t}}\hat{u} _{it}^{2}\text{ and } \\ \hat{\Xi}_{\cdot t} &=\frac{1}{N_{1}}\sum_{j=1}^{N_{1}}[\hat{\nu}_{\cdot t,x\lambda }^{\prime }\left(\boldsymbol{\hat{s}}_{\cdot t}\right) ^{-1}(x_{jt}- \mathbf{x}_{\cdot t}^{\prime }\hat{\Lambda}(\hat{\Lambda}^{\prime }\hat{ \Lambda})^{-1}\hat{\lambda}_{j})-\hat{\mu}_{\lambda }^{\prime }(N_{1}^{-1} \hat{\Lambda}^{\prime }\hat{\Lambda})^{-1}\hat{\lambda}_{j}]^{2}\hat{u} _{jt}^{2}, \end{align*} $$

where $\hat{\boldsymbol{s}}_{\cdot t}=\frac {1}{N_{1}}\mathbf {x}_{\cdot t}^{\prime }\mathbb {M}_{\hat {\Lambda }}\mathbf {x}_{\cdot t}, \hat {\mu } _{\lambda }=\frac {1}{N_{1}}\sum _{j=1}^{N_{1}}\hat {\lambda }_{j}$ , and $\hat { \nu }_{\cdot t,x\lambda }=\frac {1}{\left \vert \mathcal {G}_{t}\right \vert } \sum _{i\in \mathcal {G}_{t}}[x_{it}-$ $\mathbf {x}_{\cdot t}^{\prime }\hat {\Lambda }(\hat {\Lambda }^{\prime }\hat {\Lambda })^{-1}\hat {\lambda }_{i}].$ Then we can estimate the asymptotic variance of $\sqrt {\left \vert \mathcal {G} _{t}\right \vert }\left (\hat {\tau }_{\cdot t}-\tau _{\cdot t}\right ) $ by

$$ \begin{align*} \widehat{\text{avar}}(\sqrt{\left\vert \mathcal{G}_{t}\right\vert }\hat{\tau} _{\cdot t})=\hat{\sigma}_{\cdot t,u}^{2}+\frac{\left\vert \mathcal{G} _{t}\right\vert }{N_{1}}\hat{\Xi}_{\cdot t}. \end{align*} $$

Case (3). We can estimate $\sigma _{u}^{2}$ and $\Xi $ by the panel HAC estimators $\hat {\sigma }_{u}^{2}$ and $\hat {\Xi },$ respectively, where

$$ \begin{align*} \hat{\sigma}_{u}^{2} &=\frac{1}{N_{2}\bar{T}_{0}} \sum_{i=N_{1}+1}^{N}T_{0,i} \hat{\sigma}_{i\cdot ,u}^{2},\text{ }\hat{\Xi}= \frac{1}{N_{2}} \sum_{i=N_{1}+1}^{N}\hat{\alpha}_{i\cdot ,T}\hat{\Omega} _{i,zu}\hat{\alpha} _{i\cdot ,T}^{\prime }, \\ \bar{T}_{0} &=\frac{1}{N_{2}}\sum_{i=N_{1}+1}^{N}T_{0,i},\text{ and }\hat{ \alpha}_{i\cdot ,T}=\frac{T_{1,i}}{\sqrt{T_{0,i}\bar{T}_{1}}}(\hat{\mu} _{i\cdot ,x}^{\prime }\mathbb{S}_{1}+\hat{\mu}_{f}^{\prime }\mathbb{S}_{2}) \hat{S}_{i,T_{0,i}}^{-1}. \end{align*} $$

As in Case (1), here we use the covariance stationarity of $\left \{ \left (u_{it},\bar {z}_{it}u_{it}\right ) ,t\geq 1\right \} $ to obtain $\hat {\sigma } _{u}^{2}.$ Then we can estimate the asymptotic variance of $\sqrt {N_{2}\bar { T }_{1}}\left (\hat {\tau }-\bar {\tau }\right ) $ by

$$ \begin{align*} \widehat{\text{avar}}(\sqrt{N_{2}\bar{T}_{1}}\hat{\tau})=\hat{\sigma} _{u}^{2}+\hat{\Xi}. \end{align*} $$

Despite the fact that $\hat {\tau }_{it}$ is generally inconsistent, we can still construct valid CIs for $\tau _{it}$ under certain stationarity conditions on $\left \{ u_{it}\right \} .$ This is because $\hat {\tau }_{it}=\tau _{it}+u_{it}+o_{P}\left (1\right ) $ by Lemma A.7(i) in Section A.2 of the Supplementary Material. For example, if $ \left \{ u_{it}\right \} $ for the treated units are strictly stationary over both i and t, we can establish the 100 $\left (1-\alpha \right ) \%$ CI for $\tau _{it}$ as $\left [ \hat {\tau }_{it}+\hat {q} _{u,\alpha /2},\hat {\tau }_{it}+\hat {q}_{u,1-\alpha /2}\right ] ,$ where, for example, $\hat {q}_{u,\alpha /2}$ denotes the $\alpha /2$ -sample quantile of $\{- \hat {u}_{it}, i\in \left [ N\right ] \backslash \left [ N_{1}\right ] , t\in \left [ T_{0,i}\right ] \}.$ Similarly, if we further impose strict stationarity for $\left \{ u_{it}\right \} $ for both the control units and treated units, we can use the sample quantiles of the combination of $\{- \hat {u}_{it},\,i\in \left [ N\right ] \backslash \left [ N_{1}\right ] , t\in \left [ T_{0,i}\right ] \}$ and $\{-\hat {u}_{it},i\in \left [ N_{1}\right ] ,t\in \left [ T\right ] \}.$ Alternatively, we can apply the permutation test proposed by Chernozhukov et al. (Reference Chernozhukov, Wüthrich and Zhu2021) to conduct inference for $\tau _{it} $ .

5.1 Data-Generating Processes and Implementations

To examine the finite sample performance of our estimation method, we conduct simulations with four DGPs. For all DGPs, the untreated outcome is generated as in models (1.1) and (1.2), and the treatment effect (TE), $\tau _{it}$ is randomly drawn from the $N\left (1,1\right ) $ distribution.

We first consider simple DGPs 1 and 2 to illustrate the main idea of our approach. In both DGPs, $x_{it}$ is a scalar, the number of factors ( $R^{0}$ ) is 2, $x_{it}$ and $u_{it}$ are all generated from independent and identically distributed (IID) $N\left (0,1\right ) ,$ and $\lambda _{i}^{0}$ and $f_{t}^{0}$ are drawn from IID $\ N\left (0,\mathbb {I}_{2}\right ) .$ We consider two different cases for $\beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}.$ In DGP 1 (heterogeneous case), $\beta _{i\cdot }^{0}$ and $\beta _{\cdot t}^{0}$ are IID $N\left (0,1\right ) $ and $N\left (1,1\right ) $ random variables. In DGP 2 (homogeneous case), $\beta _{i\cdot }^{0}=0$ and $ \beta _{\cdot t}^{0}=1$ for all i and $t.$ Note that for DGP 2, our proposed model is over-fitted, therefore we can examine the effect of over-fitting by comparing it with the traditional homogeneous model where $ (\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0})$ is restricted to be a constant.

DGPs 3 and 4 are more realistic. In both DGPs, there are two covariates $ x_{it}=\left (x_{1,it},x_{2,it}\right ) ^{\prime }$ and two factors $ f_{t}^{0}=\left (f_{t,1}^{0},f_{t,2}^{0}\right ) ^{\prime }$ . Correspondingly, define $\lambda _{i}^{0}=\left (\lambda _{i,1}^{0},\lambda _{i,2}^{0}\right ) ^{\prime }, \beta _{i\cdot }^{0}=\left (\beta _{1,i\cdot }^{0},\beta _{2,i\cdot }^{0}\right ) ^{\prime }$ , and $\beta _{\cdot t}^{0}=\left (\beta _{1,\cdot t}^{0},\beta _{2,\cdot t}^{0}\right ) ^{\prime }.$ The details of the specifications are as follows:

$$ \begin{align*} f_{t,1}^{0} &=\rho f_{t-1,1}^{0}+\sqrt{1-\rho ^{2}}v_{t,1},\quad f_{t,2}^{0}=\rho f_{t-1,2}^{0}+\sqrt{1-\rho ^{2}}v_{t,2}, \\ x_{1,it} &=0.25+0.25\lambda _{i}^{0\prime }f_{t}^{0}+\eta _{it},\quad \eta _{it}=\rho \eta _{i,t-1}+\sqrt{1-\rho ^{2}}v_{it,3},\quad x_{2,it}=0.5v_{it,4}, \\ u_{it} &=\sqrt{1+0.1x_{1,it}^{2}}\varepsilon _{it},\quad \text{ and } \varepsilon _{it}=\rho \varepsilon _{i,t-1}+\sqrt{1-\rho ^{2}}v_{it,5}, \end{align*} $$

$v_{t,1}, v_{t,2}, v_{it,3}, v_{it,4}, \lambda _{i,1}^{0},$ and $ \lambda _{i,2}^{0}$ are drawn from IID $N\left (0,\mathbb {I}_{6}\right ) $ and $v_{it,5}$ is IID random variable from the distribution of $\chi ^{2}\left (1\right ) -1. \beta _{1,i\cdot }^{0}$ and $\beta _{1,\cdot t}^{0} $ are IID $N\left (0,1\right ) $ and $N\left (1,1\right ) $ random variables, and $\beta _{2,i\cdot }^{0}=0$ and $\beta _{2,\cdot t}^{0}=1$ for all i and $t.$ Therefore, the slope coefficient of $x_{1,it}$ is heterogeneous and that of $x_{2,it}$ is homogeneous. The difference between DGPs 3 and 4 lies in $\rho .$ In DGP 3, $\rho =0,$ therefore, there is no serial correlation in $x_{1,it}, f_{t,1}^{0}, f_{t,2}^{0}$ , and $u_{it}.$ In DGP 4, $\rho =1/3,$ so $x_{1,it}, f_{t,1}^{0}, f_{t,2}^{0}$ , and $ u_{it}$ are all serially correlated. Note that DGP 4 is the most complicated where the error term is conditionally heteroskedastic and serially correlated.

For DGPs 1–3, we use the White-type estimator for the variance terms to conduct inference. To improve the finite-sample performance, we impose the condition that the treated data are stationary along both cross-sectional and time dimensions. For DGP 4, we use an HAC-type estimator for the variance. It is well known that the performance of the standard HAC estimator can be poor when the time-series dimension is short (e.g., see the survey paper by Lazarus et al., Reference Lazarus, Lewis, Stock and Watson2018). To partially address this issue, we impose the stationarity of both treated and control data, and also apply a pre-whitening (PW) procedure as in Andrews and Monahan (Reference Andrews and Monahan1992) and Lu and Su (Reference Lu and Su2023). We use the growth ratio method of Ahn and Horenstein (Reference Ahn and Horenstein2013) to determine the number of factors in the initial nuclear norm regularized estimation procedure, as it outperforms singular value thresholding in our simulations and does not depend on any tuning parameter other than the maximum number of factors allowed. We set the tuning parameter as $\gamma _{NT}=0.5\max (N_{1}^{-1/2},T^{-1/2}),$ following Miao et al. (Reference Miao, Phillips and Su2023). We consider different combinations of $N_{1}$ and T: $\left (N_{1},T\right ) =\left (40,40\right ) , \left (40,80\right ) , \left (80,40\right ) , \left (80,80\right ), $ and $\left (120,120\right ) .$ We let $N_{2}=0.25N_{1}$ and $T_{0i}=T_{1i}=0.5T.$ The number of replications is 250.

5.2 Simulation Results

We first report the performance of our key Step 1 in Algorithm 2, where we use the data of control units to estimate $\{\beta _{i\cdot }^{0},\lambda _{i}^{0}\}_{i\in \left [ N_{1}\right ] }$ and $\{\beta _{\cdot t}^{0},f_{t}^{0}\}_{t\in \left [ T\right ] }.$ We define the true slope coefficient as $\beta _{it}^{0}=\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}$ and the true factor component as $\theta _{it}^{0}=\lambda _{i}^{0\prime }f_{t}^{0}.$ Correspondingly, define their estimators as $\hat {\beta }_{it}= \hat {\beta }_{i\cdot }+\hat {\beta }_{\cdot t}$ and $\hat {\theta }_{it}=\hat { \lambda }_{i}^{\prime }\hat {f}_{t},$ respectively. For ease of reporting, we only present the results for the slope coefficient of the first covariates for DGPs 3 and 4, where there are two covariates. Table 1 reports the mean squared errors (MSEs) of $\hat {\beta }_{1,it}$ (the first element of $\hat {\beta }_{it}$ ) and the correlation between $\hat {\beta }_{1,it}$ and $ \beta _{1,it}^{0}$ .Footnote ⁴ Specifically, we define

$$ \begin{align*} \text{MSE}_{\beta }=\frac{1}{N_{1}T}\sum_{i=1}^{N_{1}}\sum_{t=1}^{T}\left(\hat{\beta}_{1,it}-\beta _{1,it}^{0}\right) ^{2}\text{ and }\rho _{\beta }= \text{ }\widehat{\text{corr}}\left(\hat{\beta}_{1,it},\beta _{1,it}^{0}\right) , \end{align*} $$

and average both measures over 250 replications. Table 1 also reports the MSEs of $\hat {\theta }_{it}$ and the correlation between $\hat { \theta }_{it}$ and $\theta _{it}^{0}$ , which are defined as

$$ \begin{align*} \text{MSE}_{\theta }=\frac{1}{N_{1}T}\sum_{i=1}^{N_{1}}\sum_{t=1}^{T}\left(\hat{\theta}_{it}-\theta _{it}^{0}\right) ^{2}\text{ and }\rho _{\theta }= \widehat{\text{corr}}\left(\hat{\theta}_{it},\theta _{it}^{0}\right), \end{align*} $$

averaging over 250 replications. Table 1 shows that our estimation procedure in Step 1 works well. The correlations between $\hat { \beta }_{1,it}$ and $\beta _{1,it}^{0}$ and between $\hat {\theta }_{it}$ and $ \theta _{it}^{0}$ are all above 0.97 for DGPs 1 and 2, and all above 0.85 for DGPs 3 and 4. The MSEs of $\hat {\beta }_{1,it}$ and $\hat {\theta }_{it}$ decrease with $N_{1}$ and $T,$ as expected. Our method can determine the correct number of factors with a high probability. For simple DGPs 1 and 2, we determine the number of factors with a correct rate of 100% for all sample sizes considered. For more complicated DGPs 3 and 4, the correct rates are all above 90% and close to 100% when the sample size is large.

Table 1

Performance of Step 1 in Algorithm 2

Table 1 Long description

The table is divided into Section A: Performance of beta hat sub 1, i t and Section B: Performance of theta hat sub i t.

Both sections share the same vertical axis of sample sizes (N sub 1, T): (40, 40), (40, 80), (80, 40), (80, 80), and (120, 120).

Horizontal columns for both sections are grouped by D G P 1, D G P 2, D G P 3, and D G P 4. Each D G P contains two metrics: M S E sub beta (or M S E sub theta) multiplied by 100, and Corr: rho sub beta (or rho sub theta).

In Section A:

- For D G P 1, M S E sub beta times 100 decreases from 6.02 to 1.77 as sample size increases; correlation increases from 0.985 to 0.996.

- D G P 2 shows similar M S E values (6.02 to 1.79), but correlation is undefined.

- D G P 3 and D G P 4 show higher M S E values, starting at 16.33 and 17.53 respectively for (40, 40), and decreasing to 4.42 and 4.98 for (120, 120).

In Section B:

- M S E sub theta values are generally higher than M S E sub beta.

- For D G P 1, M S E sub theta times 100 decreases from 10.78 to 3.39; correlation increases from 0.975 to 0.992.

- D G P 3 and D G P 4 show the highest error rates, with M S E sub theta times 100 reaching 33.12 and 37.64 respectively at the smallest sample size (40, 40).

Note: Corr stands for correlation. For the homogeneous DGP 2, $\rho _{\beta }$ (the correlation between $\hat {\beta }_{it}$ and $\beta _{it}^{0})$ is not defined.

We next examine the performance of our estimator of various treatment effects, namely, $\hat {\tau }_{it}, \hat {\tau }_{i\cdot }, \hat {\tau } _{\cdot t}$ , and $\hat {\tau }.$ Table 2 reports their MSEs, which are defined as

$$ \begin{align*} \text{MSE}_{\hat{\tau}_{it}} &=\frac{1}{N_{2}T_{1}}\sum_{i=N_{1}+1}^{N} \sum_{t=T_{0,i}+1}^{T}\left(\hat{\tau}_{it}-\tau _{it}\right) ^{2},\text{ MSE}_{\hat{\tau}_{i\cdot }}=\frac{1}{N_{2}}\sum_{i=N_{1}+1}^{N}\left(\hat{ \tau}_{i\cdot }-\tau _{i\cdot }\right) ^{2},\text{ } \\ \text{MSE}_{\hat{\tau}_{\cdot t}} &=\frac{1}{\bar{T}_{1}} \sum_{t=T_{0}+1}^{T}\left(\hat{\tau}_{\cdot t}-\tau _{\cdot t}\right) ^{2} \text{ and MSE}_{\hat{\tau}}=\left(\hat{\tau}-\bar{\tau}\right) ^{2}, \end{align*} $$

averaging over 250 replications. Note that $\hat {\tau }_{it}$ is not a consistent estimator of $\tau _{it},$ and MSE $_{\hat {\tau }_{it}}$ converges to Var( $u_{it}$ ) in probability. MSE $_{\hat {\tau }_{i\cdot }}$ decreases at the rate $T_{1,i}^{-1}$ , MSE $_{\hat {\tau }_{\cdot t}}$ decreases at the rate $ \left \vert \mathcal {G}_{t}\right \vert ^{-1}$ , and MSE $_{\hat {\tau }}$ decreases at the rate of $(N_{2}\bar {T}_{1})^{-1}.$ Our simulation results confirm all these theoretical results. We also compare the four MSE measures between our heterogeneous (hetero) and the conventional homogeneous (homo) models which assume that the slope coefficient $\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}$ is a constant. Consider the simple DGPs 1 and 2 first. Note that the homo model is misspecified for DGP 1 and correctly specified for DGP 2. For DGP 1, our hetero model performs much better than the homo model, as expected. For example, MSE $_{\hat {\tau }_{it}}$ based on our hetero model and that based on the homo models are 1.08 and 3.17, respectively, when $ \left (N_{1},T\right ) =\left (120,120\right ) .$ The large difference is easy to understand. For DGP 1, the homo model effectively estimates

$$ \begin{align*} y_{it}=\beta _{0}^{0\prime }x_{it}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it}^{o},\text{ where }u_{it}^{o}=(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}-\beta _{0}^{0})^{\prime }x_{it}+u_{it}\text{ and }\beta _{0}^{0}=E(\beta _{\cdot t}^{0}). \end{align*} $$

That is, the homo model treats $\left (\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}-\beta _{0}^{0}\right ) ^{\prime }x_{it}$ as part of the unpredictable error, which is actually predictable. In contrast, our hetero model can capture all of the predictable components. Note that for this DGP, Var $ \left (u_{it}\right ) =1$ and Var $\left (u_{it}^{o}\right ) =3,$ which explains the large difference in MSE $_{\hat {\tau }_{it}}.$ For the same reason, MSE $_{\hat {\tau }_{i\cdot }},$ MSE $_{\hat {\tau }_{\cdot t}}$ , and MSE $_{ \hat {\tau }}$ based on our hetero models are much smaller than those based on the homo models. For DGP 2, where the true DGP is homogeneous, we find that the performance of our hetero model is similar to that of the homo model especially when the sample size is large, despite the fact that our hetero model is vastly over-fitted. For example, MSE $_{\hat {\tau }_{it}}$ based on our hetero model and the homo model are 1.08 and 1.05, respectively, when $ \left (N_{1},T\right ) =\left (120,120\right ) .$ This comparison in the simulations suggests that even when the true DGP is homogeneous (DGP 2), the estimation efficiency loss based on the overfitted hetero model is minimal. Conversely, when the true DGP is heterogeneous (DGP 1), our hetero model significantly outperforms the homo model. For DGPs 3 and 4, the performance of our hetero model is much better than the homo model. We also compare with the SCM and find that our hetero model outperforms SCM substantially for all DGPs and all sample sizes.Footnote ⁵ In our setup, SCM can only predict the factor components. Consider the simple DGPs 1 and 2 where the covariates $x_{it}$ is cross-sectionally independent. SCM cannot predict $\left (\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0}\right ) ^{\prime }x_{it}$ by taking the average of untreated outcomes. Therefore, SCM is essentially based on the model

$$ \begin{align*} y_{it}=c_{0}+\lambda _{i}^{0\prime }f_{t}^{0}+u_{it}^{\dagger }, &\text{ where }u_{it}^{\dagger }=[(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0})^{\prime }x_{it}-c_{0}]+u_{it}\text{ and } \\ & c_{0}=E[(\beta _{i\cdot }^{0}+\beta _{\cdot t}^{0})^{\prime }x_{it}], \end{align*} $$

and treats $u_{it}^{\dagger }$ as an unpredictive error. MSE $_{\hat {\tau } _{it}}$ based on SCM is at least as large as Var $(u_{it}^{\dagger }),$ which is equals to $4$ and $2$ for DGPs 1 and 2, respectively. The main advantage of our approach is that we are able to utilize the information of covariates. In the Supplementary Material, we compare our approach with SCM for pure factor structures without any covariates and show that both approaches work well. We emphasize that the superior performance of our approach crucially depends on the correct specification of our model. In the Supplementary Material, we design DGPs with grouped factor structure and show that SCM can outperform our method in this special case.

Table 2

MSE of estimators of TE and ATT

Table 2 Long description

A complex data table titled M S E of estimators of T E and A T T.

Header Structure:

* The top row identifies four Data Generating Processes: D G P 1, D G P 2, D G P 3, and D G P 4.

* Below each D G P are three sub-columns: hetero (heterogeneous model), homo (constant slope model), and synth (synthetic control method).

* The leftmost column lists parameter sets in the format (N sub 1, T, N sub 2, T sub 0 i, T sub 1 i).

Section A: M S E of tau-hat sub i t times 100 (T E for treated individuals at a given time):

* Values range from 105 to 835.

* For (40, 40, 10, 20, 20), D G P 1 hetero is 129, homo is 364, and synth is 592.

* For (120, 120, 30, 60, 60), D G P 4 hetero is 253, homo is 464, and synth is 695.

Section B: M S E of tau-hat sub i dot times 100 (A T T for individuals):

* Values range from 1.72 to 81.5.

* For (40, 40, 10, 20, 20), D G P 1 hetero is 6.42, homo is 17.9, and synth is 29.0.

* For (120, 120, 30, 60, 60), D G P 4 hetero is 8.68, homo is 19.6, and synth is 25.5.

Section C: M S E of tau-hat sub dot t times 100 (A T T for time periods):

* Values range from 3.44 to 110.

* For (40, 40, 10, 20, 20), D G P 1 hetero is 13.6, homo is 38.6, and synth is 78.1.

* For (120, 120, 30, 60, 60), D G P 4 hetero is 8.58, homo is 21.4, and synth is 31.8.

Section D: M S E of tau-hat times 100 (Aggregate A T T):

* Values range from 0.06 to 12.4.

* For (40, 40, 10, 20, 20), D G P 1 hetero is 0.76, homo is 1.84, and synth is 3.75.

* For (120, 120, 30, 60, 60), D G P 4 hetero is 0.30, homo is 0.77, and synth is 1.48.

General Trend: Across all sections and D G P s, the hetero estimator consistently yields the lowest M S E, followed by homo, with synth typically yielding the highest M S E.

Note: hetero refers to the heterogeneous model proposed in this article, homo refers to the model where the slope coefficients are imposed to be constant, and synth refers to SCM (synthetic control method).

Table 3 presents the actual coverage of the 95% CI for the TE, $ \tau _{it}$ and ATT, $\tau _{i\cdot }, \tau _{\cdot t},$ and $\tau .$ The CI of $\tau _{it}$ is based on the quantiles of the residuals and the CIs of $\tau _{i\cdot }, \tau _{\cdot t},$ and $\tau $ are based on our asymptotic normality result in Theorem 4.1. The actual coverages are slightly under-covered when the sample size is small, however they are all close to 95% when the sample size is relatively large. This suggests the validity of our inference procedure.

Table 3

Empirical coverage of 95% confidence intervals

Table 3 Long description

The table is organized into 17 columns. The first column lists input coordinates in the format (N sub 1, T, N sub 2, T sub 0 i, T sub 1 i). The remaining 16 columns are grouped under four main headers: D G P 1, D G P 2, D G P 3, and D G P 4. Each D G P group contains four sub-columns for the parameters tau sub i t, tau sub i dot, tau sub dot t, and tau.

Data rows by input coordinates:

* (40, 40, 10, 20, 20): D G P 1 (0.90, 0.93, 0.89, 0.91); D G P 2 (0.90, 0.92, 0.90, 0.92); D G P 3 (0.88, 0.92, 0.87, 0.88); D G P 4 (0.87, 0.87, 0.86, 0.85).

* (40, 80, 10, 40, 40): D G P 1 (0.93, 0.95, 0.92, 0.95); D G P 2 (0.93, 0.94, 0.92, 0.94); D G P 3 (0.92, 0.93, 0.91, 0.92); D G P 4 (0.91, 0.89, 0.90, 0.89).

* (80, 40, 20, 20, 20): D G P 1 (0.90, 0.92, 0.90, 0.92); D G P 2 (0.90, 0.93, 0.90, 0.95); D G P 3 (0.89, 0.92, 0.89, 0.93); D G P 4 (0.88, 0.89, 0.88, 0.91).

* (80, 80, 20, 40, 40): D G P 1 (0.93, 0.94, 0.93, 0.92); D G P 2 (0.93, 0.94, 0.93, 0.96); D G P 3 (0.92, 0.94, 0.92, 0.96); D G P 4 (0.91, 0.92, 0.92, 0.93).

* (120, 120, 30, 60, 60): D G P 1 (0.94, 0.94, 0.93, 0.95); D G P 2 (0.94, 0.95, 0.94, 0.94); D G P 3 (0.93, 0.94, 0.93, 0.94); D G P 4 (0.93, 0.92, 0.93, 0.94).

In Section S.7.3 of our Supplementary Material, we show that when $\beta _{\cdot t}$ is nonstationary, the permutation test proposed in Chernozhukov et al. (Reference Chernozhukov, Wüthrich and Zhu2021) applied to SCM may have severe size distortion. The nonstationarity of $ \beta _{\cdot t}$ causes nonstationarity of the error term in SCM, which invalidates many existing inference procedures, such as Chernozhukov et al. (Reference Chernozhukov, Wüthrich and Zhu2021), Masini and Medeiros (Reference Masini and Medeiros2021), and Masini and Medeiros (Reference Masini and Medeiros2022). Our approach estimates all $\beta _{\cdot t}$ as unknown parameters, and therefore is immune to this problem.