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.
114,307 characters · 18 sections · 49 citation commands
Prediction Intervals for Synthetic Control Methods
Keywords: causal inference, synthetic controls, prediction intervals, non-asymptotic inference.
\thispagestyle{empty}
\doublespacing \setcounter{page}{1} \pagestyle{plain}
\pagestyle{plain}
The synthetic control (SC) method was first introduced by Abadie-Gardeazabal_2003_AER as an approach to study the causal effect of a treatment affecting a single aggregate unit that is observed both before and after the treatment occurs. The authors originally motivated the method with a study of the effect of terrorism in the Basque Country on its GDP per capita. The Basque Country was one of the three richest regions in Spain before the outset of terrorism around the mid-1970s, but the region became relatively poorer in the decades that followed. The question is whether this relative decline can be attributed to terrorism. Their analysis covers the 1955-2000 period and places the beginning of intense terrorism in 1975, thus defining a “pre-treatment” period when terrorism is not salient (roughly 1955-1975), and a “post-treatment” period that starts when terrorism intensifies (roughly 1975 onward). The time series data allows for a comparison of Basque GDP before and after the onset of terrorism, but to interpret this change as the causal effect of terrorism would require assuming the absence of time trends. Instead, Abadie-Gardeazabal_2003_AER proposed to use other regions in Spain, whose GDP is also observed before and after the onset of terrorism in the Basque Country, to build an aggregate or “synthetic” control unit that captures the GDP trajectory that would have occurred in the Basque Country if terrorism had never occurred. The synthetic control is built as a weighted average of all units in the control group (the “donor pool”), where the weights are chosen so that the synthetic control's outcome in the pre-treatment period closely matches the treated unit's trajectory while also satisfying some constraints such as being non-negative, adding up to one, and accounting for other pre-treatment covariates. For a contemporaneous review of this literature, see Abadie_2021_JEL and the references therein.
The SC method has received increasing attention since its introduction, and is now a popular component of the methodological toolkit for causal inference and program evaluation Abadie-Cattaneo_2018_ARE. Methodological and theoretical research concerning SC methods has mostly focused either on expanding the SC causal framework (e.g., to dissagregated data or staggered treatment adoption settings) or on developing new implementations of the SC prediction (e.g., via different penalization constraints or matrix completion methods). Recent examples include Abadie-LHour_2021_JASA, Agarwal-Shah-Shen-Song_2021_JASA, Athey-et-al_2021_JASA, Bai-Ng_2021_JASA, BenMichael-Feller-Rothstein_2021_JASA, Chernozhukov-Wuthrich-Zhu_2021_ttest, Ferman_2021_JASA, Kellogg-Mogstad-Pouliot-Torgovitsky_2021_JASA, and Masini-Medeiros_2021_JASA; see their references for many more. In contrast, considerably less effort has been devoted to develop principled statistical inference procedures for uncertainty quantification within the SC framework. In particular, Abadie-Diamond-Hainmueller_2010_JASA propose a design-based permutation approach under additional assumptions, Li_2020_JASA relies on large-sample approximations for disaggregated data under correct specification, Chernozhukov-Wuthrich-Zhu_2021_JASA develop time-series permutation-based inference methods, and Shaikh-Toulis_2021_JASA discuss cross-sectional permutation-based inference methods in semiparametric duration-type settings. See also Feng_2021_wp, and references therein, for related large sample inference methods employing local principal component analysis based on nearest-neighbor approximations in possibly non-linear factor model settings.
We develop conditional prediction intervals for the SC framework, offering an alternative (conditional) inference method to assess statistical uncertainty. Our proposed approach builds on ideas from the literature on conditional prediction intervals Vovk_2012_ACML,Chernozhukov-Wuthrich-Zhu_2021_distributional and non-asymptotic concentration Vershynin_2018_Book,Wainwright_2019_Book in probability and statistics. As a consequence, the resulting (conditional) prediction intervals are conservative but formally shown to offer probability guarantees. We focus on uncertainty quantification via (conditional) prediction intervals because, in the SC framework, the treatment effect estimator is a random variable emerging from an out-of-sample prediction problem, based on the estimated SC weights constructed using pre-treatment data. Our inference procedures are not confidence intervals in the usual sense (i.e., giving a region in the parameter space for a non-random parameter of interest), but rather intervals describing a region on the support of a random variable where a new realization is likely to be observed.
Our construction begins by noting that the statistical uncertainty of the SC prediction is governed by two distinct sources of randomness: one due to the construction of the (likely misspecified) SC weights in the pre-treatment period, and the other due to the unobservable stochastic error in the post-treatment period when the treatment effect is analyzed. Accordingly, our proposed prediction intervals are constructed taking into account both sources of randomness. For the first source of uncertainty, we propose a simulation-based approach that is justified via non-asymptotic probability concentration and hence enjoys probability guarantees. This approach takes into account the specific construction of the SC weights. For the second source of uncertainty, which comes from out-of-sample prediction due to the unobservable error in the post-treatment period, we discuss several approaches based on nonparametric and parametric probability approximations as a framework for principled sensitivity analysis. This second uncertainty source is harder to handle nonparametrically, and hence its contribution to the overall prediction intervals should be considered with care. Our approach in this paper is to employ an agnostic sensitivity analysis, but future work will consider other approaches.
Our results are obtained under high-level conditions, but we provide primitive conditions for three examples: an outcomes-only setting with i.i.d. data, a multi-equation setting allowing for stationary weakly dependent data where the weights are obtained by not only matching the pre-treatment trends of the outcome of interest but also approximating the trajectories of additional variables such as important covariates or secondary outcomes; and a non-stationary cointegration setting. All three settings allow the weights to be covariate-adjusted in each equation. We also showcase our methods numerically, using both simulated and real data. The methods perform well in finite samples.
The rest of the paper proceeds as follows. Section (ref) provides a formal introduction to the SC framework and defines the basic quantities of interest. Section (ref) introduces the prediction intervals we focus on, and provides basic intuition for their decomposition in terms of the SC weights estimation error and the unobservable post-treatment error. Section (ref) develops a simulation-based method to account for the first source of uncertainty, and Section (ref) discusses how to (model and) account for the second source of uncertainty. Section (ref) illustrates the performance of our proposed prediction intervals with a Monte Carlo experiment and two empirical examples from the SC literature. Section (ref) concludes. Appendix (ref) provides an extension of our main in-sample uncertainty quantification approach to the case of weakly dependent ($\beta$-mixing) stationary time series data. All the proofs of our technical results, as well as additional numerical evidence, are collected in the online supplemental appendix. We provide companion replication codes in R, and a general-purpose software package is underway \citep*{Cattaneo-Feng-Palomba-Titiunik_2021_scpi}.
We consider the standard synthetic control framework with a single treated unit and several control units, allowing for both stationary and non-stationary data. The data may include only the outcome of interest, or the outcome of interest plus other variables. The researcher observes $N+1$ units for $T_0+T_1$ periods of time. Units are indexed by $i = 1,2,\ldots N, N+1$, and time periods are indexed by $t=1,2, \ldots,T_0, T_0+1, \ldots, T_0+T_1$. During the first $T_0$ periods, all units are untreated. Starting at $T_0 + 1$, unit $1$ receives treatment but the other units remain untreated. Once the treatment is assigned at $T_0+1$, there is no change in treatment status: the treated unit continues to be treated and the untreated units remain untreated until the end of the series, $T_1$ periods later.
Each unit $i$ at period $t$ has two potential outcomes, $Y_{it}(1)$ and $Y_{it}(0)$, respectively denoting the outcome under treatment and the outcome in the absence of treatment (which we call the control or the untreated condition). This notation imposes two additional implicit assumptions that are standard in this setting: no spillovers (the potential outcomes of unit $i$ depend only on $i$'s treatment status) and no anticipation (the potential outcomes at $t$ depend only on the treatment status of the same period).
Attention is restricted to the impact of the treatment on the treated unit. By treatment impact, we mean the difference between the outcome path taken by the treated unit, and the path it would have taken in the absence of the treatment. The quantity of interest is
where $\tau_{t}$ may be regarded as random or non-random depending on the framework considered. In this paper, we view $\tau_{t}$ as a random variable.
For each unit, we only observe the potential outcome corresponding to the treatment status actually received by the unit. We denote the observed outcome by $Y_{it}$, which is defined as \[ Y_{it}=
. \]
This means that, in $\tau_t$, the treated unit's potential outcome $Y_{1t}(0)$ is unobservable for all $t > T_0$. The idea of the synthetic control method is to use an appropriate combination of the post-treatment observed outcomes of the untreated units to approximate the treated unit's counterfactual post-treatment outcome, $Y_{1t}(0)$ for $t>T_0$. This idea has been formalized in different ways since it was originally proposed by Abadie-Gardeazabal_2003_AER.
In all SC frameworks, the formalization chooses a set of weights $\mathbf{w} = (w_2,w_3,\dots,w_{N+1})'$ such that a given loss function is minimized under constraints. Given a set of estimated weights $\widehat{\mathbf{w}}$, the treated unit's counterfactual predicted outcome is then calculated as $\widehat{Y}_{1t}(0) = \sum_{i=2}^{N+1} \widehat{w}_{i}Y_{it}(0) $ for $t>T_0$. The weighted average $\widehat{Y}_{1t}(0)$ is often referred to as the synthetic control of the treated unit, as it represents how the untreated units can be combined to provide the best counterfactual for the treated unit in the post-treatment period.
When the data contains only information on the outcome of interest, $\mathbf{w}$ is chosen such that the weighted average of the outcomes of the untreated units approximates well the outcome trajectory of the treated unit in the period before the treatment. That is, the weights $\mathbf{w}$ are chosen so that \[\sum_{i=2}^{N+1} w_{i}Y_{it}(0) \approx Y_{1t}(0), \qquad \text{for}\quad t=1,2,\dots,T_0,\] where the meaning of the symbol “$\approx$” varies depending on the specific framework considered. A leading example constrains the weights to be non-negative and sum to one, and estimates $\mathbf{w}$ by constrained least squares:
where $r$ denotes the intercept, and $\mathcal{W}$ and $\mathcal{R}$ denote the corresponding constraint (or feasibility) sets---we give formal definitions in the next subsection.
When the weights are chosen according to ((ref)), the resulting synthetic control will reproduce as closely as possible the outcome trajectory of the treated unit in the pre-treatment period. For example, in the Basque terrorism application, this procedure would lead to a synthetic Basque Country that would have a similar per capita GDP to the Basque Country's per capita GDP in the 1955-1975 period when terrorism is not salient.
This outcomes-only version of the SC method, however, cannot guarantee that the resulting synthetic control unit will be similar to the treated unit in any characteristics other than the (pre-treatment) outcome. In some applications, this feature may be undesirable, as researchers may have access to additional characteristics such as baseline covariates or secondary outcomes and may want to also ensure that the synthetic control approximates the treated unit in terms of these additional characteristics. The SC framework can handle this case by including additional equations for these additional characteristics and minimizing the combined loss. In this case, letting $l=1,2,\ldots, M$ index the variables that will be “matched” to produce the weights, the minimization problem above can be generalized as
where $\widehat{\mathbf{r}}=(\widehat{r}_1,\dots,\widehat{r}_M)'$ and $\{\upsilon_{t,l}\}_{1\leq t\leq T_0, 1\leq l\leq M}$ are positive constants reflecting the relative importance of different equations and periods.
For example, in the original Basque terrorism example, Abadie-Gardeazabal_2003_AER show that the Basque country differs from the rest of Spain in terms of population density, and they are concerned that pre-terrorism differences in population density may affect economic growth in the post-treatment period. In this case, we can choose the weights $\widehat{\mathbf{w}}$ to ensure not only that the per capita GDP trajectory is similar between the treated unit and the synthetic control unit, but also to ensure that the synthetic control is similar to the treated unit in terms of population density. To implement this multi-equation SC method, we fit equation ((ref)) with two variables ($M=2$) where $Y_{it,1}$ ($l=1$) is per capita GDP for region $i$ in year $t$ and $Y_{it,2}$ ($l=2$) is population density for region $i$ in year $t$. When $\widehat{\mathbf{w}}$ is chosen this way, the resulting synthetic control will resemble (to the extent that the data allows) the treated unit in terms of both per capita GDP and population density.
Equation (ref) can be viewed as a (weighted) combination of $M$ optimization problems in (ref), satisfying an additional constraint that the weights $\mathbf{w}$ must be the same across the $M$ equations. For simplicity, we let $\upsilon_{t,l}=1$ for all $t$ and $l$, but the analysis below can be applied to the more general case if additional regularity conditions are imposed on $\{\upsilon_{t,l}\}_{1\leq t\leq T_0, 1\leq l\leq M}$.
The two cases just discussed (outcomes-only and multi-equation SC frameworks) allow for weakly dependent and cointegrated data, and they also can be generalized further by including covariates in a linear and additive way in ((ref)) or ((ref)). This covariate adjustment would introduce additional parameters to the fit that would not be of primary interest; rather, they would be included to “partial out” the effect of additional covariates.
We now introduce a general framework and further notation that encompass and formalize the two particular examples discussed above as well as other synthetic control approaches in the literature. Our general framework includes the outcomes-only fit and the multi-equation fit (i.e., outcome plus other variables) as particular cases, allowing for covariate adjustment and non-stationary data in a unified way.
Consider synthetic control weights constructed simultaneously for $M$ features of the treated unit, denoted by $\mathbf{A}_l=(a_{1,l}, \cdots, a_{T_0,l})'\in\mathbb{R}^{T_0}$, with index $l=1,\cdots, M$. For each feature $l$, there exist $J+K$ variables that can be used to predict or “match” the $T_0$-dimensional vector $\mathbf{A}_l$. These $J+K$ variables are separated into two groups denoted by $\mathbf{B}_l=(\mathbf{B}_{1,l}, \mathbf{B}_{2,l}, \cdots, \mathbf{B}_{J, l})\in\mathbb{R}^{T_0\times J}$ and $\mathbf{C}_l=(\mathbf{C}_{1,l}, \cdots, \mathbf{C}_{K,l})\in\mathbb{R}^{T_0\times K}$, respectively. More precisely, for each $j$, $\mathbf{B}_{j,l}=(b_{j1,l}, \cdots, b_{jT_0,l})'$ corresponds to the $l$th feature of the $j$th unit observed in $T_0$ pre-treatment periods and, for each $k$, $\mathbf{C}_{k,l}=(c_{k1,l}, \cdots, c_{kT_0,l})'$ is another vector of control variables also possibly used to predict $\mathbf{A}_l$ over the same pre-intervention time span. For ease of notation, we let $d=J+KM$.
The goal of the synthetic control method is to search for a vector of common weights $\mathbf{w}\in\mathcal{W}\subseteq\mathbb{R}^{J}$ across the $M$ features and a vector of coefficients $\mathbf{r}\in\mathcal{R}\subseteq\mathbb{R}^{KM}$, such that the linear combination of $\mathbf{B}_l$ and $\mathbf{C}_l$ “matches” $\mathbf{A}_l$ as close as possible, for all $1\leq l\leq M$. This goal is typically achieved via the following optimization problem:
where \[ \mathbf{A}=
,\quad \mathbf{B}=
, \quad \mathbf{C}=
, \] and where the feasibility sets $\mathcal{W}$ and $\mathcal{R}$ capture the restrictions imposed. (For simplicity we do not introduce an explicit re-weighting of the $M$ equations, but recall that this extension is possible.) This framework encompasses multiple prior synthetic control formalizations in the literature, which differ in whether they include additional covariates, whether the data is assumed to be stationary, and the particular choice of constraint sets $\mathcal{W}$ and $\mathcal{R}$ used, among other possibilities.
The following list provides some examples of different constraint sets used in practice, where $\|\cdot\|_p$ denotes the $L_p$ vector norm and $Q$ and $\alpha$ are tuning parameters.
In some applications the intercept in (ref) is removed by demeaning the data before the analysis. Section (ref) discusses in detail the outcomes-only case, as well as the multi-equation case where the researcher “matches” on pre-treatment characteristics and pre-intervention outcomes simultaneously. That section also deals with stationary weakly dependant data, and non-stationary data (i.e., cointegration system).
For example, the outcomes-only setup can be obtained as a particular case of (ref) with $M=1$ (there is only one feature to match on), $J=N$ (there are $N$ units in the donor pool), and $K=1$ (there is an intercept). Then, $\mathbf{A}_1=(Y_{11}, Y_{12}, \cdots, Y_{1T_0})'$, $\mathbf{B}_{j,1}=(Y_{(j+1)1}, Y_{(j+1)2},\cdots,Y_{(j+1)T_0})'$, $\mathbf{C}_{j,1}=(1, 1,\cdots,1)'$, and (ref) reduces to the (possibly constrained) optimization problem (ref). The multi-equation setup with one outcome and one covariate can be obtained similarly by setting $M=2$ (there are two features to match on), $J=N$ ( $N$ units in the donor pool), and $K=1$ (there is an intercept), which reduces to ((ref)).
To further understand our proposed inference approach, we define the pseudo-true values $\mathbf{w}_0$ and $\mathbf{r}_0$ relative to a sigma field $\mathscr{H}$:
and thus write
where $\mathbf{U}=(u_{1,1}, \cdots, u_{T_0,1}, \cdots, u_{1,M}, \cdots, u_{T_0, M})'\in\mathbb{R}^{T_0M}$ is the corresponding pseudo-true residual relative to a sigma field $\mathscr{H}$. That is, $\mathbf{w}_0$ and $\mathbf{r}_0$ are the mean square error estimands associated with the (possibly constrained) best linear prediction coefficients $\widehat{\mathbf{w}}$ and $\widehat{\mathbf{r}}$ conditional on $\mathscr{H}$. Importantly, we do not attach any structural meaning to equation (ref). The population vectors $\mathbf{w}_0$ and $\mathbf{r}_0$ are (conditional) pseudo-true values whose meaning should be understood in context, and are determined by the assumptions imposed on the data generating process. In particular, with strong parametric functional form assumptions or rich enough nonparametric basis expansions, equation (ref) may be viewed as a representation (or approximation) of $\mathbb{E}[\mathbf{A}|\mathbf{B},\mathbf{C},\mathscr{H}]$. In such cases, $\mathbb{E}[\mathbf{U}|\mathbf{B},\mathbf{C},\mathscr{H}]=\bm{0}$ or, at least, $\mathbb{E}[\mathbf{U}|\mathbf{B},\mathbf{C},\mathscr{H}]$ is taken to be “small”. Alternatively, if the (population, conditional) linear projection coefficients lie on $\mathcal{W}\times\mathcal{R}$, that is, the constraints imposed by $\mathcal{W}$ and $\mathcal{R}$ in (ref) are not binding, then (ref) represents the best linear approximation of $\mathbf{A}$ based on $(\mathbf{B},\mathbf{C})$, conditional on $\mathscr{H}$. In this scenario, $\mathbf{U}$ is uncorrelated with $(\mathbf{B},\mathbf{C})$, conditional on $\mathscr{H}$. Most importantly, in general, $\mathbf{U}$ may not be mean zero due to the (binding) constraints imposed in the construction of $\widehat{\mathbf{w}}$ and $\widehat{\mathbf{r}}$.
Given estimated weights $\widehat{\mathbf{w}}$ and coefficients $\widehat{\mathbf{r}}$, the post-treatment counterfactual outcome for the treated unit is predicted by \[\widehat{Y}_{1T}(0) = \mathbf{x}_T'\widehat{\mathbf{w}}+\mathbf{g}_T'\widehat{\mathbf{r}} = \mathbf{p}_T'\widehat{\boldsymbol{\beta}}, \qquad \mathbf{p}_T := (\mathbf{x}_T', \mathbf{g}_T')', \qquad T>T_0, \] where $\mathbf{x}_T\in\mathbb{R}^{N}$ is a vector of predictors for control units observed in time $T$ and $\mathbf{g}_T\in\mathbb{R}^{KM}$ is another set of user-specified predictors observed at time $T$. Variables included in $\mathbf{x}_T$ and $\mathbf{g}_T$ need not be the same as those in $\mathbf{B}$ and $\mathbf{C}$, but will be part of the sigma field $\mathscr{H}$, as explained in more detail in the next section. Therefore, from this perspective, our focus is on conditional inference. We decompose the potential outcome of the treated unit accordingly:
where $e_T$ is defined by construction. In our analysis, $\mathbf{w}_0$ and $\mathbf{r}_0$ are assumed to be possibly random elements around which $\widehat{\mathbf{w}}$ and $\widehat{\mathbf{r}}$ are concentrating in probability, respectively, which is why we called them pseudo-true values.
The distance between the estimated treatment effect on the treated and the target population one is
Within the synthetic control framework, we view the quantity of interest $\tau_{T}$ as a random variable, and hence we refrain from calling it a “parameter”. Consequently, we call $\widehat{\tau}_{T}$ a prediction of $\tau_{T}$ rather than an “estimator” of it, and focus on building prediction intervals rather than confidence intervals.
Given the generic framework introduced in the previous section, we now present our proposed prediction intervals for $\tau_T$. See Vovk_2012_ACML, Chernozhukov-Wuthrich-Zhu_2021_distributional,Chernozhukov-Wuthrich-Zhu_2021_JASA, and references therein, for recent papers on (conditional) prediction intervals and related methods. Let $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ be random quantities defined on a probability space $(\Omega, \mathscr{F}, \mathbb{P})$, and $\mathscr{H}\subseteq\mathscr{F}$ be a sub-$\sigma$-field. For some $\alpha, \pi\in(0,1)$, we say a random interval $\mathcal{I}$ is an $(\alpha,\pi)$-valid $\mathscr{H}$-conditional prediction interval for $\tau_T$ if
If $\mathscr{H}$ is the trivial $\sigma$-field over $\Omega$, then $\mathcal{I}$ reduces to an unconditional prediction interval for $\tau_T$. In the general case, $\mathcal{I}$ is an $\mathscr{H}$-conditionally $(\alpha,\pi)$-valid prediction interval: the conditional coverage probability of $\mathcal{I}$ is at least $(1-\alpha)$, which holds with probability over $\mathscr{H}$ at least $ (1-\pi)$. In practice, $(1-\alpha)$ is a desired confidence level chosen by users, say $95\%$, and $\pi$ is a “small” number that depends on the sample size and typically goes to zero in some asymptotic sense. In this paper, all the results are valid for all $T_0$ large enough, with the associated probability loss $\pi$ characterized precisely. Thus, we say that the conditional coverage of the prediction interval $\mathcal{I}$ is at least $(1-\alpha)$ with high probability, or that the conditional prediction interval offers finite-sample probability guarantees. Our results imply $\pi\to 0$ as $T_0\to\infty$, but no limits or asymptotic arguments are used in this paper.
An asymptotic analogue to the above definition (ref) would be $\mathbb{P}(\tau_T\in\mathcal{I}|\mathscr{H})\geq 1-\alpha-o_\mathbb{P}(1)$ or, perhaps, $\mathbb{P}(\tau_T\in\mathcal{I}|\mathscr{H})\to_\mathbb{P} 1-\alpha$, where the probability limit is taken as the sample size grows to infinity (e.g., as $T_0\to\infty$). In this case, we say $\mathcal{I}$ is an $\mathscr{H}$-conditional prediction interval for $\tau_T$ that is asymptotically valid with coverage probability (at least) $(1-\alpha)$. This is a weaker property because it does not offer any finite-sample probability guarantees for the (conditional) coverage of the prediction interval.
We employ the following lemma to construct valid, conditional prediction intervals in the sense of (ref). This lemma follows from the union bound applied to $\widehat{\tau}_T - \tau_{T} = \mathbf{p}_T'(\boldsymbol{\beta}_0-\widehat{\boldsymbol{\beta}}) + e_T$.
This lemma provides a simple way to construct an $\mathscr{H}$-conditional prediction interval enjoying $(\alpha,\pi)$-validity with $\alpha=\alpha_1+\alpha_2$ and $\pi=\pi_1+\pi_2$: \[\mathcal{I} = \Big[\widehat{\tau}_T-M_{1,\mathtt{U}}-M_{2,\mathtt{U}} \, , \, \widehat{\tau}_T-M_{1,\mathtt{L}}-M_{2,\mathtt{L}}\Big],\] for appropriate choices of $M_{1,\mathtt{L}}$, $M_{1,\mathtt{U}}$, $M_{2,\mathtt{L}}$ and $M_{2,\mathtt{U}}$ and conditioning sigma field. In this paper, we consider conditional prediction intervals with $\mathscr{H}= \sigma(\mathbf{B}, \mathbf{C}, \mathbf{x}_T, \mathbf{g}_T)$ and focus on building a probability bound for each of the two terms, $\mathbf{p}_T'(\boldsymbol{\beta}_0-\widehat{\boldsymbol{\beta}})$ and $e_T$, separately, and then combine them to build an overall probability bound via Lemma (ref). In the decomposition leading to the prediction interval construction, we interpret $\mathbf{p}_T'(\boldsymbol{\beta}_0-\widehat{\boldsymbol{\beta}})$ as capturing the in-sample uncertainty coming from constructing the SC weights using pre-treatment information, and $e_T$ the out-of-sample uncertainty coming from misspecification along with any additional noise occurring at the post-treatment period $T>T_0$. The next two subsections are devoted to handle each of these terms, respectively.
We first quantify the in-sample uncertainty coming from $\mathbf{p}_T'(\boldsymbol{\beta}_0-\widehat{\boldsymbol{\beta}})$, thereby providing methods to determine $(M_{1,\mathtt{L}},M_{1,\mathtt{U}})$ and their probability guarantees $(\alpha_1,\pi_1)$ in Lemma (ref). Let $\mathbf{Z}=(\mathbf{B}, \mathbf{C})$, $\mathbf{D}$ be a non-negative diagonal (scaling) matrix of size $d$, possibly depending on the pre-treatment sample size $T_0$, and recall that $\mathscr{H} = \sigma(\mathbf{B}, \mathbf{C}, \mathbf{x}_T, \mathbf{g}_T)$. Because $\widehat{\boldsymbol{\beta}}$ solves (ref), we can define $\widehat{\boldsymbol{\delta}}:=\mathbf{D}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_0)$ as the optimizer of the centered criterion function: \[ \widehat{\boldsymbol{\delta}}=\underset{\boldsymbol{\delta}\in\Delta}{\operatorname*{arg\,min}}\; \big\{\boldsymbol{\delta}'\widehat{\mathbf{Q}}\boldsymbol{\delta}-2\widehat{\boldsymbol{\gamma}}'\boldsymbol{\delta}\big\}, \] where $\widehat{\mathbf{Q}}=\mathbf{D}^{-1}\mathbf{Z}'\mathbf{Z}\mathbf{D}^{-1}$, $\widehat{\boldsymbol{\gamma}}'=\mathbf{U}'\mathbf{Z}\mathbf{D}^{-1}$, and $\Delta=\{\bm{h}\in\mathbb{R}^d: \bm{h}=\mathbf{D}(\boldsymbol{\beta}-\boldsymbol{\beta}_0),\, \boldsymbol{\beta}\in\mathcal{W}\times \mathcal{R}\}$.
The following lemma, which holds whether or not $\boldsymbol{\gamma}:=\mathbb{E}[\widehat{\boldsymbol{\gamma}}|\mathscr{H}] = \mathbf{0}$, is a key building block for our prediction interval construction.
This lemma does not involve probabilistic statements, but rather follows from basic features of constrained least squares optimization. In particular, simple bounds on $\mathbf{p}_T'\mathbf{D}^{-1}\widehat{\boldsymbol{\delta}}$ can be deduced based on the basic inequality from optimization $\widehat{\boldsymbol{\delta}}'\widehat{\mathbf{Q}}\widehat{\boldsymbol{\delta}}-2(\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma})'\widehat{\boldsymbol{\delta}}\leq 0$, and the fact that any solution must satisfy the constraints imposed in (ref), i.e., $\widehat{\boldsymbol{\delta}}\in\Delta$. The second part of the lemma establishes that the set of possible localization values ($\boldsymbol{\xi}$) determining the feasibility set ($\mathcal{M}_{\boldsymbol{\xi}}$) of the bounding (random) quantities $\varsigma_{\mathtt{L}}$ and $\varsigma_{\mathtt{U}}$ in Lemma (ref) (when $\boldsymbol{\xi}=\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}$) form (random) convex sets.
Conditional on $\mathscr{H}$, the set $\mathcal{M}_{\boldsymbol{\xi}}$ is not random due to $\widehat{\mathbf{Q}}$, which is random only unconditionally. As a consequence, conditional on $\mathscr{H}$, both $\mathcal{M}_{\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}}$ and $\{\mathbf{p}_T'\mathbf{D}^{-1}\boldsymbol{\delta}:\boldsymbol{\delta}\in\mathcal{M}_{\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}}\}$ are random sets only because $\widehat{\boldsymbol{\gamma}}$ is a random quantity and, accordingly, $\varsigma_{\mathtt{L}}$ and $\varsigma_{\mathtt{U}}$ are random variables defined by a random set. If the conditional distributions of $\varsigma_{\mathtt{L}}$ and $\varsigma_{\mathtt{U}}$ were known, we could take their quantiles as lower and upper bounds for the quantiles of the conditional distribution of $\mathbf{p}_T'\mathbf{D}^{-1}\widehat{\boldsymbol{\delta}}=\mathbf{p}_T'(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_0)$, thereby transforming the first conclusion of Lemma (ref) into a probabilistic statement. However, this approach requires knowledge of the conditional (on $\mathscr{H}$) distribution of the bounding random variables $\varsigma_{\mathtt{L}}$ and $\varsigma_{\mathtt{U}}$. The convexity properties also established in Lemma (ref) allow us to provide precise bounds on the desired conditional distribution of the bounding random variables using Berry-Esseen bounds for convex sets Raivc_2019_Bernoulli.
The following theorem formalizes our first main result based on Lemma (ref). We only present the result for the upper bound to conserve space, but the analogous result holds for the lower bound. See Remarks SA-2.1, SA-2.2 and SA-2.3 in the supplemental appendix for more details. Let $\boldsymbol{\Sigma}=\mathbb{V}[\widehat{\boldsymbol{\gamma}}|\mathscr{H}]$ and $\boldsymbol{\Sigma}^{-1/2}\mathbf{D}^{-1}\mathbf{Z}'=(\tilde{\mathbf{z}}_{1,1}, \cdots, \tilde{\mathbf{z}}_{T_0,1}, \cdots, \tilde{\mathbf{z}}_{1,M}, \cdots, \tilde{\mathbf{z}}_{T_0,M})$. In addition, let $\|\cdot\|$ denote the spectral matrix norm (so that $\|\cdot\|=\|\cdot\|_2$ for vectors).
This theorem is established under two high-level conditions. Condition (T1.i) imposes independence across time for the pseudo-residuals underlying the population analogue construction of the synthetic control weights in (ref). In Appendix (ref) we relax this requirement by allowing for weak dependence across time via a $\beta$-mixing condition Doukhan_2012_Book, but to avoid untidy conditions we focus on the independent case here. Importantly, even in this case, Theorem (ref) covers non-stationarity in the outcome variable (via a cointegration relationship). See Section (ref) for different examples with independent, weakly stationary, and non-stationary data.
The second high-level requirement in Theorem (ref) helps control the (distributional) distance between $\widehat{\boldsymbol{\gamma}}-\boldsymbol{\gamma}$ and the Gaussian random vector $\mathsf{N}(\mathbf{0}, \boldsymbol{\Sigma})$, conditionally on $\mathscr{H}$, as well as the unconditional probability loss $\pi_\gamma$. Condition (T1.ii) can be verified in a variety of ways depending on the dependence structure imposed on the data and other regularity conditions, as we illustrate in Section (ref). For instance, two sufficient conditions are: $\max_{1\leq l\leq M}\max_{1\leq t\leq T_0} \mathbb{E}[|u_{t,l}-\mathbb{E}[u_{t,l}|\mathscr{H}]|^3|\mathscr{H}]\leq \eta$, a.s. on $\mathscr{H}$ for some constant $\eta>0$, and $\mathbb{P}\{\sum_{t=1}^{T_0}\sum_{l=1}^{M}\|\tilde{\mathbf{z}}_{t,l}\|^3\leq\epsilon_\gamma(42(d^{1/4}+16)\eta M^2)^{-1}\}\geq 1-\pi_\gamma$.
Since $\mathbf{p}_T'\mathbf{D}^{-1}\widehat{\boldsymbol{\delta}}=\mathbf{p}_T'(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_0)$, Theorem (ref) could immediately be applied to construct valid $M_{1,\mathtt{L}}$ and $M_{1,\mathtt{U}}$ in Lemma (ref) if $\boldsymbol{\Sigma}=\mathbb{V}[\widehat{\boldsymbol{\gamma}}|\mathscr{H}]$ was known. Thus, to finalize the in-sample uncertainty quantification we discuss a feasible simulation-based approximation for the critical value $\mathfrak{c}^\dagger(1-\alpha)$. To describe such approach, define a simulation-based criterion function conditional on the data \[\ell^{\star}(\boldsymbol{\delta}) = \boldsymbol{\delta}'\widehat{\mathbf{Q}}\boldsymbol{\delta}-2(\mathbf{G}^{\star})'\boldsymbol{\delta},\qquad \mathbf{G}^{\star}\thicksim\mathsf{N}(0,\widehat{\boldsymbol{\Sigma}}),\] where $\widehat{\boldsymbol{\Sigma}}$ is some estimate of $\boldsymbol{\Sigma}$. The form of $\widehat{\boldsymbol{\Sigma}}$ depends on the specific dependence structure underlying the data and related regularity conditions, as we illustrate in Section (ref). Naturally, the important high-level requirement is that $\widehat{\boldsymbol{\Sigma}}$ should concentrate around $\boldsymbol{\Sigma}$ with known probability; see Theorem (ref) below for the precise statement. In addition, the constraint set used in the simulation has to be properly defined to account for the parameters being possibly near or at the boundary, so that it mimics the local geometry of $\Delta$. Specifically, let $\Delta^{\star}$ denote the constraint set used in simulation. We require that
where $\mathcal{B}(\mathbf{0}, \varepsilon)$ is an $\varepsilon$-neighborhood around zero. We say $\Delta^{\star}$ is locally equal to $\Delta$ if (ref) is satisfied. Consequently, searching for the desired region under constraints in $\Delta^{\star}$ is almost equivalent to doing so under constraints in $\Delta$. We discuss below more implementation details.
The next theorem establishes the validity of our proposed simulation-based inference method and provides the associated probability guarantees, under high-level conditions. Let $\|\cdot\|_\mathtt{F}$ denote the Frobenius matrix norm (so that $\|\cdot\|_\mathtt{F}=\|\cdot\|=\|\cdot\|_2$ for vectors), and $\mathbf{I}_q$ the identity matrix of size $q$ for an integer $q>0$.
This theorem gives a feasible, simulation-based approach to determine valid $M_{1,\mathtt{L}}$ and $M_{1,\mathtt{U}}$ in Lemma (ref), with precise coverage probability guarantees. The first high-level Condition (T2.i) in Theorem (ref) takes as a starting point the conclusion of Theorem (ref) or, alternatively, the conclusion of Theorem (ref) in the appendix when the data is assumed to exhibit weak dependence via a $\beta$-mixing condition. The other three high-level conditions in Theorem (ref) are intuitive. Conditions (T2.ii) and (T2.iii) control the local geometry of the simulation feasibility set, as discussed earlier, while Condition (T2.iv) requires $\widehat{\boldsymbol{\Sigma}}$ to be a “good” approximation of $\boldsymbol{\Sigma}$, in the sense that $\widehat{\boldsymbol{\Sigma}}$ concentrates in probability around $\boldsymbol{\Sigma}$ with well-controlled errors. Importantly, Theorem (ref) is carefully crafted to accommodate both Theorem (ref) (independent data) and Theorem (ref) in the appendix (weakly dependent time series data) in a unified way. The next section illustrates different cases with practically relevant examples, and gives precise primitive conditions.
We consider the standard synthetic control constraints $\mathcal{W}=\{\mathbf{w}\in\mathbb{R}^N_+: \|\mathbf{w}\|_1=1\}$ and $\mathcal{R}=\mathbb{R}^{KM}$. For simulation-based inference, we define explicitly a relaxed constraint set based on the original estimated coefficients $\widehat{\boldsymbol{\beta}}$: $\Delta^{\star}= \{\mathbf{D}(\boldsymbol{\beta}-\widehat{\boldsymbol{\beta}}^{\star}):\boldsymbol{\beta}=(\mathbf{w}', \mathbf{r}')', \mathbf{w}\in\mathbb{R}^N_+,\; \|\mathbf{w}\|_1=\|\widehat{\mathbf{w}}^{\star}\|_1\}$, where $\widehat{\boldsymbol{\beta}}^{\star}=(\widehat{\mathbf{w}}^{\star '},\widehat{\mathbf{r}}')'$, $\widehat{\mathbf{w}}^{\star}=(\widehat{\omega}_2^{\star}, \cdots, \widehat{\omega}_{N+1}^{\star})'$, $\widehat{\omega}_j^{\star}=\widehat{\omega}_j\mathds{1}(|\widehat{\omega}_j|>\varrho)$, and $\varrho$ is a tuning parameter that ensures the constraint set in the simulation world preserves the local geometry of $\Delta$. Moreover, we set $\mathbf{x}_T=(Y_{2T}(0), \cdots, Y_{(N+1)T}(0))'$ as it is common in the SC literature. Other synthetic control methods that vary these choices, including the other constraint sets $\mathcal{W}$ discussed previously, can be handled analogously, but we do not discuss them in this paper due to space limitations. Finally, in the remaining of this paper, we let $\mathfrak{C}$, $\mathfrak{C}^{\star}$ and $\mathfrak{c}$, with various sub-indexes, denote non-negative finite constants not depending on $T_0$. In simple cases, we give the exact expression of these constants, while in other cases they can be characterized from the proofs of the results. Let $\lambda_{\min}(\mathbf{M})$ and $\lambda_{\max}(\mathbf{M})$ be the minimum and the maximum eigenvalues of a generic square matrix $\mathbf{M}$.
We start with the simplest possible example already introduced in Section (ref). The SC weights are constructed based on past outcomes only, and the model allows for an intercept. Thus, the working model simplifies to \[a_{t}=\mathbf{b}_{t}'\mathbf{w}_{0} + r_0 + u_t, \qquad t=1, \cdots, T_0,\] where $a_{t}:=Y_{1t}(0)$, $\mathbf{b}_{t} :=(Y_{2t}(0), Y_{3t}(0), \ldots, Y_{(N+1)t}(0))'$, and with $M=1$, $K=1$, and $d=N+1$. Recall that $\mathbf{w}_0=(w_{0,1},w_{0,2},\dots,w_{0,J})'$ is defined in (ref), and let $\mathbf{z}_t=(\mathbf{b}_t', 1)'$, $\boldsymbol{\beta}_0=(\mathbf{w}_0', r_0)'$. We further assume independent sampling across time, and thus set $\mathbf{D}=T_0^{1/2}\mathbf{I}_d$. A natural variance estimator is \[\widehat{\boldsymbol{\Sigma}}=\frac{1}{T_0}\sum_{t=1}^{T_0}\mathbf{z}_t\mathbf{z}_t' (\widehat{u}_{t}-\widehat{\mathbb{E}}[u_{t}|\mathbf{b}_t])^2,\] where $\widehat{u}_t=a_t-\mathbf{z}_t'\widehat{\boldsymbol{\beta}}$, and $\widehat{\mathbb{E}}[u_{t}|\mathbf{b}_t]$ denotes some estimate of the conditional mean of the pseudo-residuals.
Theorem SA-1 in the supplemental appendix gives precise primitive conditions to verify the high-level conditions of Theorems (ref) and (ref). In particular, assuming that $\{\mathbf{z}_t, u_t\}_{t=1}^T$ is i.i.d. over $t=1, \cdots, T_0$, and that $\max_{1\leq t\leq T_0}\mathbb{E}[|u_t|^3|\mathbf{B}]\leq\bar{\eta}_1$ a.s. on $\sigma(\mathbf{B})$ and $\mathbb{E}[\|\mathbf{z}_t\|^6]\leq \bar{\eta}_2$, $\min_{1\leq t\leq T_0}\mathbb{V}[u_{t}|\mathbf{B}]\geq\underline{\eta}_1$ a.s. on $\sigma(\mathbf{B})$, and $\lambda_{\min}(\mathbb{E}[\mathbf{z}_t\mathbf{z}_t'])\geq \underline{\eta}_2$, for finite non-negative constants $\bar{\eta}_1$, $\bar{\eta}_2$, $\underline{\eta}_1$ and $\underline{\eta}_2$, we show that the conditions of Theorem (ref) hold with $\pi_{\gamma}=\mathfrak{C}_{\pi}T_0^{-1}$ and $\epsilon_{\gamma}=\mathfrak{C}_{\epsilon}T_0^{-1/2}$, where $\mathfrak{C}_{\pi}=\frac{d}{\bar{\eta}_2}+\frac{4d^4\bar{\eta}_2}{\underline{\eta}_2^2}$ and $\mathfrak{C}_{\epsilon}=42(d^{1/4}+16) \frac{2^{5/2}d^{3/2}\bar{\eta}_1\bar{\eta}_2}{(\underline{\eta}_1\underline{\eta}_2)^{3/2}}$. Furthermore, under additional primitive conditions, we also show that the conditions of Theorem (ref) hold with precise non-asymptotic probability bounds characterized in the proof.
Theorem SA-1 characterizes precisely the probability guarantees for the in-sample prediction---that is, the precise values of $\alpha_1$ and $\pi_1$ in Lemma (ref) obtained via Theorems (ref) and (ref). The conditions imposed are primitive (e.g., moment bounds and rank conditions), with perhaps the exception of conditions (SA-1.iii) and (SA-1.v) in Theorem SA-1 in the supplemental appendix. Specifically, Condition (SA-1.iii) requires $\varrho=\varpi_\delta^{\star}/\sqrt{T_0}$ and $\mathbb{P}(\min\{|w_{0,j}|:w_{0,j}\neq 0\}\geq\varrho)\geq 1-\pi_{w}^{\star}$, for non-negative constants $\varpi_\delta^{\star}$ and $\pi_{w}^{\star}$, which is also primitive insofar it relates to the separation from zero of the non-zero (possibly random) coefficients $\mathbf{w}_0$ entering the best linear approximation (ref), which is a standard (sparsity-type) assumption in the literature of constrained least squares estimation. On the other hand, Condition (SA-1.v) requires $\mathbb{P}[\mathbb{P}(\max_{1\leq t\leq T_0}|\widehat{\mathbb{E}}[u_t|\mathbf{b}_t]-\mathbb{E}[u_t|\mathbf{b}_t]|\leq \varpi_u^{\star}|\mathscr{H})\geq 1-\epsilon_u^{\star}]\geq 1-\pi_u^{\star}$, for non-negative constants $\varpi_u^{\star}$, $\epsilon_u^{\star}$ and $\pi_u^{\star}$, which is purposely not as primitive (but still easily interpretable) because it is meant to cover many different approximation approaches for $\mathbb{E}[u_t|\mathbf{b}_t]$. In practice, researchers may assume $\mathbb{E}[u_t|\mathbf{b}_t]=0$ or, alternatively, employ flexible-parametric/nonparametric approaches to form the estimator $\widehat{\mathbb{E}}[u_t|\mathbf{b}_t]$. Since the latter approaches are setting-specific and technically well-understood, we chose to present our results using the generic condition (SA-1.v) rather than providing primitive conditions for a specific example of $\widehat{\mathbb{E}}[u_t|\mathbf{b}_t]$.
The second example is the multi-equation setup introduced in Section (ref), where we incorporate pre-intervention covariates in the construction of the SC weights and allow for stationary weakly dependent time series data. See Kilian-Lutkepohl_2017_Book and references therein for an introduction to time series analysis. We let $M=2$ (two features) and $K=0$ (no additional controls) for simplicity, which gives the working model
$t=1, \cdots, T_0$. The first equation could naturally correspond to pre-intervention outcomes as in the previous example, i.e., $a_{t,1}:=Y_{1t}(0)$ and $\mathbf{b}_{t,1} :=(Y_{2t}(0), Y_{3t}(0), \ldots, Y_{(N+1)t}(0))'$, while the second equation could correspond to some other covariate (such as population density in the Basque terrorism application) also used to construct $\widehat{\mathbf{w}}$ in (ref). Let $\mathbf{b}_{t,l}=(b_{1t,l}, \cdots, b_{Jt,l})'$, for $l=1,2$. To provide interpretable primitive conditions, we also assume $\mathbf{u}_t=(u_{t,1},u_{t,2})'$ and $\mathbf{b}_t=(\mathbf{b}_{t,1}',\mathbf{b}_{t,2}')'$ follow independent first-order stationary autoregressive (AR) processes:
where $\bm{\zeta}_{t,u}$ and $\bm{\zeta}_{t,b}$ are i.i.d. over $t$, independent of each other, and $\operatorname*{diag}(\cdot)$ denotes a diagonal matrix with the function arguments as the corresponding diagonal elements. Let $\mathbf{D}=T_0^{1/2}\mathbf{I}_d$, and note that $\mathbf{U}=(u_{1,1}, \cdots, u_{T_0,1}, u_{1,2}, \cdots, u_{T_0, 2})'$ in this case. A natural, generic variance estimator is \[\widehat{\boldsymbol{\Sigma}}=\frac{1}{T_0}\mathbf{Z}'\widehat{\mathbb{V}}[\mathbf{U}|\mathscr{H}]\mathbf{Z},\] where $\widehat{\mathbb{V}}[\mathbf{U}|\mathscr{H}]$ is an estimate of $\mathbb{V}[\mathbf{U}|\mathscr{H}]$. In this example, $\boldsymbol{\Sigma}$ corresponds to the (conditional) long-run variance, and naturally $\widehat{\boldsymbol{\Sigma}}$ can be chosen to be any standard estimator thereof.
Theorem SA-2 in the supplemental appendix gives primitive conditions that verify the high-level conditions of Theorem (ref) in the appendix, and the high-level conditions of Theorem (ref) for implementation. Note that because of the time dependence in this example, the primitive conditions are for Theorem (ref) instead of Theorem (ref). In particular, we show that under standard conditions guaranteeing $\beta$-mixing and moment and rank conditions (similar to those imposed in the previous example), the conditions of Theorem (ref) hold with $\pi_{\gamma}=\mathfrak{C}_\pi T_0^{-\mathfrak{c}_\pi}$ and $\epsilon_{\gamma}=\mathfrak{C}_\epsilon T_0^{-\mathfrak{c}_\epsilon}$ for non-negative constants $\mathfrak{C}_\pi$ and $\mathfrak{C}_\epsilon$, and some positive constants $\mathfrak{c}_\pi$ and $\mathfrak{c}_\epsilon$, which are characterized precisely in the supplemental appendix. Theorem (ref) is also verified using the primitive conditions imposed in Theorem SA-2 in the supplemental appendix, and the associated non-asymptotic constants are characterized in its proof.
As in the previous example, Theorem SA-2 in the supplemental appendix illustrates the kind of primitive conditions needed to quantify in-sample uncertainty using our proposed methods. In this example, we accommodate multiple covariates (equations) in the construction of the SC weights and also allow for AR(1) dependent (stationary) time series data. The only intentionally high-level condition imposed is (SA-2.v), $\mathbb{P}(\mathbb{P}(\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\|\leq \epsilon_{\Sigma,1}^{\star}|\mathscr{H})\geq 1-\epsilon_{\Sigma,2}^{\star})\geq 1-\pi_{\Sigma}^{\star}$ for non-negative constants $\epsilon_{\Sigma,1}^{\star}$, $\epsilon_{\Sigma,2}^{\star}$ and $\pi_{\Sigma}^{\star}$, which requires a concentration probability bound for the long-run variance estimator $\widehat{\boldsymbol{\Sigma}}$ used to approximate the quantiles of the conditional (on $\mathscr{H}$) distribution of the bounding random variables $\varsigma_{\mathtt{L}}$ and $\varsigma_{\mathtt{U}}$ via simulations (Theorem (ref)). This condition is not difficult to verify for specific examples.
Our third and final example illustrates how non-stationary data can also be handled within our framework. See Tanaka_2017_Book and references therein for an introduction to non-stationary time series analysis. Suppose that for each $1\leq l\leq M$, $\{a_{t,l}\}_{t=1}^T$, $\{b_{1t,l}\}_{t=1}^T, \cdots, \{b_{Jt,l}\}_{t=1}^T$ are $I(1)$ processes, and $\{c_{1t,l}\}_{t=1}^T, \cdots, \{c_{Kt,l}\}_{t=1}^T$ and $\{u_{t,l}\}_{t=1}^{T}$ are $I(0)$ processes. Therefore, $\mathbf{A}$ and $\mathbf{B}$ form a cointegrated system. For simplicity, consider the following example: for each $l=1, \cdots, M$ and $j=1, \cdots, J$,
where $u_{t,l}$ and $v_{jt,l}$ are stationary unobserved disturbances. In this scenario, $(1, -\mathbf{w}_0')'$ plays the role of a cointegrating vector such that the linear combination of $\mathbf{A}$ and $\mathbf{B}$ is stationary. The normalizing matrix $\mathbf{D}=\mathrm{diag}\{T_0, \cdots, T_0, \sqrt{T_0}, \cdots, \sqrt{T_0}\}$, where the first $J$ elements are $T_0$ and the remaining ones are $\sqrt{T_0}$. Let $\check{\mathbf{Z}}_t=(\check{\mathbf{z}}_{t,1}, \cdots, \check{\mathbf{z}}_{t,M})$ where $\check{\mathbf{z}}_{t,l}$ is the $((l-1)T_0+t)$th column of $\operatorname*{diag}\{T_0^{-1/2}\mathbf{I}_J, \mathbf{I}_{KM}\}\mathbf{Z}'$, for $l=1, \cdots, M$. Recall that $\mathbf{u}_t=(u_{t,1}, \cdots, u_{t,M})'$. Write $\mathbf{v}_{t,l}=(v_{1t,l}, \cdots, v_{Jt,l})'$, $\mathbf{v}_t=(\mathbf{v}_{t,1}', \cdots, \mathbf{v}_{t, M}')'$, and $\mathbf{c}_{t,l}=(c_{1t,l}, \cdots, c_{kt,l})'$. We allow some elements in $\mathbf{v}_t$ to be used in $\{\mathbf{c}_{t,l}\}_{l=1}^M$. Let $\mathbf{q}_t$ collect all distinct variables in $\mathbf{u}_t$, $\mathbf{v}_t$, $\mathbf{c}_{t,1}$, $\cdots$, $\mathbf{c}_{t,M}$. As in the previous example, a generic variance estimator is \[\widehat{\boldsymbol{\Sigma}}=\frac{1}{T_0}\sum_{t=1}^{T_0}\check{\mathbf{Z}}_t\widehat{\mathbb{V}}[\mathbf{u}_t|\mathscr{H}]\check{\mathbf{Z}}_t',\] where $\widehat{\mathbb{V}}[\mathbf{u}_t|\mathscr{H}]$ is an estimate of $\mathbb{V}[\mathbf{u}_t|\mathscr{H}]$.
Theorem SA-3 in the supplemental appendix gives more primitive conditions and verifies the high-level conditions of Theorems (ref) and (ref) in the cointegration scenario. More precisely, it provides conditions so that Theorem (ref) holds with $\pi_{\gamma}=\mathfrak{C}_{\pi,1}T_0^{-\psi\nu}+\mathfrak{C}_{\pi,2}T_0^{-1}+\pi_{Q,1}+\pi_{Q,2}$ and $\epsilon_{\gamma}=\mathfrak{C}_\epsilon (\log T_0)^{\frac{3}{2}(1+\mathfrak{c}_Q)}T_0^{-1/2}$ for non-negative constants $(\psi,\nu,\pi_{Q,1},\pi_{Q,2},\mathfrak{c}_Q)$ specified in the assumptions of the theorem and non-negative constants $(\mathfrak{C}_{\pi,1},\mathfrak{C}_{\pi,2},\mathfrak{C}_\epsilon)$ characterized in the proof. Similarly, Theorem (ref) is also verified under more primitive conditions, including a higher-level condition of the form $\mathbb{P}(\mathbb{P}(\|\widehat{\boldsymbol{\Sigma}}-\boldsymbol{\Sigma}\|\leq \epsilon_{\Sigma,1}^{\star}|\mathscr{H})\geq 1-\epsilon_{\Sigma,2}^{\star})\geq 1-\pi_{\Sigma}^{\star}$ for non-negative constants $\epsilon_{\Sigma,1}^{\star}$, $\epsilon_{\Sigma,2}^{\star}$ and $\pi_{\Sigma}^{\star}$, as in the previous examples.
When $\mathbf{C}$ is excluded, $\widehat{\mathbf{w}}$ is a least squares estimator of the cointegrating vector, which is typically biased due to the potential correlation between $\mathbf{v}_t$ and $\mathbf{u}_t$. In Theorem SA-3 in the supplemental appendix, we include $\mathbf{C}$ and allow it to include contemporary $\mathbf{v}_t$ to correct this bias. More generally, one may augment the regression with $\mathbf{v}_t$ and its leads and lags, which is termed dynamic OLS in the time series literature. The results for this general case may be established using a similar strategy.
The unobserved random variable $e_T$ in (ref) is a single error term in period $T$, which we interpret as the error from out-of-sample prediction, conditional on $\mathscr{H}= \sigma(\mathbf{B}, \mathbf{C}, \mathbf{x}_T, \mathbf{g}_T)$. Naturally, in order to set appropriate $M_{2,\mathtt{L}}$ and $M_{2,\mathtt{U}}$ in Lemma (ref), it is necessary to determine certain features of the conditional distribution $\mathbb{P}[e_T \leq \cdot | \mathscr{H}]$. In turn, determining those features would require strong distributional assumptions between pre-treatment and post-treatment periods, or perhaps across units. In this section we propose principled but agnostic approaches to quantify the uncertainty introduced by the post-treatment unobserved shock $e_T$. Since formalizing the validity of our methods requires strong assumptions, in this paper we recommend a generic sensitivity analysis to incorporate out-of-sample uncertainty to the prediction intervals. In particular, we propose employing three distinct methods for quantifying the uncertainty introduced by $e_T$ as a starting point, and then assessing more generally whether the additional uncertainty would render the prediction intervals large enough to eliminate any statistically significant treatment effect.
Our starting point is a non-asymptotic probability bound on $e_T$ via concentration inequalities. Such textbook results can be found in, for example, Vershynin_2018_Book and Wainwright_2019_Book. We rely on the following lemma, which provides the desired bounds for $e_T$ under different moment-like conditions.
This lemma gives (possibly crude) bounds on the necessary features of the conditional distribution of $e_T$ given $\mathscr{H}$. Lemma (ref)(G) corresponds to a sub-Gaussian tail assumption, while Lemma (ref)(P) exploits only a polynomial bound on moments of $e_t|\mathscr{H}$. In both cases, the only unknowns are the “center” and “scale” of the distribution: $\mathbb{E}[e_T|\mathscr{H}]$ and $\sigma_{\mathscr{H}}^2$ (or higher-moments), respectively. These unknown features can be estimated or tabulated based on (i) model assumptions and (ii) observed pre-treatment data, at least as an initial step towards a sensitivity analysis.
For practical purposes, we first outline three alternative strategies to assess the uncertainty coming from $e_T$, starting with Lemma (ref) and progressively adding more restrictions. After introducing these approaches, we turn to discussing how they can be used as an initial step towards a principled sensitivity analysis for uncertainty quantification of the synthetic control estimator. Section (ref) illustrates this idea using simulated data and empirical applications.
While the three approaches above are simple and intuitive, they are potentially unsatisfactory because their validity would require arguably strong assumptions on the underlying data generating process linking the pre-treatment and post-treatment data. Such assumptions, however, are difficult to avoid because the ultimate goal is to learn about uncertainty introduced by an unobserved random variable after the treatment began (i.e., $e_T|\mathscr{H}$ for $T>T_0$). Without additional data availability or specific modelling assumptions allowing for transferring information from the pre-treatment period into the post-treatment period, it is difficult to formally set $M_{2,\mathtt{L}}$ and $M_{2,\mathtt{U}}$ in Lemma (ref).
Nevertheless, it is possible to approach out-of-sample uncertainty quantification as a principled sensitivity analysis, using the methods above as a starting point. Given the formal and detailed in-sample uncertainty quantification developed in the previous section, it is natural to progressively enlarge the final prediction intervals by adding additional out-of-sample uncertainty to then ask the question: how large does the additional out-of-sample uncertainty contribution coming from $e_T|\mathscr{H}$ need to be in order to render the treatment effect $\tau_t$ in (ref) statistically insignificant? Using the approaches above, or similar ones, it is possible to construct natural initial benchmarks. For instance, the variability displayed by the pre-treatment outcomes or synthetic control residuals can help guide the level of “reasonable” out-of-sample uncertainty. Alternatively, in specific applications, natural levels of uncertainty for the outcomes of interests could be available, and hence used to tabulate the additional out-of-sample uncertainty. In Section (ref) we further discuss and illustrate this idea numerically.
We revisit the three examples considered in Section (ref) and illustrate how the implementation of the three approaches outlined earlier may accommodate different assumptions on the data generating process. We discuss the outcomes-only case in more detail, which suffices to showcase our basic strategy of out-of-sample uncertainty quantification. For the other two examples, we briefly explain some important conceptual and implementational issues. As mentioned above, these methods rely on strong assumptions and should be viewed as a starting point of a general sensitivity analysis.
Recall that the data is assumed to be i.i.d. over $1\leq t\leq T$ in this case. The conditional distribution of $e_T$ given $\mathscr{H}$ then reduces to that given the contemporary covariates $\mathbf{b}_T$ only. Also, we set $\mathbf{x}_T=\mathbf{b}_T$ and $\mathbf{g}_T=(1, \cdots, 1)'$. Then, the out-of-sample error $e_T$ is equivalent to the pseudo-true residual $u_T$. By stationarity of the data, the information about the conditional distribution of $u_T$ can be learned using the pre-intervention residuals. These substantial simplifications facilitate the implementation of the proposed methods for quantifying the out-of-sample uncertainty.
When data is weakly dependent and multiple features are used in the construction of SC weights, the implementation of the three approaches is similar to that in the outcomes-only case, but two outstanding issues need to be addressed. First, as described in Section (ref), the SC weights are obtained by matching on two pre-intervention features $\{a_{t,1}\}_{t=1}^{T_0}$ and $\{a_{t,2}\}_{t=1}^{T_0}$, while in most SC applications, the final counterfactual prediction is constructed by setting $\mathbf{x}_T=\mathbf{b}_T$ (and $\mathbf{g}_T=\bm{0}$ in this case). Conceptually, the out-of-sample error $e_T=Y_{1T}(0)-\mathbf{b}_T'\mathbf{w}_0$ may or may not correspond to the pseudo-true residual $\mathbf{u}_t$ prior to the treatment. For example, if the pre-treatment outcomes are used in the first equation, then by construction, $e_t$ in this scenario is equivalent to $u_{t,1}$ for $t=1, \cdots, T$. In the pre-intervention period, the residuals $\{\widehat{u}_{t,1}\}_{t=1}^{T_0}$ from the SC fitting play the role of proxies for $\{u_{t,1}\}_{t=1}^{T_0}$. In contrast, if pre-treatment outcomes are not used in any of the two equations, then $e_t$ is generally not the same as $u_{t,1}$ or $u_{t,2}$. Nevertheless, we can manually construct $\widehat{e}_t=Y_{1t}(0)-\mathbf{b}_t'\widehat{\mathbf{w}}$ as a proxy for $e_t$ in the pre-intervention period.
Second, the dependence of $e_T$ on $\mathscr{H}$ should be appropriately characterized. Consider a simple scenario where pre-treatment outcomes are used in the construction of SC weights so that $e_t=u_{t,1}$. By our assumptions on $\bm\zeta_{t,u}$ and $\bm\zeta_{t,b}$, $\{u_{t,1}\}_{t=1}^{T}$ is independent of $\{\mathbf{b}_t\}_{t=1}^T$. If we further assume the two components of $\bm\zeta_{t,u}$ are independent of each other, the conditional distribution of $u_{T,1}$ given $\mathscr{H}$ reduces to its unconditional distribution. Therefore, to implement the three approaches, one only needs to estimate the unconditional mean, variance or quantiles using the residuals $\{\widehat{u}_{t,1}\}_{t=1}^{T_0}$. In practice, however, the independence between $u_{T,1}$ and $\mathscr{H}$ may be unrealistic. Assuming an appropriate weak dependence structure, we can still characterize or approximate the conditional mean, variance or quantiles of $u_{t,1}$ by functions of $\mathbf{b}_t$ and lags thereof, which could be estimated by parametric or nonparametric regressions using the pre-intervention data.
As in the second example, we first determine the pre-intervention analogue to the out-of-sample error $e_T$. For instance, we let $a_{t,1}=Y_{1t}(0)$ and $b_{jt,1}=Y_{(j+1)t}(0)$ for $j=1, \cdots, N$. In practice, the final counterfactual prediction is often constructed by setting $\mathbf{x}_T=(Y_{2t}(0), \cdots, Y_{(N+1)t}(0))'$ and $\mathbf{g}_T=\bm{0}$, i.e., no additional control variables are used in the out-of-sample prediction. Then, we have $e_t=u_{t,1}+\sum_{k=1}^Kc_{kt,1}r_{0,k,1}$ for $t=1, \cdots, T$. In view of the assumptions imposed in Theorem SA-3 in the supplemental appendix, the conditional distribution of $e_t$ given $\mathscr{H}$ reduces to that given the contemporary variables $\{\mathbf{v}_t,\mathbf{c}_{t,1}, \cdots, \mathbf{c}_{t,M}\}$. As in the previous examples, we can estimate its conditional mean, variance or quantiles using various parametric or nonparametric methods.
The assumption that $\{\mathbf{q}_t\}_{t=1}^T$ is i.i.d. in Theorem SA-3 may be too strong. In practice, as mentioned previously, we may want to augment the regression of the residual $\widehat{e}_t$ by lags (and leads) of $\mathbf{v}_t$ and $\{\mathbf{c}_{t,l}\}_{l=1}^M$ and transformations thereof to take into account potential time series dependence.
We illustrate the performance of the proposed prediction intervals with a Monte Carlo experiment and two empirical examples. To implement the methods described in Section (ref) and (ref), we take a simple “plug-in” estimator $\widehat{\boldsymbol{\Sigma}}$ of the long-run variance $\boldsymbol{\Sigma}$ and employ parametric polynomial regressions to estimate the conditional mean, variance and quantiles of $e_T$ given $\mathscr{H}$ whenever needed. In addition, the choice of the tuning parameter $\varrho$ can be based on a bound implied by optimization. Specifically, since $\widehat{\boldsymbol{\delta}}$ must satisfy the basic inequality specified in the definition of $\mathcal{M}_\xi$ (see Lemma (ref)), we can construct a threshold $\varrho$ for $\widehat{\boldsymbol{\beta}}$ based on some estimates of the variance of $\widehat{\boldsymbol{\gamma}}$ and the eigenvalues of $\widehat{\mathbf{Q}}$. More details are discussed below, and we also provide complete replication codes. Last but not least, in view of the small sample size in many SC applications, these practical choices play the role of a reasonable starting point for a principled sensitivity analysis, as illustrated in this section.
We conduct a Monte Carlo investigation of the finite sample performance of our proposed methods. We consider the outcomes-only case where $M=1$, $J=N$, $K=0$, and only the outcome variable is used. Then, $\mathbf{A}_1=(Y_{11}, Y_{12}, \cdots, Y_{1T_0})'$, $\mathbf{b}_{j,1}=(Y_{(j+1)1}, Y_{(j+1)2}, \cdots, Y_{(j+1)T_0})'$. We set $T_0=100$, $T_1=1$, and $N=10$. We consider three data generating processes for $b_{jt}:=b_{jt,1}$: $b_{jt}=\rho b_{j(t-1)}+v_{jt}$ with $\rho\in\{0,0.5,1\}$, and where $\mathbf{v}_t=(v_{1t}, \cdots, v_{Nt})'\thicksim i.i.d. \;\mathsf{N}(0, \mathbf{I}_N)$.
To examine the conditional coverage, we first generate a sample of $\{b_{jt}:1\leq t\leq T_0+T_1, 1\leq j\leq N\}$ using one of the three models. We set $5$ evaluation points in the post-treatment period: $\tilde{b}_{1(T_0+1)}:=b_{1(T_0+1)}+\mathtt{c}\cdot \mathsf{sd}(b_{1t})$ where $\mathtt{c}\in\{-1, -0.5, 0, 0.5, 1\}$ and $\mathsf{sd}(b_{1t})$ is the sample standard deviation of $\{b_{1t}\}_{t=1}^{T_0}$. In other words, we construct $5$ designs by varying the value of the first conditioning variable in the last period. Taking each of them as given, we generate the treated unit $a_t:=Y_{1t}=\mathbf{b}_t'\mathbf{w}+u_t$ by randomly drawing the error $u_t\thicksim i.i.d.\;\mathsf{N}(0,0.5)$ independent of $\{\mathbf{b}_t\}_{t=1}^T$. We set $\mathbf{w}=(0.3, 0.4, 0.3, 0, \cdots, 0)'$. By construction, $\mathbf{w}$ is exactly equivalent to the pseudo-true weight $\mathbf{w}_0$ defined in (ref) and satisfies the positivity and sum-to-one constraints in the standard SC. We consider $5,000$ simulated datasets. Note that the design $\{\mathbf{b}_t\}_{t=1}^{T_0}$ is fixed throughout the simulation study, and we only draw a new realization of the error term $u_t$ at each repetition.
For comparison we also investigate the unconditional coverage of the prediction intervals. The models for $a_t$ and $\mathbf{b}_t$ are the same as described above, but in this case the design $\{\mathbf{b}_t\}_{t=1}^T$ is not fixed across the simulation study. Instead, at each repetition we randomly draw $\{(a_t, \mathbf{b}_t)\}_{t=1}^{T_0+1}$. Therefore, the coverage probability obtained in this alternative exercise is unconditional, i.e., not conditional on a fixed design.
In the above construction, the error term $u_t$ is independent of the design and has mean zero. To check the performance of the proposed method in models with misspecification error ($\mathbb{E}[u_t|\mathscr{H}]\neq 0$), we also consider several other data generating processes: for models with $\rho=0$ and $\rho=0.5$, we generate $u_t=0.2b_{1t}+\zeta_{t}$, while for the model with $\rho=1$, $u_t=0.9(b_{1t}-b_{1(t-1)})+\zeta_t$, where $\zeta_t\thicksim i.i.d. \;\mathsf{N}(0, 0.5)$. Notably, in these misspecified models, the weight $\mathbf{w}$ defined previously is no longer equivalent to the pseudo-true weight $\mathbf{w}_0$.
We focus on several versions of our proposed prediction intervals for the counterfactual outcome $Y_{1T}(0)$ of the treated unit with $90\%$ nominal (conditional) coverage probability: “M1” denotes the prediction interval based on assuming a Gaussian bound and applying the concentration inequality in Lemma (ref)(G) (“approach 1”); “M1-S” is the same as “M1” except that we increase the (estimated) conditional standard deviation $\widehat{\sigma}_{\mathscr{H}}$ by a factor of $2$; “M2” denotes the prediction interval based on the location-scale model (“approach 2”); and “M3” denotes the prediction interval based on linear quantile regression (“approach 3”). For comparison, we also include the prediction interval “CONF”, which is based on the conformal method developed in Chernozhukov-Wuthrich-Zhu_2021_JASA. In the supplemental appendix, we also report the performance of the prediction interval based on the cross-sectional permutation method proposed in Abadie-Diamond-Hainmueller_2010_JASA.
As mentioned previously, we choose the tuning parameter $\varrho$ based on the basic optimization inequality $\varrho=\widehat{\sigma}_u(\log T_0)^{\mathtt{c}}/(\min_{1\leq j\leq J}\widehat{\sigma}_{b_j}T_0^{1/2})$, where $\widehat{\sigma}_u$ is the estimated (unconditional) standard deviation of $u_t$, $\widehat{\sigma}_{b_j}$ is the estimated (unconditional) second moment of $b_{jt}$, and $\mathtt{c}=1$ if $\rho=1$ and $\mathtt{c}=0.5$ if $\rho=0$ or $0.5$. This strategy accommodates the simple data generating processes considered in simulations, and it can be tailored to better suit the assumptions in more complex statistcal models as well. In addition, we use polynomial regression to estimate various features of the conditional distribution of $u_t$ given $\mathscr{H}$. To avoid overfitting, we only use a subset of the $J$ control outcomes which have non-zero weights in $\widehat{\mathbf{w}}^{\star}$. Regarding the long-run variance $\boldsymbol{\Sigma}$, we take a simple plug-in estimator $\widehat{\boldsymbol{\Sigma}}=\mathbf{D}^{-1}\mathbf{Z}'\operatorname*{diag}\{\tilde{u}_1^2, \cdots, \tilde{u}^2_{T_0}\}\mathbf{Z}\mathbf{D}^{-1}$, where $\tilde{u}_t=\widehat{u}_t-\widehat{\mathbb{E}}[u_t|\mathscr{H}]$ and $\widehat{\mathbb{E}}[u_t|\mathscr{H}]$ is the estimate of the conditional mean of $u_t$ given $\mathscr{H}$.
Panels A and B of Table (ref) summarize the results for models with and without misspecification error, respectively, where we use linear regression methods to estimate the conditional mean, variance or quantiles whenever needed. The proposed prediction intervals exhibit good coverage properties throughout different data generating processes and evaluation points, though they are conservative in several cases. In contrast, the actual coverage probability of conformal prediction intervals developed in Chernozhukov-Wuthrich-Zhu_2021_JASA is lower than the target nominal level in general. See Section SA-3 of the supplement for additional simulation evidence.
We showcase our methods by reanalyzing two empirical examples from the synthetic control literature. The first example concerns the effect of California's tobacco control program, known as Proposition 99, on per capita cigarette sales Abadie-Diamond-Hainmueller_2010_JASA. The second example corresponds to the economic impact of 1990 German reunification on West Germany Abadie_2021_JEL. To conserve space, the results for the second example are reported in Section SA-4 of the supplement.
The outcome variable of interest is per capita cigarette sales in California, which is arguably non-stationary. We consider both the raw data and the first-differenced data, corresponding to the analysis of levels and growth rates of per capita sales, respectively. In each scenario, we construct (1) the synthetic control prediction $\mathbf{x}_T'\widehat{\mathbf{w}}$ as commonly done in the literature; (2) the prediction interval for the (conditionally non-random) “synthetic control component” $\mathbf{x}_T'\mathbf{w}_0$ only; and (3) three distinct prediction intervals for the counterfactual $Y_{1T}(0)$. More specifically, the prediction intervals for $Y_{1T}(0)$ are implemented using the three methods outlined in Section (ref): (i) conditional subgaussian bound using Lemma (ref)(G), labeled as approach 1; (ii) conditional bound based on location-scale model, labeled as approach 2; (iii) conditional bound using conditional quantile regression of residuals, labeled as approach 3.
The three methods used to quantify the out-of-sample uncertainty can be viewed as particular instances of a more general sensitivity analysis. In other words, varying the additional uncertainty contribution of $e_T$ in a principled way, researchers can better understand its impact on the construction of the prediction intervals. We also illustrate this approach (“sensitivity analysis”): focusing on approach 1 for concreteness, we rely on Gaussian bounds in Lemma (ref)(G) to assess how the prediction intervals behave as the variance of $e_T$ varies.
Accordingly, we present six plots: the SC prediction $\mathbf{x}_T'\widehat{\mathbf{w}}$, the prediction interval (PI) only for the synthetic unit $\mathbf{x}_T'\mathbf{w}_0$ with at least $95\%$ nominal coverage probability, the three different constructions of PIs for the counterfactual $Y_{1T}(0)$ with at least $90\%$ nominal coverage probability, and a sensitivity analysis for one chosen post-treatment period. Because the size of the donor pool is larger than the number of available pre-treatment periods, our procedure may rely on loose bounds for Gaussian approximation errors in high-dimensional settings.
We first consider the raw data of per capita cigarette sales. Figure (ref)(a) shows the trajectory of per capita sales of the synthetic California (dashed blue) and the actual California (solid black). After 1988, the synthetic California series is above the observed one, suggesting a negative shock of Proposition 99 on cigarette sales in California. Figure (ref)(b) adds a $95\%$ conservative prediction interval for the synthetic control component of California that takes into account the in-sample uncertainty due to the estimated SC weights. We add the uncertainty associated with $e_T$ in Figures (ref)(c)-(e). The observed sequence is generally below the prediction intervals for the counterfactual outcome of California, suggesting statistically significant effects of Proposition 99. Figure (ref)(f) shows the sensitivity analysis of the effect in 1989. The result is robust: the corresponding PIs are well separated from the observed outcome of California if we vary the estimated (conditional) standard deviation of $e_T$ in a relatively wide range.
We also analyze the (log) growth rate of per capita cigarette sales. The result is reported in Figure (ref). We can see that the observed growth rate during the post-treatment period is generally lower than the SC prediction, but throughout the three constructions of PIs for the counterfactual outcome, the observed series is within the PIs for most post-treatment periods except in 1989. These empirical findings suggest a statistical significant effect of the tobacco control program on the growth rate of per capita cigarette sales only in year 1989. The sensitivity analysis in Figure (ref)(f) shows that the significance of the effect in 1989 is robust.
The synthetic control method is part of the standard program evaluation toolkit. Despite its popularity, many important methodological and theoretical developments remain outstanding. We focus on quantifying the uncertainty of the SC method in predicting the main quantity of interest, $\tau_T=Y_{1T}(1)-Y_{1T}(0)$, in the standard SC framework. This quantity is the difference between the observed outcome of the treated unit in a post-treatment period $T$, and the outcome that the treated unit would have had in the same period in the absence of treatment. Because we view $\tau_T$ as a random variable and there is a single treated unit, we propose conditional prediction intervals that offer finite-sample probability guarantees regarding the realization of the counterfactual treated outcome. Our approach takes the SC constrained least squares optimization approach as the starting point. We model the counterfactual of the treated unit in period $T$ as the weighted sum of the untreated units' features at $T$ (with weights estimated with pre-treatment data), and an error term. This decomposition highlights two sources of uncertainty, one from the in-sample estimation of the SC weights in the pre-treatment period, and the other from the post-treatment error that arises due to the unavoidable out-of-sample prediction involved in the SC method, which may include potential misspecification errors from the SC weights. Using finite-sample concentration bounds, we derive prediction intervals that incorporate both sources of uncertainty. Because the uncertainty stemming from the out-of-sample post-treatment error term is hard to handle (specially under general misspecification), we recommend combining the prediction interval for the SC outcome with a principled sensitivity analysis for the post-treatment error. Our empirical illustrations show that our methods perform well using both simulated and real data. A general-purpose software package is underway \citep*{Cattaneo-Feng-Palomba-Titiunik_2021_scpi}.
\onehalfspacing