Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,212 characters · 29 sections · 120 citation commands
An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls
\if11 {
} \fi
\if01 {
} \fi
\spacingset{1.5} \pagenumbering{arabic}
We consider the problem of making inferences on the causal effect of a policy intervention in an aggregate time series setup with a single treated unit. The treated unit is observed for $T_0$ periods before and $T_\ast$ periods after the intervention occurs. Often, there is additional information in the form of possibly very many untreated units, which can serve as controls. Such settings are ubiquitous in applied research, and there are many different approaches for estimating the causal effect of the policy. Examples include difference-in-differences methods, synthetic control (SC) approaches, factor, matrix completion, and interactive fixed effects (FE) models, and times series models.\footnote{We refer to DI16, gobillon2016regional, and abadie2019jel for excellent comparative overviews and reviews.} We refer to these methods as counterfactual and synthetic control (CSC) methods.
This paper provides generic and robust procedures for making inferences on policy effects estimated by CSC methods. We propose a general counterfactual modeling framework that nests and generalizes many traditional and new methods for counterfactual analysis. We focus on methods that are able to generate mean-unbiased proxies, $P_t^N$, for the counterfactual outcomes of the treated unit in the absence of the policy intervention, $Y_{1t}^N$: \[ Y_{1t}^N=P_t^N+u_t, \quad E\left(u_t \right)=0, \quad t=1,\dots,T_0+T_\ast. \] The policy effect in period $t$ is $\theta_t=Y_{1t}^I-Y_{1t}^N$, where $Y_{1t}^I$ is the counterfactual outcome of the treated unit with the policy intervention. We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, $\theta=\{\theta_t\}_{t=T_0+1}^{T_0+T_\ast}$. Specifically, we postulate a trajectory $\theta^0= \{ \theta^0_t\}_{t=T_0+1}^{T_0+T_\ast}$ and test the sharp null hypothesis that $\theta=\theta^0$. We also consider the problem of testing hypotheses about per-period effects $\theta_t$ and propose a simple algorithm for constructing pointwise confidence intervals via test inversion.
We recast the inference problem as a (counterfactual) prediction and a structural breaks testing problem. This allows us to build on the literature on conformal prediction vovk2005algorithmic and end-of-sample stability testing dufour1994generalized,andrews2003end to construct inference procedures that are provably robust against misspecification and accommodate many classical and modern high-dimensional methods for estimating $P_t^N$.
The basic idea of our testing procedures is as follows. Suppose that there is only one post-treatment period and that $P_t^N$ is known. Under the sharp null that $\theta_{T_0+1}=\theta_{T_0+1}^0$, we can compute $Y_{1t}^N$ and $u_t=Y_{1t}^N-P_t^N$ for all time periods. If the stochastic shock process $\{u_t\}$ is stationary and weakly dependent, and its distribution is invariant under the intervention, the distribution of the error in the post-treatment period, $u_{T_0+1}$, should be the same as the distribution of the errors in the pre-treatment period, $\{u_t\}_{t=1}^{T_0}$. We operationalize this idea by proposing inference methods in which $p$-values are obtained by permuting blocks of estimated residuals across the time series dimension.
The proposed methods are valid under two different sets of conditions:
The main theoretical results in this paper are finite sample (non-asymptotic) bounds on the size accuracy of our methods; these bounds imply that our methods are exact as $T_0\rightarrow \infty$. Unlike traditional asymptotic results, which are only informative when the sample size is large enough, our non-asymptotic bounds show how different factors affect the finite sample performance. This feature is relevant in CSC applications where sample sizes are often small.
A key feature of our conformal inference methods is that $P_t^N$ is estimated under the null hypothesis based on data from all $T_0+T_\ast$ periods. Estimation under the null guarantees the exact finite sample validity of our procedures if the data are iid or exchangeable. Even when exchangeability fails, imposing the null for estimation is essential for a good performance in CSC applications where $T_0$ is often rather small. Figure (ref) plots the empirical rejection probabilities for testing the null that $\theta_{T_0+1}=0$ when $T_0=19$, $J=50$ (as in our empirical application), $\{u_t\}$ is an AR(1) process, and $P_t^N$ is estimated using SC. The size properties of our method are excellent. By contrast, estimating $P_t^N$ based on the $T_0$ pre-treatment periods without imposing the null yields substantial size distortions. Figure (ref) suggests that imposing the null continues to improve size accuracy even when exchangeability fails and that these improvements can be substantial in small samples.
We make two additional contributions that may be of independent interest. First, we introduce the $\ell_1$-constrained least squares estimator or constrained Lasso raskutti2011minimax as an essentially tuning-free alternative to existing penalized regression estimators and study its theoretical properties. Constrained Lasso nests SC and difference-in-differences, providing a unifying approach for the regression-based estimation of the mean proxies $P_t^N$. Second, we obtain theoretical consistency results for SC estimators in settings with potentially very many control units.
We develop three extensions of our main results. First, we show that our method can be modified to test hypotheses about average effects over time. Second, we extend our method to settings with multiple treated units. Third, we propose easy-to-implement placebo tests for assessing the credibility of inferences based on our method.
Monte Carlo simulations suggest that our procedures exhibit excellent size properties and are robust to misspecification. We find that imposing additional constraints (e.g., using SC instead of the more general constrained Lasso) does not improve power when these additional restrictions are correct but can cause power losses when they are not.
Finally, we re-analyze the causal effect of decriminalizing indoor prostitution on sexually transmitted infections. Following cunningham2018decriminalizing, we exploit the unanticipated decriminalization of indoor prostitution in Rhode Island in 2003. We find that decriminalizing indoor prostitution significantly decreased the incidence of female gonorrhea.
\paragraph{Related Literature.} We contribute to the literature on inference procedures for CSC methods with few treated units. A popular method is the finite population permutation approach of abadie10sc, see also firpo18synthetic and abadie2019jel. This approach permutes the policy assignment and relies on permutation distributions for inference. It corresponds to conventional randomization inference fisher1935design under random assignment of the policy abadie10sc,abadie2019jel. However, random assignment is not plausible in typical CSC applications, and assignment mechanisms are difficult to model and estimate when there are only few treated units abadie2019jel. shaikh2019randomization propose randomization tests for settings with staggered treatment adoption, which encompass the approach of abadie10sc. The main assumption of shaikh2019randomization's approach is that policy adoption follows a Cox proportional hazards model. We do not model the assignment mechanism. Instead, we exploit stationarity and weak dependence of the errors across time in a repeated sampling framework. One advantage of exploiting the time series dimension is that we only require a suitable model for the potential outcome of the treated unit. By contrast, cross-sectional approaches often require estimating models for all units. That is, we only require a good “local” instead of a good “global” fit, which reduces the risk of model misspecification. On the other hand, our approach requires a large number of pre-treatment periods and relies on invariance of the error distribution under the intervention.
There is also an active literature on asymptotic inference methods for CSC models. Several papers focus on testing hypotheses about average or expected effects over time, requiring $T_0$ and $T_\ast$ to be large. li2017estimation, carvalho2018arco, chernozhukov2019ttest, and li2020statistical introduce inference methods based on penalized and constrained regression methods. arkhangelsky2018synthetic propose inference methods for a version of SC with time and unit weights, which admits a weighted regression formulation. Asymptotic inference methods based on factor and interactive FE models are proposed by hsiao2012panel, gobillon2016regional, chan2016policy, li2017estimation, xu2017generalized, and li2018inference. Here we focus on sharp null hypotheses and permutation distributions and provide non-asymptotic performance guarantees. Our approach is generic and valid with many different methods, including constrained regression, factor models, and interactive FE estimators. conley2011inference propose inference methods for difference-in-difference settings with few treated units. They exploit the cross-sectional dimension, relying on weak dependence and stationarity of the error terms across units, which may be hard to justify in typical CSC settings. By contrast, our procedures rely on stationarity and weak dependence of the errors over time. On the other hand, exploiting the time series dimension, our approach requires $T_0$ to be large, whereas conley2011inference allow $T_0$ to be fixed. In related work, cattaneo2021prediction provide prediction intervals for per-period effects estimated by SC methods. Their key observation is that there is both randomness from estimating the SC weights and from the prediction error. They propose a sampling-based inference method based on non-asymptotic probability bounds that accounts for both types of randomness and is valid with stationary and non-stationary data.
We recast the causal inference problem as a (counterfactual) prediction problem and build on the literature on conformal prediction vovk2005algorithmic,vovk2009online,lei2013distribution,lei2014distribution,lei2017distributionfree and on the literature on permutation tests romano1990behavior,lehmann2005testing, which was started by fisher1935design in the context of randomization; see rubin1984bayesianly for a Bayesian justification. On a more general level, our approach is also connected to transformation-based approaches to model-free prediction politis2015modelfree. Let us discuss in more detail the relationship to chernozhukov2018exact, who extend classical conformal prediction to time series settings. Besides a different focus (prediction intervals for future outcome values vs.\ inference on policy effects), there are several important differences. First, we rely on permuting residuals, whereas CWZ18 permute blocks of data. Second, we theoretically analyze different types of permutations. In particular, we study the set of all permutations, which yields precise $p$-values in small samples. This set of permutations cannot be used in the framework of CWZ18 unless the data are iid. Third, we allow for non-stationary data, whereas the prediction methods in CWZ18 are strictly limited to stationary data. Forth, CWZ18 rely on abstract high-level conditions on the test statistics and do not provide any primitive conditions. By contrast, we develop transparent sufficient conditions that can be verified for many traditional and modern CSC methods, and we provide explicit primitive conditions for a large selection of popular approaches. Finally, we establish the validity of our methods with time series data under misspecification and stability, whereas the theoretical results for weakly dependent data in CWZ18 require correct specification and consistency.
Finally, we show that the problem of making inferences on policy effects can be recast as a structural breaks testing problem with a known break date. Therefore, we build on and contribute to the literature on testing for structural breaks and, in particular, to the literature on structural breaks testing using permutation approaches antoch2001permutation,zeileis2013toolbox. Besides a different focus (inference on policy effects vs.\ testing for structural breaks), our paper differs from the existing literature in that we specifically focus on testing at the end of the sample, allow for a very general class of estimators, including modern high-dimensional methods, and provide non-asymptotic performance guarantees under correct specification and misspecification. Our paper is also related to tests for structural breaks at the end of the sample dufour1994generalized,andrews2003end.\footnote{hahn2017synthetic informally suggest applying a variant of andrews2003end's end-of-sample stability test in the context of SC, and ferman2019inference use a version of this test in the context of difference-in-differences approaches with few treated groups.} Let us discuss the differences to andrews2003end's end-of-sample instability test based on subsampling in more detail. First, we focus on causal inference, whereas andrews2003end is concerned with structural breaks testing. Second, our procedures are exactly valid under exchangeability, and we obtain finite sample bounds under weak conditions on the estimators, while the theoretical properties of andrews2003end's test rely on asymptotic analyses. Third, our methods are valid under misspecification, whereas andrews2003end assumes correct specification. Forth, our results under correct specification only require stationarity and weak dependence of $\{u_t\}$, while andrews2003end's test assumes stationarity of the data.\footnote{ andrews2003end briefly comments on page 1681 (comment 4) that his test can be shown to be asymptotically valid under stationary errors but does not provide a formal result.} Finally, our procedures work in conjunction with many modern high-dimensional estimators, whereas andrews2003end focuses on low-dimensional GMM models.
\paragraph{Notation.} For $q\geq 1$, the $\ell_q$-norm of a vector is denoted by $\| \cdot\|_{q}$. We use $\|\cdot\|_0$ to denote the number of nonzero entries of a vector; $\|\cdot\|_{\infty}$ is used to denote the maximal absolute value of entries of a vector. We use the notation $a\lesssim b$ to denote $a \leq cb$ for some constant $c > 0$ that does not depend on the sample size. We use the notation $a \asymp b $ to denote $a\lesssim b$ and $b\lesssim a$. For a set $A$, $|A|$ denotes the cardinality of $A$. For any $a\in\mathbb{R} $, we define $\lfloor a \rfloor =\max\{z\in\mathbb{Z}:z\leq a\} $ and $\lceil a \rceil =\min\{z\in\mathbb{Z}:z\geq a\} $, where $\mathbb{Z}$ is the set of integers. We use $\mathbb{N}$ to denote the set of natural numbers.
We consider a time series of $T$ outcomes for a treated unit, labeled $j=1$. During the first $T_0$ periods, the unit is not treated by a policy and, during the remaining $T-T_0=T_\ast$ periods, it is treated by the policy. Extensions to more than one treated unit are discussed in the Appendix. Our typical setting is where $T_\ast$ is short compared to $T_0$. There may be other units that are not exposed to the policy, and they will be introduced below. We denote the observed outcome of the treated unit by $Y_{1t}$. We employ the potential (latent) outcomes framework neyman1923application,rubin1974estimating and denote potential outcomes with and without the policy as $Y_{1t}^I$ and $Y_{1t}^N$. The effect of the policy intervention in period $t$ is $\theta_t=Y_{1t}^I-Y_{1t}^N$.
Our conformal inference method will rely on the following counterfactual modeling framework, which nests many traditional and new methods for counterfactual policy analysis; see Sections (ref)--(ref) for examples.
Assumption (ref) introduces the potential outcomes, but also postulates an identifying assumption in the form of the existence of mean-unbiased proxies $P^N_t$ such that $E \left( P_t^N\right) = E \left( Y_{1t}^N \right)$. Assumption (ref) allows $\{P_t^N\}$ to be fixed or random and does not impose any restrictions on the dependence between $\{P_t^N\}$ and $\{u_t\}$. In Sections (ref)--(ref), we will discuss specific panel data and time series models that postulate (and identify) what $P_t^N$ is under a variety of conditions. Additional assumptions on the stochastic shock process $\{u_t\}$ will be introduced later, in essence requiring $\{u_t\}$ to be either iid or, more generally, a stationary and weakly dependent process.
Assumption (ref) also postulates that the stochastic shock sequence is invariant under the intervention. This is the fundamental identifying assumption. It requires that the timing of the policy intervention is independent of factors that change the distribution of $\{u_t\}$.\footnote{In principle, we can relax this assumption by specifying, for example, the scale and quantile shifts in the stochastic shocks that result from the policy, and then working with the resulting model; we leave this extension to future work.} If the policy changes the distribution of $\{u_t\}$, one can either interpret our method as a structural breaks test or view the policy effect as random in which case our method yields valid prediction sets; see the Appendix for details.
Often, there is additional information in the form of untreated units, which can serve as controls. Specifically, suppose that there are $J \geq 1$ control units, indexed by $j=2,\dots,J+1$. We assume that we observe all units for all $T$ periods, although this assumption can be relaxed. Let $Y_{jt}$ denote the observed outcome for these untreated units. This observed outcome is equal to the outcome in the absence of the policy intervention, i.e., $Y_{jt} = Y_{jt}^N$ for $2\le j\le J+1$ and $1\le t\le T$. For each unit, we may also observe a vector of covariates $X_{jt}$. This motivates a variety of strategies for modeling and identifying $P_t^N$ as discussed below.
We are interested in testing hypotheses about the trajectory of policy effects in the post-treatment period, $\theta=\left(\theta_{T_0+1},\dots,\theta_{T}\right)'$. Our main hypothesis of interest is
where $\theta^0=\left(\theta^0_{T_0+1},\dots,\theta^0_{T}\right)'$ is a postulated policy effect trajectory. Hypothesis (ref) is a sharp null hypothesis. It fully determines the value of the counterfactual outcome in the absence of the intervention in the post-treatment period since $Y_{1t}^N=Y^I_{1t}-\theta_t=Y_{1t}-\theta_t$. In the Appendix, we show that our method can also be used to test hypotheses about average effects.
To describe our procedure, we write the data under the null hypothesis as $\mathbf{Z}:=\mathbf{Z}(\theta^0)=(Z_1,\dots,Z_{T})'$, where \[ Z_{t}=
\]
Using one of the methods described below, we will obtain a counterfactual proxy estimate, $\hat{P}^N_t$, based on $\mathbf{Z}$, and compute the residuals $\hat{u}=\left(\hat{u}_1,\dots,\hat{u}_T \right)'$, where $\hat{u}_t=Y_{1t}^N-\hat{P}^N_t$ for $1\le t\le T$. Since $P_t^N $ is computed using $\mathbf{Z}=\mathbf{Z}(\theta^0)$, $P_t^N$ is estimated under the null hypothesis, which is essential for a good small sample performance. In “ideal” settings where the data are iid, imposing the null guarantees the model-free exact finite sample validity of our method; see the Appendix for details. By contrast, when $P_t^N$ is estimated based on the pre-treatment data $\{Z_t\}_{t=1}^{T_0}$ without imposing the null, permuting blocks of residuals does not yield procedures with exact finite sample validity, not even with iid data.
Definition of Test Statistic $\mathbf{S}$. We consider the following test statistic: \[ S (\hat u) = S_q(\hat u) = \left ( \frac{1}{\sqrt{T_{\ast}}} \sum_{t=T_0+1}^T|\hat{u}_t|^q \right)^{1/q}. \]
Note that $S$ is constructed such that high values indicate rejection. Different choices of $q$ lead to power against different alternatives. For instance, if the intervention has a large but only temporary effect (i.e., if $|\theta_t|$ is large for few periods), choosing $q=\infty$ yields high power. On the other hand, if the intervention has a permanent effect (i.e., if $\theta_t$ is non-zero for many post-treatment periods), tests using $S_1$ or $S_2$ exhibit good power properties. In our application, we will be using $S_1$, which behaves well under heavy-tailed data. Throughout the paper, when the nature of the statistic is not essential, we write $S = S_q$.
We use (block) permutations to compute $p$-values. A permutation $\pi$ is a one-to-one mapping $\pi:\{1,\dots,T\}\mapsto\{1,\dots,T\}$. We denote the set of permutations under study as $\Pi$ and assume that $\Pi$ contains the identity map $\mathbb{I}$. We focus on two different sets of permutations: (i) the set of all permutations, which we call iid permutations, $\Pi_{\text{all}}$, and (ii) the set of all (overlapping) moving block permutations, $\Pi_{\to}$.\footnote{We can also consider other types of permutations; for example, the “iid block” permutations. Specifically, let $\{b_1, \dots, b_K\}$ be a partition of $\{1,\dots,T\}$, then we collect all the permutations $\pi$ of these blocks, forming the “iid m-block” permutations $\Pi_{mb}$. In our context, choosing $m=T_*$ is natural, though other choices should work as well, similarly to the choice of block size in the time series bootstrap. We refer to CWZ18 for more results on block permutations.} The elements of $\Pi_{\to}$ are indexed by $j \in \{0,1, \dots, T-1\}$, and the permutation $\pi_{j}$ is defined as \[ \pi_{j}(i)=
\] Figure (ref) presents a graphical illustration of $\Pi_{\text{all}}$ and $\Pi_{\to}$.
The choice of $\Pi$ does not matter for the exact finite sample validity of our procedures if the residuals are exchangeable. However, $\Pi_{\text{all}}$ has more elements than $\Pi_{\to}$, allowing us to compute more precise $p$-values and to test at lower significance levels. For the asymptotic validity under estimator consistency, the choice of $\Pi$ depends on the assumptions that we are willing to impose on the stochastic shock sequence $\{u_t\}$ (cf. Section (ref)).
For each $\pi \in \Pi$, let $\hat{u}_\pi=(\hat{u}_{\pi(1)},\dots,\hat{u}_{\pi(T)})'$ denote the vector of permuted residuals.\footnote{If the estimator of $P_t^N$ is invariant under permutations of the data $\{Z_t\}$ across the time series dimension (which is the case for many estimators in Section (ref)), permuting the residuals $\{\hat{u}_t\}$ is equivalent to permuting the data $\{Z_t\}$.} The permutation $p$-value is defined as follows.
Definition of $p$-Value. The $p$-value is
We are often interested in testing pointwise hypotheses about $\theta_t$, $H_0:\theta_{t}=\theta_{t}^0$, and in constructing pointwise confidence intervals for $\theta_{t}$. Pointwise hypotheses can be tested by defining the data under the null as $\mathbf{Z}=\left(Z_1,\dots,Z_{T_0},Z_{t}\right)'$, provided that $P_t^N$ can be estimated based on $\mathbf{Z}$. Pointwise $(1-\alpha)$ confidence intervals for $\theta_{t}$ can be constructed via test inversion as described in Algorithm (ref).
The availability of control units motivates several strategies for modeling the counterfactual mean proxies $P_t^N$. We estimate $P_t^N$ based on the imputed data under the null hypothesis, $\mathbf{Z}(\theta^0)$, and write $Y_{1t}^N$ instead of $Y^I_{1t}-\theta_t^0$ to alleviate the exposition.
The difference-in-differences method postulates the following model for the counterfactual mean proxy DI16: $ P_t^N=\mu+\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}. $ This model automatically embeds the identifying information. The counterfactual mean proxy can be estimated as $ \hat{P}^N_t=\frac{1}{T}\sum_{s=1}^{T}\left(Y^N_{1s}-\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{js}\right)+\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}. $
The canonical SC method abadie2003economic,abadie10sc,abadie2015comparative postulates the following model:
We need to impose an identification condition that allows us to identify the weights $w$, for example:\footnote{More generally, other exclusion restrictions and identifying assumptions could be used. See also abadie10sc, ferman2019synthetic and ferman2019properties, who study the behavior of SC when the data are generated by a factor model.}
The counterfactual is estimated as $\hat{P}^N_t=\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$. We focus on the following canonical SC estimator for $w$:\footnote{This formulation of canonical SC without covariates is due to DI16, who refer to the estimator (ref) as “constrained regression”. Note that unlike DI16, we estimate $w$ under the null hypothesis based on all the data. We focus on the canonical problem (ref) for concreteness. abadie10sc,abadie2015comparative consider a more general version that also includes covariates into the estimation of the weights. Our inference method also works in conjunction with more recently proposed modified versions of SC, such as the augmented SC estimator of benmichael2018augmented.}
As an alternative, we can consider the more flexible model\footnote{The idea to relax the non-negativity constraint on the weights is not new. It first appeared in hsiao2012panel, who compared their factor model approach to SC, and also in valero2015synthetic, who used the cross-validated Lasso to estimate the weights, and in DI16, who used cross-validated Elastic Net for estimation of weights. They do not establish the formal properties of these estimators. Here we emphasize another version of relaxing SC, model (ref), which leads to constrained Lasso (ref). Constrained Lasso demonstrates an excellent theoretical and practical performance: it is tuning-free, performs very well empirically and in simulations, and we prove that it is consistent for dependent data without any sparsity conditions on the weights and that it satisfies the estimator stability condition required for validity under misspecification. We emphasize that this estimator generally differs from the cross-validated Lasso estimator.}
maintaining the same identifying assumption (SC). The counterfactual is estimated as $\hat{P}^N_t=\hat\mu+\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$ by the $\ell_1$-constrained least squares estimator, or constrained Lasso raskutti2011minimax:
The advantage over other penalized regression methods discussed next is that constrained Lasso is essentially tuning-free, does not rely on any sparsity conditions, and is valid for dependent data under weak assumptions. Moreover, constrained Lasso encompasses both difference-in-differences and canonical SC as special cases (by setting $w=(1/J,\dots,1/J)'$ and $\mu=0,w\ge 0$, respectively) and, thus, provides a unifying approach for the regression-based estimation of $P^N_t$.
Section (ref) provides primitive conditions that guarantee that the SC and the constrained Lasso estimators are valid in our framework in settings with potentially many control units (large $J$). Finally, we note that it is straightforward to incorporate (transformations of) covariates $X_{jt}$ into the estimation problems (ref) and (ref).
Consider a linear model for $P_t^N$: $ P^N_t=\mu+\sum_{j=2}^{J+1}w_jY^N_{jt}. $ We maintain the identifying assumption (SC). The counterfactual is estimated by $\hat{P}^N_t=\hat\mu+\sum_{j=2}^{J+1}\hat{w}_jY^N_{jt}$, where
and $\mathcal{P}(w)$ is a penalty function that penalizes deviations away from zero. If it is desired to penalize deviations away from other focal points $w^0$, for example, $w^0 =(1/J, \dots , 1/J)$ used in the difference-in-differences approach, we may always use instead: $\mathcal{P}(w) \leftarrow \mathcal{P}(w- w^0)$. Note that it is straightforward to incorporate covariates $X_{jt}$ into the estimation problem (ref).
Different variants of $\mathcal{P}(w)$ can be considered. Examples include: Lasso tibshirani96regression, where $\mathcal{P}(w)=\lambda \|w\|_1$ and $\lambda$ is a tuning parameter; Elastic Net Zou2005, where $\mathcal{P}(w)=\lambda \left( (1-\alpha)\|w\|_2^2+\alpha \|w\|_1 \right)$ and $\lambda$ and $\alpha$ are tuning parameters; Lava chernozhukov2017lava, where $\mathcal{P}(w)=\inf_{a + b = w}\lambda \left((1-\alpha) \|a\|_2^2+\alpha \|b\|_1 \right)$ and $\lambda$ and $\alpha$ are tuning parameters.
In the context of CSC methods, Lasso was used by valero2015synthetic, li2017estimation, and carvalho2018arco, while DI16 proposed to use Elastic Net. We will impose only weak requirements on the performance of the estimators (pointwise consistency and consistency in prediction norm), which implies that these estimators are valid in our framework under any set of sufficient conditions that exists in the literature.
Consider the following interactive FE model for treated and untreated units:
where $F_t$ are unobserved factors, $\lambda_j$ are unit-specific factor loadings, and $\beta$ is a vector of common coefficients. Model (ref) nests the classical factor model when $\beta=0$ and also covers the traditional linear FE model, in which $\lambda_j'F_t = \lambda_j+F_t$. Consider the following assumption.
The model leads to the following proxy:
Counterfactual proxies are estimated by $\hat{P}_t^N=\hat\lambda_1'\hat{F}_t + X_{1t}'\hat \beta$, where $\hat{\lambda}_1$ and $\hat{F}_t$, and $\hat \beta$ are obtained using the alternating least squares method applied to the model (ref); see, for example, bai2009panel and hansen2019factor for a version with high-dimensional covariates.
hsiao2012panel appears to the be first work that proposed the use of factor models for predicting the (missing) counterfactual responses specifically in SC settings. gobillon2016regional and xu2017generalized employ bai2009panel's estimator in this setting, albeit provide no formally justified inference methods. Formal inference results for interactive FE and factor models in SC designs are developed in chan2016policy and li2018inference among others.\footnote{Factor models are widely used in macroeconomics for causal inference and prediction; see, for example, stock2016factor and the references therein. In microeconometrics, factor models are used for estimation of treatment/structural effects; see, for example, hansen2019factor who use interactive FE models to estimate the effect of gun prevalence on crime.}
Other recent applications to predicting counterfactual responses include amjad2018robust and athey2018matrix (using, respectively, singular value thresholding and the nuclear norm penalization).\footnote{Note that athey2018matrix's analysis applies to a broader collection of problems with general missing data patterns, nesting SC and difference-in-difference problems as special cases.} Our method delivers a way to perform valid inference for policy effects using any of the factor model estimators used in these proposals applied to the complete data under the null.\footnote{Note that in our case the sharp null allows us to impute the missing counterfactual response and apply any of the factor estimators to estimate the factor model for the entire data, which is then used for conformal inference. Hence our inference approach does not provide inference for the counterfactual prediction methods given in those papers. Indeed, there, the missing data entries are being predicted using factor models, whereas in our case the missing data entries are known under the null, and we use any form of low-rank approximation or interactive FE model to estimate the model for the entire data under the null hypothesis.} We shall be focusing on bai2009panel's alternating least squares estimator\footnote{We choose to focus on PCA/SVD and the alternating least squares estimator for the following reasons: (1) they are by far the most widely used in practice, (2) the alternating least squares estimator is computationally attractive and easily accommodates unbalanced data.} and on matrix completion via nuclear norm penalization when verifying our conditions.
If no control units are available, one can use time series models for the single unit exposed to the intervention. For example, consider the following autoregressive model:\footnote{We can also add a moving average component for the errors, but we do not do so for simplicity.}
In model (ref), the mean unbiased proxy is given by $ P_t^N = \mu + \rho (Y_{1{(t-1)}}^N- \mu). $ Note that the policy effect here is transitory, namely it does not feed-forward itself on the future values of $Y_{1t}^I$ beyond the current values.\footnote{We leave the model with persistent feed-forward effects, $Y_{1t}^I = \rho (Y_{1(t-1)}^I) + \theta_t + u_t $, to future work. } Under the null hypothesis, we can impute the unobserved counterfactual as $Y_{1t}^N = Y_{1t} - \theta_t$ and estimate the model using traditional time series methods, and we can conduct inference by permuting the residuals.
The simplest form of the autoregressive model is the AR($K$) process, where the $\rho(\cdot)$ take the form: $\rho ( \cdot ) = \sum_{k=0}^K \rho_k \mathrm{L}^k (\cdot ),$ where $\mathrm{L}$ is the lag operator. There are many identifying conditions for these models, see, for example, Hamilton1994 or brockwell2013time. More generally, we can use a nonlinear function of lag operators, $ \rho (\cdot ) = m( \cdot, \mathrm{L}^1 (\cdot), \dots, \mathrm{L^k} (\cdot ))$, as, for example, when applying neural networks to time series data chen1999improved,chen2001semiparametric, and we refer to the latter for identifying conditions.
A simple and generic way to combine the insights from the panel data and time series models is as follows. Consider the system of equations:
where $C_t^N$ is a panel model proxy for $Y_{1t}^N$, identified by one of the panel data methods. Note that the model has the autoregressive formulation: $ Y_{1t}^N = C_t^N + \rho (Y_{1(t-1)}^N - C_{t-1}^N) + u_t, $ thereby generalizing the previous model.
Here the mean unbiased proxy for $Y_{1t}^N$ is given by $P_t^N = C_t^N + \rho (\varepsilon_{t-1})$. $P_t^N$ is a better proxy than $C_t^N$ because it provides an additional noise reduction through prediction of the stochastic shock by its lag. The model combines any favorite panel model $C_t^N$ for counterfactuals with a time series model for the stochastic shock model in a nice way: we can identify $C_t^N$ under the null by ignoring the time series structure, and then identify the time series structure of the residuals $Y_{1t}^N - C_t^N$. Estimation can proceed analogously. This approach will often improve the size accuracy of our inferential procedures.
When the data are iid (or exchangeable), our procedure is exactly valid in finite samples as shown in the Appendix. In this section, we establish the validity of our inference methods with time series data. Our results are non-asymptotic in nature and, hence, hold in {\it finite samples}. Finite sample bounds are provided for the size properties of our procedure; these bounds imply that our approach is exact as $T_0\rightarrow \infty$. In Section (ref), we establish the validity of our procedure when the estimator of $P_t^N$ satisfies weak and easy-to-verify small error conditions (pointwise consistency and consistency in the prediction norm). This result accommodates non-stationary data and only requires stationarity and weak dependence of the stochastic shock process $\{u_t\}$. In Section (ref), we consider a setting that accommodates misspecification and inconsistent estimators. We show that if the data are stationary and weakly dependent, our procedure is valid, provided that the estimators are stable.
The main condition underlying the results in this section is the following assumption on the stochastic shock process.
Assumption (ref) allows the data to be non-stationary and exhibit general dependence patterns. Assumption (ref).(ref) of iid shocks is our first sufficient condition. Under this condition, we will be able to use iid permutations, giving us a precise estimate of the $p$-value. The iid assumption can be replaced by Assumption (ref).(ref), which holds for many commonly encountered stochastic processes such as ARMA and GARCH. It can be easily replaced by an even weaker ergodicity condition, as can be inspected in the proofs. Under this assumption, we will have to rely on the moving block permutations.
We also impose the following condition on the estimation error under the null hypothesis. Let $P^N=(P^N_1,\dots,P^N_T)'$ and $\hat P^N=(\hat{P}^N_1,\dots,\hat{P}^N_T)'$.
Assumption (ref) imposes weak and easy-to-verify conditions on the performance of the estimators $\hat P_t^N$ of the counterfactual mean proxies $P_t^N$. These conditions are readily implied by the existing results for many estimators discussed in Section (ref). In Section (ref), we provide explicit primitive conditions and references to primitive conditions implying Assumption (ref).\footnote{While our general results in this section are non-asymptotic, some of the analysis in Section (ref) will not be non-asymptotic in nature.}
The above bound is non-asymptotic, allowing us to claim uniform validity with respect to a rich variety of data generating processes. Using simulations and empirical examples, we verify that our tests have good power and generate meaningful empirical results. There are other considerations that also affect power. For example, the better the model for $P_t^N$, the less variance the stochastic shocks will have, subject to assumed invariance to the policy. The smaller the variance of the shocks, the more powerful the testing procedure will be.
Misspecification is an important practical concern, and consistency of the estimators of the counterfactual mean proxies $P_t^N$ may be questionable in certain settings. The classical analysis of misspecification focuses on convergence to pseudo-true values white1996estimation. If it is possible to show that the estimator of the counterfactual mean proxy, $\hat{P}_t^N$, is consistent for some pseudo-true value $P_t^{N\ast}$ and that $\left\{Y_{1t}^N-P_t^{N\ast}\right\}_{t=1}^T$ is stationary and weakly dependent, the theoretical results in Section (ref) imply the validity of our procedure. Pseudo-true consistency can often be verified for low-dimensional models, but consistency results under misspecification remain elusive in high-dimensional settings. Therefore, we consider a notion of approximate exchangeability, which only requires the estimator to be stable instead of consistent for a pseudo-true value. This stability condition does not require $\hat{P}_t^N$ to be consistent for anything, nor does it rely on correct specification of the counterfactual mean proxies. In the Appendix, we illustrate the difference between consistency and stability based on the analytically tractable example of Ridge regression.
The basic idea underlying the theoretical analysis here is as follows. If the estimators are non-random or independent of the data, then stationarity and weak dependence of the data would mean that $\hat{p}$ based on moving block permutations approximately has a uniform distribution under the null. This result follows from uniform laws of large numbers for dependent data. However, in practice, the estimators are computed using the data and are thus not independent of the data. Our key insight is that stable estimators are approximately independent of individual observations.
We now formalize the notion of stability of an estimator. To emphasize the dependence of $S(\hat{u}) $ on the estimator, with a slight abuse of notation, we write $S(\mathbf{Z},\beta)=\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\beta)$. Let $\{\tilde{Z}_{t}\}_{t=1}^{T}$ be iid from the distribution of $Z_{1}$ and independent of $\mathbf{Z}$. For any $H\subset\{1,\dots,T\}$, let $Z_{t,H}=Z_{t}\mathbf{1}\{t\notin H\}+\tilde{Z}_{t}\mathbf{1}\{t\in H\}$, and $\mathbf{Z}_{H}=\{Z_{t,H}\}_{t=1}^{T}$. Hence, $\mathbf{Z}_{H}$ is a perturbed version of $\mathbf{Z}$ under $H$, i.e., $\mathbf{Z}$ with elements in $H$ replaced by $\{\tilde{Z}_{t}\}_{t\in H}$.
By stability, we mean that the estimator computed using $\mathbf{Z}$ is similar to that computed using $\mathbf{Z}_H$ for $H\in \mathbb{H}$. Let $R\in\mathbb{N}$ and define $m=\left\lfloor T_{0}/R\right\rfloor $. The class $\mathbb{H}=\{\widetilde{H}_{1},\dots,\widetilde{H}_{R}\}$ contains $R$ members with $|\widetilde{H}_j|\leq 3m$ elements. The plan is to require stability under $R\asymp T_0/\log(T_0)$ (so $|\widetilde{H}_j|\asymp \log(T_0)$). Since $\log(T_0)\ll T_0$, swapping out $O(\log(T_0))$ out of $T_0+T_*$ data points should not cause a large change in the estimator for reasonable estimators.
We now give precise definitions of sets in $\mathbb{H}$. For $j\in\{1,\dots,R\}$, let $H_{j}=\{(j-1)m+1,\dots,jm\}$. Since the test statistic depends on $T_*$ data points after obtaining the estimator, defining $\mathbb{H}$ to be $\{H_1,\dots,H_R \}$ is not enough for technical arguments; we need a “wedge” to ensure that these $T_*$ data points do not cause a problem. To do so, we enlarge $H_j$ as follows. Let $k\in\mathbb{N}$ satisfy $T_{*}<k<m$. We let $\widetilde{H}_{j}$ denote the $k$-enlargement of $H_{j}$, i.e., $\widetilde{H}_{j}=\{s: \min_{t\in H_{j}}|s-t|\leq k\}$. Note that $\widetilde{H}_{j}=\{(j-1)m+1-k,\dots,jm+k\}$ for $2\leq j\leq R-1$, $\widetilde{H}_{1}=\{1,\dots,m+k\}$ and $\widetilde{H}_{R}=\{(R-1)m+1-k,\min\{Rm+k,T\}\}$.
Assumption (ref) specifies the estimator stability condition. It strengthens the perturb-one sensitivity of lei2017distributionfree. When the model is misspecified, Assumption (ref) holds whenever the estimator $\hat{\beta}(\mathbf{Z})$ is consistent to a pseudo-true parameter value. However, it is more general in that the estimator $\hat{\beta}(\mathbf{Z})$ need not converge to any non-random quantity as long as it is stable under perturbations in a few observations. This feature is crucial in our setting as it allows us to accommodate high-dimensional CSC methods for many of which consistency results under misspecification are not available. Primitive sufficient conditions for Assumption (ref) are provided in the Appendix.
Let $\Psi(x;\beta)=P(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\beta)\leq x)$. Our strategy is to show that, under the null hypothesis, $\hat{F}(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z})))$ is approximately uniform on $(0,1)$. We exploit the stability condition in Assumption (ref) and show that $\hat{F}(\phi(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z})))$ can be approximated by $\Psi\left(\phi(\bar{Z}_{T_{0}+1},\dots,\bar{Z}_{T_{0}+T_{*}};\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}}));\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}})\right) $, which has the uniform distribution on (0,1). Here $(\bar{Z}_{T_{0}+1},\dots,\bar{Z}_{T_{0}+T_{*}})$ has the same distribution as $(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}})$ and is independent of $\mathbf{Z}_{\widetilde{H}_{R}}$. This essentially confirms the above intuition that for stable estimators, $\hat{\beta}(\mathbf{Z})$ is almost independent of the last few observations $(Z_{T_{0}+1},\dots,Z_{T_{0}+T_{*}})$.
We impose the following regularity conditions on the data.
Stationarity and $\beta$-mixing are commonly imposed conditions on time series data. For a large class of Markov chains, GARCH and various stochastic volatility models, $D_{3}=1$ carrasco2002mixing. Let $(\dot{Z}_{T_{0}+1},\ldots,\dot{Z}_{T_{0}+T_{*}})$ be an independent copy of $ (Z_{T_{0}+1},\ldots,Z_{T_{0}+T_{*}})$ and also independent of $(\mathbf{Z},\{\tilde{Z}\}_{t=1}^{T})$. The bounded derivative of $\Psi\left(x;\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}})\right)$ condition says that the density of $\phi(\dot{Z}_t,\ldots,\dot{Z}_{t+T_{*}-1};\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}}))$ conditional on $\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{j}})$ is bounded by $\xi_{T}$ with high probability. The bounded density condition states that the distribution of the residual does not collapse into a degenerate one or one with point mass. In many cases, $\xi_{T}=O(1)$ for continuous distributions. For example, if $(Y_t,X_t)$ is jointly Gaussian and the variance of $Y_t$ given $X_t$ is bounded below by a constant, then for any $w$ satisfying the SC restrictions, the density of $Y_t-X_t'w$ is bounded by a constant that does not depend on $w$.
The following result states the approximate validity of our testing procedure.
In the theoretical arguments, we actually show a stronger result. The above bound holds for $E|P(\hat{p}\leq\alpha\mid \hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}}))-\alpha|$. Since the stability condition states that $\hat{\beta}(\mathbf{Z}_{\widetilde{H}_{R}})\approx\hat{\beta}(\mathbf{Z})$, this means that $\hat{p} $ conditional on $\hat{\beta}(\mathbf{Z})$ almost has a uniform distribution on $(0,1)$; with iid or exchangeable data, $\hat{p} $ conditional on $\hat{\beta}(\mathbf{Z})$ has an exact uniform distribution. Therefore, we can view Theorem (ref) as a result for approximate exchangeability.
Due to the exponential decay of $\beta_{{\rm mixing}}(\cdot)$, the bound in Theorem (ref) tends to zero if we choose $k$ to be a slowly growing sequence and $T_{0}/R$ to be of the same order. For example, we can choose $k$ and $R$ such that $k\asymp T_{0}/R \asymp \log T_0$. Since $|\widetilde{H}_{j}|=\left\lfloor T_{0}/R\right\rfloor +2k $, Assumption (ref) only requires that the changes to $S(\mathbf{Z}^{\pi},\hat{\beta}(\mathbf{Z}))$ are small if we replace only $\log T_{0}$ observations in computing $\hat{\beta}(\mathbf{Z})$. Under finite dependence, it suffices to choose $k$ and $T_0/R$ to be large enough constants. Note that $R$ is only needed in the theoretical arguments; we do not need to choose $R$ when implementing the proposed procedure.
The theoretical analysis in this section suggests that allowing for both unrestricted patterns of non-stationarity and misspecification is not possible in general. To obtain valid inferences with non-stationary data, one has to either rely on correct specification and consistency or impose assumptions on the particular structure of the non-stationarity, which allow for pre-processing the data to make them stationary.
In this section, we revisit the representative models of counterfactual proxies introduced in Section (ref). Primitive conditions are provided to guarantee that the estimation of the counterfactual mean proxies is accurate enough for the asymptotic validity of the proposed procedure. In particular, these conditions can be used to verify Assumption (ref). The regularity conditions (e.g., bounded moments, weak serial dependence) for different models are stated in the Appendix and are commonly imposed in the literature for these models. The counterfactual mean proxies $P_t^N$ are estimated based on the imputed data under the null, $\mathbf{Z}(\theta^0)$, and we write $Y_{1t}^N$ instead of $Y^I_{1t}-\theta_t^0$ to alleviate the exposition. All the results in this section assume that $T_0\rightarrow \infty$ and $J\rightarrow \infty$ (if $J$ is present in the model).
In Section (ref), we have seen that the counterfactual mean proxies implied by the canonical difference-in-differences model are: $ P_t^N=\mu+J^{-1}\sum_{j=2}^{J+1}Y^N_{jt}. $ We consider the following estimator: $ \hat{P}^N_t=\hat{\mu}+\frac{1}{J}\sum_{j=2}^{J+1}Y_{jt},$ where $ \hat{\mu}=\frac{1}{T}\sum_{t=1}^{T}\left(Y^N_{1t}-\frac{1}{J}\sum_{j=2}^{J+1}Y^N_{jt}\right)=\mu+ \frac{1}{T}\sum_{t=1}^{T} u_t. $ Since $\hat{P}_t^N-P_t^N=\hat \mu -\mu$, Assumption (ref) holds for the simple difference-in-differences model provided that $ T^{-1}\sum_{t=1}^{T} u_t = o_P(1)$, which is true under very weak conditions.
Several models in Section (ref) (including SC and constrained Lasso) imply a structure in which the counterfactual proxy is a linear function of observed outcomes of untreated units.
To provide a unified framework for these models, we use $Y$ to denote a generic vector of outcomes and $X$ to denote the design matrix throughout this section. For example, in Section (ref), we set $Y=Y^N_1$ and $X=(Y^N_2,\ldots,Y^N_{J+1}) $, where $Y^N_j=(Y^N_{j1},\ldots,Y^N_{jT})' \in \mathbb{R}^T$ for $1\leq j\leq J+1 $. These models can be written as
where $u=(u_1,\dots,u_{T})' \in \mathbb{R}^T$. Identification is achieved by requiring that $X$ and $u$ be uncorrelated (cf.\ Condition (SC)).
Under the framework in (ref), different models correspond to different specifications for the weight vector $w$. For the SC model in Section (ref), $w $ is an unknown vector whose elements are nonnegative and sum up to one. More generally, one can simply restrict $w$ to be any vector with bounded $\ell_1 $-norm. This is the constrained Lasso estimator.
Since $P_t^N $ is the $t$-th element of the vector $Xw$, the natural estimator is $\hat{P}_t^N $ being the $t$-th element of $X\hat{w} $, where $\hat{w}$ is an estimator for $w$. The estimation of $w$ depends on the specification. Let $\mathcal{W} $ be the parameter space for $w$. We consider the following version of the original SC estimator
The constrained Lasso estimator is
where $K$ is bounded and $K>0$. In light of the estimator (ref), a natural choice is $K=1$.
In general, we choose the parameter space $\mathcal{W} $ to be an arbitrary subset of an $\ell_1 $-ball with bounded radius. The following result gives very mild conditions under which the constrained least squares estimators are consistent and satisfy Assumption (ref).\footnote{To simplify the exposition, we do not include an intercept in Lemma (ref). Similar arguments could be used to prove an analogous result with an unconstrained intercept.}
Lemma (ref) provides several features that are important for counterfactual inference in our setup. First, we allow $J$ to be large relative to $T$. To be precise, we only require $\log J=o(T^c) $, where $c>0$ is a constant depending only on the $\beta$-mixing coefficients; see the Appendix for details. This is particularly relevant for settings in which the number of (potential) control units and the number of time periods have a similar order of magnitude as in our empirical application in Section (ref). Second, Lemma (ref) does not rely on any sparsity assumptions on $w$, allowing for dense vectors. Third, compared to typical high-dimensional estimators (e.g., Lasso or Dantzig selector), our estimator does rely on tuning parameters that can be difficult to choose in times series settings. Finally, Lemma (ref) provides new theoretical consistency results for the canonical SC estimator in settings with time series data and potentially very many control units.
The models for counterfactual proxies introduced in Section (ref) have factor structures. We provide estimation results for pure factor models (without regressors), factor models with regressors (interactive FE models), and matrix completion models. In this subsection, following standard notation, we let $N=J+1$.
Recall from Section (ref) the standard factor model $ Y_{jt}^{N}=\lambda_{j}'F_{t}+u_{jt}, $ where $F=(F_1,\ldots,F_{T})'\in \mathbb{R}^{T\times k} $ and $\Lambda=(\lambda_{1},\ldots,\lambda_{N})'\in \mathbb{R}^{N\times k}$ represent the $k$-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for $Y_{1t}^N $ is $P_t^N=\lambda_{1}'F_{t} $. We identify $P_t^N $ by imposing the condition that the idiosyncratic terms and the factor structure are uncorrelated (cf.\ Condition (FE)).
We use the standard principal component analysis (PCA) for estimating $P_t^N $.\footnote{Note that PCA amounts to singular value decomposition, which can be computed using polynomial time algorithms, trefethen1997numerical. } Let $Y^{N}\in\mathbb{R}^{T\times N}$ be the matrix whose $(t,j)$ entry is $Y_{jt}^{N}$. We compute $\hat{F}=(\hat{F}_1,\ldots,\hat{F}_T)'\in\mathbb{R}^{T\times k}$ to be the matrix containing the eigenvectors corresponding to the largest $k$ eigenvalues of $Y^{N}(Y^{N})'$ with $\hat{F}'\hat{F}/T=I_{k}$. Let $\hat{\lambda}_{j}'$ denote the $j$-th row of $\hat{\Lambda}=(Y^{N})'\hat{F}/T$. Let $\hat{F}_{t}'$ denote the $t$-th row of $\hat{F}$. Our estimate for $P_{t}^N $ is $\hat{P}_t^N= \hat{\lambda}_{1}'\hat{F}_t $. The following lemma guarantees the validity of this estimator in our context under mild regularity conditions.
The only requirement on the sample size is that both $N$ and $T$ need to be large. Similar to Theorem 3 of bai2003inferential, we do not restrict the relationship between $N$ and $T$. This is flexible enough for a wide range of applications in practice as the number of units is allowed to be much larger than, much smaller than, or similar to the number of time periods.
Now we study the general form of panel models with interactive FEs. Following Section (ref), these models take the form $ Y_{jt}^{N}=\lambda_{j}'F_{t}+X_{jt}'\beta+u_{jt}, $ where $X_{jt}\in\mathbb{R}^{k_{x}}$ are observed covariates and $F=(F_1,\ldots,F_{T})'\in \mathbb{R}^{T\times k} $ and $\Lambda=(\lambda_{1},\ldots,\lambda_{N})'\in \mathbb{R}^{N\times k}$ represent the $k$-dimensional unobserved factors and their loadings, respectively. The counterfactual proxy for $Y_{1t}^N $ is $P_t^N=\lambda_{1}'F_{t}+X_{1t}'\beta $. In this model, we identify the counterfactual proxy through the condition that the idiosyncratic terms are independent of the factor structure and the observed covariates (cf.\ Condition (FE)).
The two most popular estimators are the common correlated effects (CCE) estimator by pesaran2006estimation and the iterative least squares estimator by bai2009panel. We focus on the iterative least squares approach, but analogous results can be established for CCE estimators. The notations for $F_t $, $\lambda_{j} $, $\hat{F}_t $ and $\hat{\lambda}_{j} $ are the same as before. We compute $$ (\hat{F},\hat{\Lambda},\hat{\beta})=\underset{F,\Lambda,\beta}{\arg \min } \sum_{t=1}^T \sum_{j=1}^{N} (Y^N_{jt}-X_{jt}'\beta-F_t'\lambda_{j} )^2\quad \text{ s.t. }\quad F'F/T=I_k \quad \Lambda'\Lambda = \text{Diagonal}_k. $$
The estimate for $P_t^N $ is $\hat{P}_t^N=\hat{\lambda}_{1}'\hat{F}_t+X_{1t}'\hat{\beta} $. The following result states the validity of applying this estimator in conjunction with our inference method.
Under the conditions in Theorem 3 of bai2009panel, $N$ is of the same order as $T$ so that rate is really $T^{-1/2}$; however, the stated bound should hold more generally.
Suppose that
where $M_{jt}$ is the $(j,t)$-element of an unknown matrix $M\in\mathbb{R}^{(J+1)\times T}$ satisfying $\|M\|_{*}\leq K$, where $\|\cdot\|_{*}$ denotes the nuclear norm (the sum of singular values). We observe $Y_{jt}^N$ for $(j,t)\in \{1,\dots,T\}\times \{1,\dots,J+1\} \backslash \{(1,t):T_0+1\leq t\leq T\} $. The identifying condition is that $E(u\mid M)=0$ and that conditional on $M$, $\{u_{j}\}_{j=1}^{J+1}$ is independent across $j$, where $u_{j}=(u_{j1},\dots,u_{jT})'\in\mathbb{R}^{T}$. The counterfactual proxy is $P_{t}^N=M_{1t} $ for $1\leq t\leq T$.
The main challenge is to recover the entire matrix $M$ despite the missing entries $\{Y_{1t}^N:T_0+1\leq t\leq T\} $. The literature on matrix completion considers the model (ref) under the assumption of missingness at random and exploits the assumption that the rank of $M$ is low.\footnote{See, for example, candes2009exact, recht2010guaranteed, candes2011tight, koltchinskii2011nuclear, negahban2011estimation, rohde2011estimation, and chatterjee2015matrix.} Recently, athey2018matrix introduce this method to study treatment effects in panel data models and point out the unobserved counterfactuals correspond to entries that are missing in a very special pattern, rather than at random. Assuming the usual low rank condition on $M$, they employ the nuclear norm penalized estimator and provide bounds on the estimation error in the typical setup of causal panel data models.
We take a different approach here since our main goal is hypothesis testing instead of estimation. The key observation is that under the null hypothesis, there are no missing entries in the data. By imposing the null hypothesis, we replace the missing entries with the hypothesized values and obtain a dataset that contains $\{Y_{jt}^N:1\leq j \leq J+1,\ 1\leq t \leq T \}$. The estimator for $M$ we examine here is closely related to existing nuclear norm regularized estimators and is defined as
where $K>0$ is the bound on the nuclear norm of the true matrix. In principle, it can be a sequence that tends to infinity. When $M$ represents a factor structure with strong factors, $K$ can be shown to grow at the rate $\sqrt{NT} $. Clear guidance on how to choose $K$ is still unavailable, but following athey2018matrix one can use cross-validation.\footnote{The properties of cross-validation remain unknown in these settings.} Alternatively one can use a pilot thresholded SVD estimator to get a sense of what $K$ is, and use a somewhat larger value of $K$. The following result guarantees the validity of this estimator in our context under mild regularity conditions.
The result is notable because no sub-Gaussian assumptions are required. The estimator in ((ref)) does not explicitly require a low-rank condition on $M$. Instead, we impose a growth restriction on $K$. When $M$ is generated by a strong factor structure and the null hypothesis contains full information on the missing entries, we can choose $K\asymp \sqrt{NT}$ and our consistency result holds as long as $N,T\rightarrow \infty $ and $E(|u_{jt}|^{2+c}\mid M) $ is uniformly bounded for some $c>0$. In the case of weak factors, we can choose $K \ll \sqrt{NT} $ and obtain consistency.
As pointed out in Section (ref), time series models can be used to model counterfactual proxies with or without control units. We now discuss low-level conditions under which fitting these models yields estimates good enough for the purpose of our conformal inference approach.
The linear autoregressive model with $K$ lags can be written as $ Y_{1t}^{N}=\rho_{0}+\sum_{j=1}^{K}\rho_{j}Y_{1t-j}^{N}+u_{t}, $ where $\{u_{t}\}_{t=1}^{T}$ is an iid sequence with $E(u_{t})=0$.\footnote{Here the model seems different, but Section (ref)'s model implies this one with $\rho_0 = \mu (1 - \sum_{j=1}^{K}\rho_{j})$.} The counterfactual proxy for $Y_{1t}^N $ is $P_t^N= \rho_{0}+\sum_{j=1}^{K}\rho_{j}Y_{1t-j}^{N}$. We write $P_t^N$ as $P_t^N = y_t'\rho$, where $y_{t}=(1,Y_{1t-1}^{N},Y_{1t-2}^{N},\dots,Y_{1t-K}^{N})'\in\mathbb{R}^{K+1}$ and $\rho=(\rho_0,\dots,\rho_K)'\in \mathbb{R}^{K+1}$. The coefficient vector $\rho$ can be estimated using least squares: $\hat{\rho}=\left(\sum_{t=K+1}^{T}y_{t}y_{t}'\right)^{-1}\left(\sum_{t=K+1}^{T}y_{t}Y_{1t}^{N}\right)$. The estimator for $P_t^N $ is $\hat P_t^N = y_t '\hat \rho$.
As mentioned in Section (ref), we can also apply nonlinear autoregressive models $ Y_{1t}^N=\rho(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N)+u_{t}, $ where $\rho $ is a nonlinear function, in which case the counterfactual proxy is $P_t^N= \rho(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N)$.
Let $\hat{\rho} $ be an estimator for $\rho$ and $\hat{P}_{t}^N= \hat{\rho}(Y_{1t-1}^N,Y_{1t-2}^N,\ldots,Y_{1t-K}^N) $. This estimator can be parametric, semiparametric, or fully nonparametric and is only required to be consistent.
The primitive regularity conditions and the definitions of the neural network estimators possessing these properties can be found, for example, in chen1999improved and chen2001semiparametric.
Here we provide generic conditions for the fused panel/time series models described in Section (ref). In particular, AR models can be used to filter the estimated residuals and obtain near iid errors. In Equation ((ref)) of Section (ref), we introduce an autoregressive structure in the error terms: $ Y_{1t}^N=C_t^N+\varepsilon_t$ and $\varepsilon_t=\rho(\varepsilon_{t-1})+u_t, $ where $C_t^N $ can be specified as a panel data model discussed before. Due to the autoregressive structure in $\varepsilon_{t}$, the counterfactual proxy is $P_t^N=C_t^N+\rho(\varepsilon_{t-1}) $.
We estimate $P_t^N $ via a two-stage procedure. In the first stage, we estimate $C_t^N $ using the techniques we considered before and obtain say $\hat{C}_t^N $. In the second stage, we estimate $\rho(\varepsilon_{t-1}) $ by fitting an autoregressive model to the estimated residuals $\{\hat{\varepsilon}_t\}_{t=1}^T $, where $\hat{\varepsilon}_t=Y_{1t}^N-\hat{C}_t^N $. For simplicity, we consider a linear model in the second stage estimation. Analogous results can be obtained for more general models. To be specific, assume that $\varepsilon_{t}=x_{t}'\rho+u_{t}$, where $x_{t}=(\varepsilon_{t-1},\varepsilon_{t-2},\dots,\varepsilon_{t-K})'\in\mathbb{R}^{K}$ and $\rho=(\rho_{1},\rho_{2},\dots,\rho_{K})'\in\mathbb{R}^{K}$.
Given $\{\hat{\varepsilon}_t \}_{t=1}^T $ from the first-stage estimation, we define $\hat{x}_{t}=(\hat{\varepsilon}_{t-1},\hat{\varepsilon}_{t-2},\dots ,\hat{\varepsilon}_{t-K})'\in\mathbb{R}^{K}$ and $\hat{\rho}=\left(\sum_{t=K+1}^{T}\hat{x}_{t}\hat{x}_{t}'\right)^{-1}\left(\sum_{t=K+1}^{T}\hat{x}_{t}\hat{\varepsilon}_{t}\right)$. To compute the $p$-value, we use $\{\hat{u}_t\}_{t=K+1}^T $ with $\hat{u}_t=\hat{\varepsilon}_t-\hat{x}_t'\hat{\rho} $ in the permutation. By the following result, this procedure is valid under very mild conditions for the first-stage estimation.
Note that the conditions in Lemma (ref) for the autoregressive part are the same as in Lemma (ref). Consistency of $\hat{C}_t^N $ can be verified using existing results, for example, those in Sections (ref)--(ref).
We revisit the analysis in cunningham2018decriminalizing who study the impact of decriminalizing indoor prostitution. They consider the case of Rhode Island, where a judge unanticipatedly decriminalized indoor sex work in July 2003 such that, until the recriminalization in November 2009, Rhode Island had decriminalized indoor and prohibited street prostitution.
We focus on the effect of legalizing indoor prostitution on female gonorrhea incidence. Our outcome of interest is log female gonorrhea incidence per 100,000. We use the data on gonorrhea cases from the Center for Disease Control (CDC)'s Gonorrhea Surveillance Program previously analyzed by cunningham2018decriminalizing; see their Section 3 for a detailed description and descriptive statistics. The female gonorrhea series date back to 1985 such that $T_0=19$ and $T_\ast=6$. Figure (ref) displays the raw data for Rhode Island and the rest of the U.S. states.
We apply three different CSC methods: difference-in-differences, canonical SC, and constrained Lasso with $K=1$. Recall that constrained Lasso nests both difference-in-differences and SC. Following cunningham2018decriminalizing, the set of potential control units includes all other U.S. states and the District of Columbia ($J=50$). We choose $S_1$ as our test statistic and report $p$-values computed based on moving block and iid permutations.\footnote{To keep computation tractable, we randomly sample 10,000 iid permutations with replacement.} All computations were performed in R R2020.
Before turning to the main results, we use the placebo tests proposed in the Appendix to assess the plausibility of the underlying assumptions. Specifically, based on the pre-treatment data, we test $H_0: \theta_{2003-\tau+1}=\dots=\theta_{2003}=0$ for $\tau\in \{1,2,3\}$. Rejections of this null undermine the credibility of the assumptions underlying our procedure and the inferences on policy effects in the post-treatment period. Table (ref) presents the results. Figure (ref) complements the formal tests with plots of the residuals from fitting the three models to the pre-treatment data. The placebo tests and the residual plots provide evidence in favor of the credibility of our inference method in conjunction with SC and, especially, constrained Lasso, but suggest that the difference-in-differences results need to be interpreted with caution.
Table (ref) reports $p$-values from testing the null hypothesis of a zero effect:
The null hypothesis (ref) is rejected at the 10% level based on both permutation schemes and all three methods.
Figure (ref) displays pointwise 90% confidence intervals. The results are similar for all three methods. While the effect was not or only marginally significant during the first three years, legalizing indoor prostitution significantly decreased the incidence of female gonorrhea thereafter, corroborating the findings by cunningham2018decriminalizing.
To investigate the robustness of our results, we perform a leave-one-out robustness check abadie2015comparative to assess whether our findings are driven by a single control state. We iteratively exclude from the control group one of the states for which either the SC or constrained Lasso weights estimated based on the pre-treatment data are non-zero and compute the $p$-values for testing hypothesis (ref). Figure (ref) displays the distribution of the resulting $p$-values. Overall, our results are robust and not driven by a single control state: except for one specification, all results are significant at the 10%-level.
\if11 {
We are grateful to Guido Imbens, Jacopo Diquigiovanni, Bruno Ferman, the Co-Editor (Matias Cattaneo), anonymous referees, and many seminar and conference participants for valuable comments. We would like to thank Scott Cunningham and Manisha Shah for sharing the data for the empirical application. W\"uthrich is also affiliated with CESifo and the Ifo Institute. Victor Chernozhukov gratefully acknowledges funding by the National Science Foundation. All errors are our own.
} \fi
\if01 {
} \fi
\spacingset{1.0} {2pt}