Counterfactuals in factor models

Counterfactuals in factor models

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\it Keywords:} counterfactuals, average marginal effects, panel data, factor models, high-dimensional data

{\it JEL codes:} C32, C33, C38 \spacingset{1.5}

Introduction

The goal is to infer the causal effect of a possibly continuous treatment variable $d_{it}\in{\mathbb{R}}$ on an outcome variable in a panel data setting. Let $y_{it}(d)$ be the potential outcome under treatment status $d$ of unit $i\in\{1,\dots,N\}$ at date $t\in\{1,\dots,T\}$. The realized outcome of unit $i$ at time $t$ is given by $y_{it}=y_{it}(d_{it})$. Our main parameters of interest are the unit-specific and time-specific average marginal effects (AME) of the treatment. The unit-specific AME of unit $i$ is $$\Delta_i ={\mathbb{E}}_t[y_{it}'(d_{it})],$$ where $y_{it}'(\cdot)$ is the derivative of $y_{it}(\cdot)$ and ${\mathbb{E}}_t$ denotes the expectation over $t$, for fixed $i$ formally, ${\mathbb{E}}_t[y_{it}'(d_{it})]=\operatorname*{plim}\limits_{T\to\infty}T^{-1}\sum_{t=1}^T y_{it}'(d_{it}) $. Note that, when we introduce derivatives or probability limits, we implicitly assume that they exist. Similarly, the time-specific AME at date $t$ is defined by $$\Delta_t={\mathbb{E}}_i[y_{it}'(d_{it})],$$ where ${\mathbb{E}}_i$ denotes the expectation over $i$, for fixed $t$, formally, ${\mathbb{E}}_i[y_{it}'(d_{it})]=\operatorname*{plim}\limits_{N\to\infty}N^{-1}\sum_{i=1}^N y_{it}'(d_{it}).$ We are also interested in the AME over the whole population, that is $$\Delta={\mathbb{E}}[\Delta_i]={\mathbb{E}}[y_{it}'(d_{it})].$$ Other important parameters such as average counterfactuals and treatment effects are also discussed in the paper.

To estimate the AMEs, we assume that the potential outcomes are related through a set of random factors $f_t\in{\mathbb{R}}^R$ such that

where $R$ is the number of factors, $\lambda_{i*}(d)\in{\mathbb{R}}^R$ is a vector of nonrandom loadings and $\tilde u_{it}(d)\in{\mathbb{R}}$ is a random error term. We have

where $\lambda_{*i}'(\cdot)$ is the derivative of $\lambda_{*i}(\cdot)$, so that knowledge of $\lambda_{i*}(\cdot)$ and $f_t$ is sufficient to identify $\Delta$. When $d_{it}$ is continuous, we only observe $y_{it}(d)$ for each $d$ with probability $0$, so that $\lambda_{i*}(\cdot)$ and $f_t$ cannot be estimated from data on $y_{it}$. To overcome this issue, we leverage a panel $\{x_{\ell t}\in{\mathbb{R}}, \ell=1,\dots,L,\ t=1,\dots,T\}$ that loads on the factors. That is, there exists (deterministic) loadings $\lambda_{\ell}\in{\mathbb{R}}^R$ and random error terms $e_{\ell t}\in{\mathbb{R}}$ such that

the factors can then be estimated from this panel using principal components analysis (or any other method such as cross-section averages). The loadings $\lambda_{*i}(\cdot)$ can then be estimated by linear regression of $y_{it}$ on the estimated factors interacted with powers of $d_{it}$ and additional controls.

The proposed approach allows for treatment effects heterogenous across $i$ and $t$ and leverages the very common modelling assumption of factor models. This set-up can model high-dimensional time series when $N=1$ and both $L$ (number of variables in the time series) and $T$ are large. For instance, in a macroeconomic application, $d_{1t}$ could correspond to the US federal funds rate, while $y_{1t}$ is the unemployment rate and the panel $\{x_{\ell t}\}_{\ell,t}$ corresponds to the well-known FRED-MD dataset of mccracken2016fred. Another interesting context corresponding to our set-up is that of a large $N,T$ panel. For instance, suppose that we observe a set of $N$ countries at $T$ dates. On top of the outcome and the treatment, for each country we observe $K$ variables at each of the $T$ dates. These variables then constitute the panel $\{x_{\ell t}\}_{\ell, t}$, with $L=KN$. \\

Contributions. As mentioned above, we formulate an estimation procedure for $\Delta_i, \Delta_t$, and $\Delta$. Initially, we estimate the factors from the panel $\{x_{\ell t}\}_{\ell, t}$ using principal components analysis (PCA). Subsequently, we perform a linear regression of $y_{it}$ on the estimated factors interacted with powers of $d_{it}$ and additional controls. This approach enables the estimation of $\lambda_{i*}(\cdot)$ and facilitates the recovery of the AMEs through plug-in methods. We establish the asymptotic properties of the estimator $\widehat{\Delta}_i$ for $\Delta_i$ in the high-dimensional time series case, where $T, L \to \infty$ while $N$ can remain fixed. Additionally, we demonstrate the asymptotic normality of the estimators $\widehat{\Delta}_t$ of $\Delta_t$ and $\widehat{\Delta}$ of $\Delta$ in the large $N, T$ panel case, where $T,N,L\to \infty$. The finite sample performance of the estimators is evaluated through simulations. An empirical application on income inequality using a large $N,T$ panel illustrates our novel methodology.

This paper introduces a comprehensive approach to model counterfactuals using a factor model. While our focus in this work centers on estimating AMEs due to their importance in nonlinear panel data models fernandez2018fixed, our methodology is versatile. It readily extends to estimating various average counterfactuals such as ${\mathbb{E}}_i[y_{it}(d)], {\mathbb{E}}_t[y_{it}(d)]$, or ${\mathbb{E}}[y_{it}(d)]$ for fixed $d$, as well as implementing instrumental variables estimation. Section (ref) delves into these extensions and provides corresponding estimators. To maintain conciseness, we refrain from explicitly deriving the asymptotic distributions of these additional estimators. However, similar proof techniques to those employed in establishing the asymptotic normality of the estimators of the average marginal effects can be applied. \\

Related literature. Modeling potential outcomes of a binary treatment using a factor structure is a well-established approach in the literature gobillon2016regional,athey2021matrix,bai2021matrix,fan2022we. These papers use a factor model for the untreated potential outcome but do not model the “treated” counterfactual. Specifically, athey2021matrix,bai2021matrix delve into the realm of causal matrix completion, employing a panel denoted as $\{y_{it}\}_{i,t}$. Given the binary nature of the treatment in their study, they are able to observe numerous entries of the panel $\{y_{it}(0)\}_{i,t}$. By assuming a factor structure on $y_{it}(0)$, they successfully recover $y_{it}(0)$ for all pairs of $i$ and $t$ by effectively “completing" the matrix using the information from $\{y_{it}(0)\}_{i,t}$. However, it is crucial to note that this strategy is not applicable in our current study, which involves a continuous treatment. In this scenario, with probability $1$, none of the entries in the matrix $\{y_{it}(d)\}_{i,t}$ for a given treatment level $d$ will be observed, making it impossible to successfully complete the matrix (the same type of issues would arise with the approaches of gobillon2016regional, fan2022we, since we would also almost never observe the untreated potential outcome). To address this challenge, we propose (i) assuming that the factors are shared across all potential outcomes and (ii) leveraging their learnability from the panel $\{x_{\ell t}\}_{\ell, t}$.

Another related strand of literature is that of panel data models with interactive fixed effects pesaran2006estimation,bai2009panel,greenaway2012asymptotic. In this set-up it is assumed that $y_{it}=\beta_d d_{it}+\beta_x^\top \tilde x_{it}+\lambda_i^\top f_t+\epsilon_{it}$, where $\tilde x_{it}$ is a set of controls, $\lambda_i$ are nonrandom loadings and $\epsilon_{it}$ is an error term. Similarly to us, some approaches also assume that $\tilde x_{it}$ follows a factor model with $f_t$ as common factors pesaran2006estimation,greenaway2012asymptotic. Translated to potential outcomes, panel data models with interactive fixed effects assume that $y_{it}(d)=\beta_d d+\beta_x^\top \tilde x_{it}+\lambda_i^\top f_t+\epsilon_{it}.$ Comparing with our model, we see that in our case the effect of $d$ directly interacts with the factors. Thanks to that, the treatment effects can be heterogenous across time and units in our framework. Instead, in panel data models with interactive fixed effects, although it is possible to let $\beta$ depend on $i$ pesaran2006estimation, one cannot have treatment effects heterogenous across both dimensions. Moreover, since this approach assumes that the regressors follow a factor model, we argue that it is more natural to also assume that the potential outcomes also load on the same factors and that $d$ affects the loadings as we do rather than assuming that the effect of $d$ on the potential outcomes is just linear. Note also that if one factor is a constant, our model directly generalizes the panel data models with interactive fixed effects of pesaran2006estimation,greenaway2012asymptotic.

Next, as mentioned earlier, our approach allows us to estimate treatment effects with high-dimensional time series. A leading approach to estimate treatment effects in time series is local projections jorda2023local. Local projection is a linear regression of the outcome on the treatment and some controls and this approach has been extended to high dimensional data babii2022high,adamek2022local. The standard local projection methodology does not allow for time-varying treatment effects. We improve on local projections in this dimension, at the cost of having to assume that the data follows an approximate factor model.

Finally, this paper is related to the literature on characteristics-based factor model connor2007semiparametric,connor2012efficient,fan2016projected. This literature considers a large N,T panel $\{\tilde y_{it}\}_{i,t}$ and assumes that $\tilde y_{it}$ follows a factor model. The loadings are modelled as nonparametric functions of some time-invariant characteristics $\tilde{x}_i$. The loadings and factors can they be estimated by either a semiparametric least squares approach or a projection strategy. The present paper differs from this strand of research in at least two respects. First, we consider a causal inference problem with counterfactuals. Second, in our case the loadings are functions of a time-varying variable $d_{it}$. Since the approach of the aforementioned papers crucially relies on the fact that the loadings are functions of time-invariant characterustics, we cannot build up on their estimation procedure. For this reason, we developed a new methodology, where the factors are first estimated from the panel $\{x_{\ell t}\}_{\ell,t}$. We note that fan2016projected provide a theory where the loadings are estimated nonparametrically using series. This is in the spirit of our approach, but to simplify the exposition of the current paper, we do not let the number of series terms go to $\infty$ in theory. \\

Outline. The paper is organized as follows. Section (ref) introduces the estimation procedure. Then, we expose the asymptotic theory in Section (ref). Extensions are discussed in Section (ref). Section (ref) contains simulations. The methodology is illustrated through an empirical application in Section (ref). Finally, we conclude in Section (ref).\\

Notation. For an integer $n$, we let $[n]=\{1,\dots,n\}$. For a $n_1\times n_2$ matrix $A$, its $k^{th}$ singular value is $\sigma_k(A)$ and $A=\sum_{k=1}^{\min(n_1,n_2)}\sigma_k(A)u_k(A)v_k(A)^\top$ is the singular value decomposition of $A$, where $\left\{u_k\left(A\right)\right\}_{k=1}^{\min(n_1,n_2)}$ is a family of orthonormal vectors of $\mathbb{R}^{n_1}$ and $\left\{v_k\left(A\right)\right\}_{k=1}^{\min(n_1,n_2)}$ is a family of orthonormal vectors of $\mathbb{R}^{n_2}$. The operator norm is $\|A\|_{\text{op}}=\sigma_1(A)$ and the Euclidian norm is $\|A\|_2=\sqrt{\sum_{i=1}^{n_1}\sum_{j=1}^{n_2}A_{ij}^2}$. For an integer $n$, $0_n$ and $1_n$ are, respectively, the $n\times 1$ vectors with all coefficients equal to 0 and 1. Moreover, $I_n$ is the identity matrix of size $n$. For integers $n_1$ and $n_2$, $0_{n_1,n_2}$ and $1_{n_1,n_2}$ are, respectively, the $n_1\times n_2$ matrices with all coefficients equal to 0 and 1. For a random variable $s_{it}$, we let ${\mathbb{E}}_i[s_{it}]=\operatorname*{plim}\limits_{N\to\infty} \frac1N \sum_{i=1}^N s_{it}$ and ${\mathbb{E}}_t[s_{it}]=\operatorname*{plim}\limits_{T\to\infty} \frac1T \sum_{t=1}^T s_{it}$, when the latter quantities exist.

Estimation

In this paper, we estimate the factor by principal components analysis (PCA) applied to the $T\times L$ matrix X such that $X_{t\ell } = x_{\ell t}$. Note that (ref) becomes $$X=F\Lambda^\top+E$$ in matrix form, where $F=(f_1,\dots,f_T)^\top$ is $T\times R$, $\Lambda=(\lambda_1,\dots,\Lambda_L)^\top$ is $L\times R$ and $E$ is the $T\times L$ matrix such that $E_{t\ell}=e_{\ell t}$. Formally, we first obtain an estimate $\widehat{R}$ of $R$ applying to $X$ one of the estimators of the number of factors available in the literature bai2002determining,onatski2010determining,ahn2013eigenvalue,bai2019rank,fan2022estimating. Then, we let the columns of $\widehat{F}/\sqrt{T}$ be the eigenvectors corresponding to the leading $\widehat{R}$ eigenvalues of $XX^\top$. The estimate $\widehat{f}_t\in{\mathbb{R}}^{\widehat{R}}$ of $f_{t}$ corresponds to the $t^{th}$ colmun of $\widehat{F}^\top$. We note that an alternative estimation procedure for the factors would be to use cross-section averages as in pesaran2006estimation.

Once, we obtained an estimate $\widehat{f}_t$ of the factors, the idea is to leverage the fact that $y_{it}=\lambda_{*i}(d_{it})^\top f_t+\tilde u_{it}$ to estimate $\lambda_{*i}$. To simplify estimation, we assume that

where $\varphi_j:{\mathbb{R}}\mapsto {\mathbb{R}}, j\in[J]$ are $J-1$ known non constant functions of $d$ and $\beta_{irj}\in{\mathbb{R}}$ are unknown and need to be estimated. For instance, one can pick $J=2$ and $\varphi_1(d)= d,\varphi_2(d)=d^2$, which yields $\lambda_{*ir}(d)=\beta_{i1r}+\beta_{i2r}d+\beta_{i2r}d^2$, allowing for nonlinearity of the effect of $d$ on $\lambda_{*ir}(\cdot)$. We note that $\lambda_{*ir}(\cdot)$ is identified without imposing (ref) and the latter is only useful for estimation. In the spirit of the literature on series estimation, it would be possible to let $J$ grow with the sample size and to assume instead that $\lambda_{i*}(d)$ is only approximated by $ \sum_{j=1}^J \beta_{jr}\varphi_{j}(d)$. We do not do so in the present paper in order to simplify the exposition and the methodology.

We obtain $$y_{it}= \sum_{r=1}^R \beta_{i0r} f_{tr} +\sum_{j=1}^J\sum_{r=1}^R \beta_{ijr}\varphi_{j}(d_{it}) f_t+ \tilde u_{it},$$ where $\tilde u_{it}=\tilde u_{it}(d_{it})$. Thanks to this, we could estimate $\beta_i=(\beta_{i11},\dots,\beta_{iR1},\dots,\beta_{i1J},\dots,\beta_{iRJ})^\top$ by regressing linearly $y_{it}$ on $f_t$ and $\varphi_j(d_{it})f_t,\ j\in[J]$. Such an approach would only be valid under the exogeneity assumption ${\mathbb{E}}_t[f_t\tilde u_{it}]= {\mathbb{E}}_t[\varphi_{j}(d_{it})f_{t}\tilde u_{it}]=0,\ j\in[J]$. This assumption may not always be credible. To weaken it, we include observed controls $c_{it}\in {\mathbb{R}}^{d_{c}}$, such that $\tilde{u}_{it} =\alpha_i^\top c_{it} +u_{it}$, with ${\mathbb{E}}_t[f_t u_{it}]= {\mathbb{E}}_t[\varphi_{j}(d_{it})f_{t} u_{it}]={\mathbb{E}}_t[c_{it}u_{it}]=0,\ j\in[J]$ and $\alpha_i\in{\mathbb{R}}^{d_c}$. The controls $c_{it}$ can belong to $X$ as is standard in panel data models with interactive fixed effects where the variables used to estimate the factors also act as controls pesaran2006estimation,greenaway2012asymptotic. The second step regression model can then be written

where $\gamma_i=(\beta_{i}^\top,\alpha_i^\top)^\top$ and $$w_{it}=(1_{R}^\top, \varphi_1(d_{it})f_t^\top, \dots,\varphi_J(d_{it})f_t^\top,c_{it}^\top)$$ is $((J+1)R+d_c)\times 1$. A natural estimator $\widehat{\gamma}_i$ of $\gamma_i$ is then the linear regression of $y_{it}$ on

a $((J+1)\widehat{R}+d_c)$-dimensional vector. Formally,

Note that, by (ref) and (ref), we have $\Delta_i=\gamma_i^\top{\mathbb{E}}_t[z_{it}]$, where $$z_{it}=(\varphi_1'(d_{it})f_t^\top,\dots,\varphi_J'(d_{it}) f_t^\top,0_{d_c}^\top)^\top$$ is $((J+1)R+d_c)\times 1$ and $\varphi_j'(\cdot)$ is the derivative of $\varphi_j(\cdot)$. By plug-in the final estimator of $\Delta_i$ is

where

is a $(J+1)\widehat{R}+d_c$-dimensional vector. When $N$ is large, we can also estimate $\Delta_t$ by

where $$\widehat{\gamma}=\frac1N \sum_{i=1}^N \widehat{\gamma}_i$$ is an estimator of $\gamma={\mathbb{E}}[\gamma_i].$ To conclude, we estimate $\Delta$ by the sample average of the $\widehat{\Delta}_i$, that is

Algorithm (ref) summarizes the estimation procedure.

Theory

We place ourselves in an asymptotic regime where $L,T\to \infty$. When we derive the properties of $\widehat{\Delta}_t$ and$\widehat{\Delta}$ , we additionally assume that $N\to\infty$. The number of factors $R$ is fixed with $T$. It would be possible to let it grow and/or allow for nonstrong factors, see, for instance, beyhum2019square,beyhum2022factor,freeman2023linear,bai2023approximate. In this section, as is standard in the literature, we suppose that we know $R$, that is $\widehat{R}=R$ almost surely. Our results would nevertheless stay valid as long as ${\mathbb{P}}(\widehat{R}=R)\to 1$.

Assumptions

We make different types of assumptions. First, there are some standard assumptions for factor models on $F,\Lambda$ and $E$. These assumptions correspond to that made in bai2023approximate and are therefore not new. These are stated in Section (ref). Additional assumptions, not directly imposed on the factor model consisting in $F,\Lambda$ and $E$ are then outlined in Section (ref)

Standard assumptions on $\bf{F,\Lambda,E}$

We let $e_t=(e_{1t},\dots, e_{Lt})^\top$ and $e_{\ell}=(e_{\ell 1},\dots, e_{\ell T})^\top$

We refer the reader to bai2023approximate for a discussion of Assumptions (ref) to (ref).

New assumptions

Let $h_{it}=(u_{it}w_{it}^\top, v_{it}^\top)^\top$, with $v_{it}= z_{it}-{\mathbb{E}}_t[z_{it}]$. The next assumption allows us to show asymptotic normality of $\Delta_i$.

Assumption (ref) (iii) contains mild regularity conditions. The rate condition in Assumption (ref) (iv) is the same as the one used in the literature on factor-augmented regression to show that the error made in estimating the factors is asymptotically negligible bai2006confidence. Assumption (ref) (v) states that certain law of large numbers and central limit theorems hold. Let $$q_{\ell t}={\mathbb{E}}_i[\gamma_i^\top b_{i\ell t}] ,$$ where $$ b_{i\ell t}= ( 0_R^\top, \varphi_{1}'(d_{it}) g_{\ell t}^\top ,\dots, \varphi_{J}'(d_{it}) g_{\ell t}^\top )\ \text{with } g_{\ell t}=\Sigma_\Lambda^{-1} \lambda_\ell e_{\ell t}.$$ We leverage the next assumption only to derive the asymptotic distribution of $\widehat{\Delta}_t$.

Assumption (ref) (ii) can only be satisfied when both $T$ and $L$ go to infinity. If $X$ consists of $K$ variables for each unit, we have $L=NK$ and the condition simplifies to $\sqrt{N}/T\to 0$. Let $$m_t=(0_R^\top,{\mathbb{E}}_i[\varphi_{1}'(d_{it}) ]f_t^\top-{\mathbb{E}}[ \varphi_1'(d_{it})f_t]^\top,\dots,{\mathbb{E}}_i[\varphi_{J}'(d_{it}) ]f_t^\top-{\mathbb{E}}[ \varphi_J'(d_{it})f_t]^\top,0_{d_c}^\top)^\top.$$ The next assumption is only used to establish asymptotic normality of $\widehat{\Delta}$.

Asymptotic distributions

To state the asymptotic distribution of our estimators, we need more notation. Let $$D=\text{diag}((\sqrt{\sigma_1(\Sigma_\Lambda\Sigma_F)},\dots,\sqrt{\sigma_R(\Sigma_\Lambda\Sigma_F)})^\top)$$ be the diagonal matrix consisting in the ordered square-root of the eigenvalues of $\Sigma_\Lambda\Sigma_F$. Moreover, we introduce $Q= D \Upsilon^\top \Sigma_{\Lambda}^{-1/2}$, where $\Upsilon=(u_1(\Sigma_F^{1/2} \Sigma_\Lambda \Sigma_F^{1/2}),\dots, u_R(\Sigma_F^{1/2} \Sigma_\Lambda \Sigma_F^{1/2}))$ consists in the eigenvectors of the matrix $\Sigma_F^{1/2} \Sigma_\Lambda \Sigma_F^{1/2}$. We also let $\Psi$ be the $((J+1)R+d_c)\times ((J+1)R+d_c)$ matrix defined by $$\Psi =\left(

\right).$$

Next, we define the $R\times R$ matrix $\widehat{D}$ whose $r^{th}$ diagonal is the $r^{th}$ singular value of $X/\sqrt{TL}$ and introduce $\widehat{H}=\left(\frac{\Lambda^\top\Lambda}{N} \right)\left(\frac{F^\top\widehat{F}}{T} \right)\widehat{D}^{-2}.$ We also denote by $\widehat{\Psi}$ the $((J+1)R+d_c)\times ((J+1)R+d_c)$ matrix defined as $\Psi$ but replacing the $Q^{-1}$ matrices by $\widehat{H}$. Finally, define the vector $\omega_i= \left( {\mathbb{E}}_t[z_{it}]^\top\Psi^\top(\Psi^\top \Sigma_{w_iw_i}\Psi)^\top, (\Psi^{-1}\gamma_i)^\top \right)^\top.$

We show in the Appendix that, under our assumptions, $\widehat{\Psi}\xrightarrow[T,L\to\infty]{{\mathbb{P}}}\Psi$ and that $\widehat{w}_{it}$, $\widehat{z}_{it}$ and $\widehat{\gamma}_{i}$ estimate consistently $\widehat{\Psi}^\top w_{it}$, $\widehat{\Psi}^\top z_{it}$ and $\widehat{\Psi}^{-1} \gamma_i$, respectively. This allows to derive the asymptotic distributions of our estimators.

The following theorems give the asymptotic distributions of our estimators. First, considering an asymptotic regime with $L,T\to \infty$, we have the following result.

Next, we consider an asymptotic regime with $N,T,L$ all jointly going to infinity.

In Theorem (ref), we assume that $N/T\to c<\infty$ just to have a well defined limiting variance. Note that $c$ can be equal to $0$, so that $N$ can be negligible with respect to $T$.

Estimation of the variance

Let us now discuss estimation of $\sigma_i,\sigma_t$ and $\sigma$. This is crucial to use the theorems for inference.\\

Estimation of $\bf{\sigma_i}$. We can estimate $\omega_i$ by $$\widehat{\omega}_i=\left( \left(\frac1T\sum_{t=1}^T \widehat{z}_{it}\right)^\top\left(\frac1T\sum_{t=1}^T \widehat{w}_{it}\widehat{w}_{it}^\top\right)^\top, \widehat{\gamma}_i^\top \right)^\top.$$ It remains to estimate $ \Sigma_{h_i}.$ First, we can define $\widehat{h}_{it}=(\widehat{u}_{it}\widehat{w}_{it}^\top, \widehat{v}_{it}^\top)^\top,$ where $\widehat{u}_{it}=y_{it}-\widehat{\gamma}_i^\top w_{it}$. There are different potential estimators of $\Sigma_{h_ih_i}$. If the $h_{it}$ are not serially correlated, one can use the “standard” heteroscedasticity consistent (HC) estimator: $$ \widehat{\Sigma}^{HC}_{h_ih_i}=\frac{1}{T} \sum_{t=1}^T \widehat{h}_{it} \widehat{h}_{it}^\top.$$ However, in case of autocorrelation in $h_{it}$ (which may be caused by serial correlation of the factors), one should use an heteroscedasticity and autocorrelation consistent (HAC) estimator à la newey1987a: $$ \widehat{\Sigma}^{HAC}_{h_ih_i}=\sum_{j=1}^Tk\left(\frac{j}{b}\right)\widehat{\Gamma}_{h_ij},$$ where $\widehat{\Gamma}_{h_ij}= \frac{1}{T} \sum_{t=1}^{T-j} \widehat{h}_{it+j} \widehat{h}_{it}^\top$, $k$ is a kernel function and $b>0$ is a bandwidth. The final estimate of $\sigma_i^2$ is given by $$ \widehat{\sigma}_i^2= \widehat{\omega}_i^\top \widehat{\Sigma}_{h_ih_i} \widehat{\omega}_i,$$ where $\widehat{\Sigma}_{h_ih_i} $ is either $\widehat{\Sigma}^{HC}_{h_ih_i}$ or $\widehat{\Sigma}^{HAC}_{h_ih_i} $.\footnote{An alternative approach is to use fixed-b critical values as in lazarus2018har. This would require working out the asymptotic distribution of the test statistic under fixed-b asymptotic, which seems challenging in the present context.}\\

Estimation of $\bf{\sigma_t}$. Our estimator of $\Lambda$ is the $L\times \widehat{R}$ matrix $\widehat{\Lambda} =(\widehat{\lambda}_1,\dots, \widehat{\lambda}_L)^\top=X^\top \widehat{F}/T$. We also let $\widehat{e}_{\ell t}= x_{\ell t}-\widehat{\lambda}_\ell^\top \widehat{f}_t$ be our estimator of $e_{\ell t}$. These are the standard PCA-based estimators. The estimators of $g_{\ell t}$, $b_{i\ell t}$ and $q_{\ell t}$ are then respectively,

Let us denote by $$\widehat{\sigma}_{q_{t}}^2= \frac{1}{L}\sum_{\ell =1}^L \widehat{q}_{\ell t} \widehat{q}_{\ell t}^\top$$ the estimator of $\sigma_{q_t}^2$.\footnote{This estimator will be consistent when $\{q_{\ell t}\}_{\ell}$ are uncorrelated across $\ell$. This will happen if the $\{e_{\ell t}\}_{\ell}$ are uncorrelated across $\ell$ and independent of the loadings. This is an assumption that could be considered natural in factor models, and the results should not be too sensitive to mild violations of this assumption. To avoid such an assumption, one could use a cross-section consistent estimator as in bai2006confidence.} We can estimate $\text{var}_i(\gamma_i^\top z_{it})$ by $$\widehat{\text{var}}_i(\gamma_i^\top z_{it})=\frac1N\sum_{i=1}^N \left(\widehat{\gamma}_i^\top \widehat{z}_{it}- \widehat{\Delta}_t\right)^2.$$ The final estimator of $\sigma_{t}^2$ is then $$\widehat{\sigma}_t^2= \left(\frac{N}{L}\right) \widehat{\sigma}_{q_{t}}^2+ \widehat{\text{var}}_i(\gamma_i^\top z_{it}).$$

Estimation of $\bf{\sigma}$. The variable $ \Psi^\top m_t$ can be estimated by $$\widehat{m}_t= \left(0_R^\top,\widehat{m}_{t1}^\top,\dots,\widehat{m}_{tJ}^\top,0_{d_c}\right)^\top,$$ where $$\widehat{m}_{tj}=\left(\sum_{i=1}^N \varphi_{j}'(d_{it}) \right)\widehat{f}_t^\top-\left(\frac{1}{NT} \sum_{i=1}^N\sum_{t=1}^T \varphi_j'(d_{it})\widehat{f}_t\right)^\top.$$ A natural estimator of $\sigma^2$ is $$\widehat{\sigma}^2=\frac{N}{T}\widehat{\gamma}^\top \widehat{\Sigma}_{mm}\widehat{\gamma}+ \widehat{\text{var}}(\Delta_i),$$ where $\widehat{\gamma}=N^{-1}\sum_{i=1}^N \widehat{\gamma}_i$ estimates $\Psi^{-1}\gamma$, $\widehat{\text{var}}(\Delta_i)=N^{-1}\sum_{i=1}^N (\widehat{\Delta}_i-\widehat{\Delta})^2$ and $ \widehat{\Sigma}_{mm}$ is either $$\widehat{\Sigma}^{HC}_{mm} = \frac1T\sum_{t=1}^T \widehat{m}_t\widehat{m}_t^\top$$ or, in case of serial correlation, $$ \widehat{\Sigma}^{HAC}_{mm}=\sum_{j=1}^Tk\left(\frac{j}{b}\right)\widehat{\Gamma}_{mj},$$ where $\widehat{\Gamma}_{mj}= \frac{1}{T} \sum_{t=1}^{T-j} \widehat{m}_{t+j} \widehat{m}_{t}^\top$, $k$ is a kernel function and $b>0$ is a bandwidth.

Extensions

In this section, we discuss extensions of our approach to other target parameters and IV estimation. This showcases the flexibility of our modelling strategy.\\

Average counterfactuals. Let us now discuss other interesting parameters. First, one may want to know what would be the average outcome if the treatment was fixed at $d\in{\mathbb{R}}$. The following parameters answer this question: $$\varphi_i(d)={\mathbb{E}}_t[y_{it}(d)],\ \varphi_t(d)={\mathbb{E}}_i[y_{it}(d)],\ \varphi(d)={\mathbb{E}}[y_{it}(d)].$$ The mappings $\varphi_i(\cdot)$, $\varphi_t(\cdot)$ and $\varphi(\cdot)$ are functional parameters and they can be estimated as follows: $$\widehat{\varphi}_i(d)=\widehat{\gamma}_i^\top\left(\frac1T\sum_{t=1}^T \widehat{w}_{it}(d)\right),\,\ \widehat{\varphi}_t(d)=\widehat{\gamma}^\top\left( \frac1N\sum_{i=1}^N \widehat{w}_{it}(d) \right),\ \varphi(d)=\frac1N\sum_{i=1}^N \widehat{\varphi}_i(d),$$ where $$ \widehat{w}_{it}(d)=(\widehat{f}_{t}^\top,\varphi_1(d)\widehat{f}_t^\top,\dots,\varphi_J(d) \widehat{f}_t^\top, 0_{d_c}^\top)^\top.$$ The mappings $\varphi_i(\cdot)$, $\varphi_t(\cdot)$ and $\varphi(\cdot)$ can also be differentiated to obtain average marginal effects at fixed d and also integrated (possibly after differentiation) with respect to some relevant distibution of $d$.\\

IV estimation. Suppose now that $d_{it}$ is endogenous even after including the controls, that is ${\mathbb{E}}_t[\varphi_{j}(d_{it})f_{t} u_{it}]\ne 0$ for some $j\in[J]$. A remedy to this problem is to rely on an IV $s_{it}\in {\mathbb{R}}$. Let us define $\widehat{r}_{it} = (\widehat{f}_{t}^\top,\varphi_1(s_{it})\widehat{f}_t^\top,\dots,\varphi_J(s_{it}) \widehat{f}_t^\top ,c_{it}^\top)^\top.$ To estimate $\gamma_i$ one can then use the IV estimator $$\widehat{\gamma}_i^{IV}=\left(\sum_{t=1}^T \widehat{r}_{it}\widehat{w}_{it}^\top \right)^{-1}\sum_{t=1}^T\widehat{r}_{it} y_{it}.$$ The estimators of the other parameters of interest can then be obtained by plug-in. This approach relies on the exogeneity assumption: ${\mathbb{E}}_t[f_t u_{it}]= {\mathbb{E}}_t[\varphi_{j}(s_{it})f_{t} u_{it}]={\mathbb{E}}_t[c_{it}u_{it}]=0$.

Simulations

In this section, we provide a Monte Carlo study which sheds light on the finite sample performance of our proposed inference procedures. We start by analyzing $\widehat{\Delta}_i$ in Section (ref) before studying $\widehat{\Delta}_t$ and $\widehat{\Delta}$ in Section (ref).

Single unit

To analyze $\widehat{\Delta}_i$, we restrict ourselves to a sample with a single unit $i=1$, that is $N=1$. We consider $T\in\{50,100, 200\}$ and $L\in\{50,100,200\}$. We generate samples with $R=2$ factors. The loadings are i.i.d. such that $\lambda_{\ell r}\sim\mathcal{U}[-1,1], \ell\in[L], r\in[R]$. The factors are generated as $f_{tr}=0.5+\rho_f^r (f_{t-1r}-0.5)+ \tilde f_{tr}$ for $t=2,\dots, T$ and $r=1,2$, where $\tilde f_{tr}$ are i.i.d. $\mathcal{N}\left(0,I_R\left(1-\rho_f^{2r}\right)\right)$. The stationary distribution of $f_{tr}$ is $\mathcal{N}(0.5,1)$, we initialize $f_{0r}$ as such. The quantity $\rho_{f}$ controls the level of serial correlation and we let $\rho_f\in \{0,0.5\}$. The idiosyncratic components $\{e_{\ell t}\}$ are i.i.d. $\mathcal{N}(0,1)$. We have two controls, $c_{1t} = (x_{1t},x_{2t})$. We let $\lambda_{*1r}(d)=0.5+\sum_{j=1}^J 0.5\times d^{j},$ where $J\in\{1,2\}$. The treatment is such that $d_{1t}=f_{t1}+0.5\times e_{1t}+0.5\times e_{2t}+\epsilon_t$, where $\epsilon_t\sim\mathcal{N}(0,1)$. The outcome is generated as $$y_{1t} =\lambda_{*1r}(d_{1t})^\top f_t-0.5\times c_{1t}- 0.5\times c_{2t} + u_{1t},$$ where $u_{1t}$ is i.i.d. $\mathcal{N}(0,1)$. We estimate $\Delta_i$ for $i=1$, which is equal to $0.5$ when $J=1$ and $2$ when $J=2$. We use the growth ratio estimator of ahn2013eigenvalue to estimate the number of factors. Following Section (ref), we consider three possible estimators of the covariance matrix $\Sigma_{h_ih_i}$. First, there is the standard estimator $\widehat{\Sigma}^{HC}_{h_ih_i}$, second, there is a HAC estimator $ \widehat{\Sigma}^{HAC}_{h_ih_i}$ with quadratic spectral kernel and third a HAC estimator $ \widehat{\Sigma}^{HAC}_{h_ih_i}$ with Parzen kernel. For the HAC estimators, we pick the bandwidth equal to $1.3T^{1/2}$, which corresponds to the recommendation of lazarus2018har.\footnote{The choice in lazarus2018har is for linear regression and fixed-b critical values, that is a different context than ours.} For each estimator of the covariance matrix, we then construct 95% confidence intervals based on the Gaussian approximation of Theorem (ref).

Tables (ref) and (ref) present the results for $J=1$ and $J=2$, respectively. These are averages over 8,000 replications. We report the bias, variance, mean-squared error (MSE) of our estimator along with the average radius of $95\%$ confidence intervals built using the three different estimators of the covariance matrix and the Gaussian approximation of Theorem (ref). We see that our estimator has low bias, variance and MSE, although they worsen when $J$ or the level of autocorrelation in the factors increase. The coverage of the confidence intervals is close to nominal when the factors are not serially correlated $\rho_f=0$. In this case, the differences in radius and coverage between the three types of confidence intervals is minimal. When $\rho_f=0.5$ the coverage worsens for all three types of CI, but the decrease is much steeper for the standard CI not taking into account time series dependence. This suggests than one should indeed use HAC estimators of the covariance matrix in practice, since they greatly improve the results when there is autocorrelation at little cost in the absence of the latter. The CIs based on the Quadratic Spectral and the Parzen Kernel seem to have similar performance. \afterpage{ \thispagestyle{empty}

}

\afterpage{ \thispagestyle{empty}

}

Large panel

Now, we study $\widehat{\Delta}_t$ and $\widehat{\Delta}$. To do so we consider a panel data with $N\in\{50,100,200\}$ units and $T\in\{50,100, 200\}$ dates. We set $L=2N$, which mimicks the case where the panel $X$ consists in $2$ auxiliary variables (corresponding to $\ell=2i-1$ and $\ell=2i$) for each unit $i$ at each date $t$. We generate the factors $f_t$ and idiosyncratic errors $e_{it}$ as in Section (ref). There are two controls $c_{it} = (x_{(2i-1)t},x_{(2i)t})^\top$. We let $\lambda_{*ir}(d)=\beta_{0i}+\sum_{j=1}^J\beta_{ji}\times d^{j},$ where $\beta_{ji}=0.5+\mathcal{U}[-0.5,0.5]$ and $J\in\{1,2\}$. The treatment is such that $d_{it}=f_{t1}+0.5\times e_{(2i-1)t}+0.5\times e_{(2i)t}+\epsilon_{it}$, where $\epsilon_{it}$ is i.i.d. $\mathcal{N}(0,1)$. The outcome is generated as $$y_{it} =\lambda_{*1r}(d_{it})^\top f_t-0.5\times x_{(2i-1)t}- 0.5\times x_{(2i)t} + u_{it},$$ where $u_{it}$ is i.i.d. $\mathcal{N}(0,1)$.

First, we estimate $\Delta_t$ for $t=1$. It is equal to $0.5(f_{11}+f_{12})$ when $J=1$ and $0.5(f_{11}+f_{12})+2\times 0.5(f_{11}^2+f_{11}f_{12})$ when $J=2$. The results are presented in Tables (ref) and (ref) and are, again, averages over 8,000 replications. The variance $\sigma_{q_t}$ is not affected by autocorrelation in the factors and therefore we only report 95% confidence intervals built using $\widehat{\sigma}_{q_t}$ and Gaussian approximation. The results confirm that the performance of the estimator is indeed not affected by autocorrelation. The bias, variance and MSE of the estimator is again low and the coverage of the confidence intervals close to nominal.

Then we estimate $\Delta$. In this design, we have $\Delta=0.5$ when $J=1$ and $\Delta=2$ when $J=2$. The results are displayed in Tables (ref) and (ref) and are, again, averages over 8,000 replications. As in Section (ref), we find that aucocorrelation or increasing $J$ worsen the results but HAC estimators of the covariance matrix can mitigate the decrease in coverage of the CIs.

\afterpage{ \thispagestyle{empty}

}

\afterpage{ \thispagestyle{empty}

}

Empirical application

To illustrate our methodology, we revisit voigtlander2014skill. This paper considers a panel data with $N=313$ sectors of the U.S. economy observed yearly from 1958 to 2005, that is $T=48$ dates. It investigates if increasing skill intensity of the economy is a cause of increasing wage inequality in the United States. To this end, voigtlander2014skill builds an input skill intensity measure $\sigma_{i,t}$ and evaluate its effect on $\ln(w_{L,i,t}/w_{H,i,t})$, the logarithm of the ratio of the average wage of low-skilled workers in sector $i$ at time $t$, $w_{L,i,t}$, over the average wage of high-skilled workers in that sector at the same date. The goal is to assess if skill intensity in a given sector causes wage inequality in that sector.

This data has recently been analyzed by yin2021focused and juodis2022regularization using common correlated effects approaches. We use the dataset of yin2021focused available at \url{https://www.tandfonline.com/doi/abs/10.1080/07350015.2019.1623044?casa_token=Jmzhpt-330cAAAAA:MrhnmhaKnNDF1uWAl6jwO09Nz6S8Mx1oONJ7OHWmdLEagE0TeZ8pL60jSYL97XcXRjGwPaqYpkY}, it corresponds to a balanced version of the original data of voigtlander2014skill, where some sectors with missing data have been deleted.

We let $y_{it}=\ln(w_{L,i,t}/w_{H,i,t})$ and $d_{it}=\ln(\sigma_{i,t})$. For the panel $\{x_{\ell t}\}_{\ell,t}$ we use 7 variables for each sector at each date including the logarithm of the ratio of high skilled workers over low-skilled workers, capital equipment per worker. These corresponds to the control variables used in yin2021focused and juodis2022regularization. This leads us to $L= 7\times 313=2191$ in the panel $\{x_{\ell,t}\}_{\ell,t}$. The control variables $c_{it}$ correspond to these $7$ variables along with a constant.

Table (ref) reports the estimates and confidence intervals of $\widehat{\Delta}$ for $J=1,2,3$. The various confidence intervals are built as in the simulations and the growth ratio estimator finds one factor. The estimates and associated confidence intervals are not very sensitive to the choice of $J$. The results are significant and the coefficient is negative: increasing skill intensity increases wage inequality. Our estimates are much more negative than that of yin2021focused (ranging from $-0.73$ to $-0.59$) and juodis2022regularization (ranging from $-0.99$ to $-0.59$).

Next, we explore the time trend of the time-specific average marginal effect $\widehat{\Delta}_t$. In Figure (ref), we plot $\widehat{\Delta}_t$ along with its (pointwise) confidence intervals estimated with $J=1$. We see that the average marginal effect of skill intensity tends to increase over time, suggesting that the wage premium of skilled workers (for a fixed value of the skill intensity measure) is decreasing with time.

Conclusion

In this paper, we present a novel factor model for potential outcomes, outlining the methodology for estimating various key parameters within the model. Notably, our approach accommodates both temporal and unit-specific variations in treatment effects. We further establish the asymptotic distribution for several pivotal estimators, paving the way for statistical inference. To enhance the simplicity of inference procedures, the development of a bootstrap technique for constructing confidence intervals across a broad spectrum of estimators in our model would be beneficial. This avenue remains open for exploration in future research endeavors. Another interesting extension is the case where the controls are high-dimensional and lasso is used for estimation. fan2022we demonstrates the importance of such high-dimensional controls in counterfactual analysis with factor models and a binary treatment.

\if11 {

Acknowledgments

Jad Beyhum gratefully acknowledges financial support from the Research Fund KU Leuven through the grant STG/23/014. He thanks Andrii Babii, Geert Dhaene and Jonas Striaukas for useful comments. } \fi