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.
90,651 characters · 18 sections · 13 citation commands
scpi: Uncertainty Quantification for Synthetic Control Methods
Keywords: program evaluation, causal inference, synthetic controls, prediction intervals, non-asymptotic inference.
\thispagestyle{empty}
\doublespacing \setcounter{page}{1}
\setcounter{page}{0}\thispagestyle{empty}
\pagestyle{plain}
The synthetic control method was introduced by \citet*{Abadie-Gardeazabal_2003_AER}, and since then it has become a popular approach for program evaluation and causal inference in observational studies. It offers a way to study the effect of an intervention (e.g., treatments at the level of aggregate units, such as cities, states, or countries) by constructing weighted averages of untreated units to approximate the counterfactual outcome that the treated unit(s) would have experienced in the absence of the intervention. While originally developed for the special case of a single treated unit and a few control units over a short time span, the methodology has been extended in recent years to a variety of other settings with longitudinal data. See \citet*{Abadie_2021_JEL} for a review on synthetic control methods, and \citet*{Abadie-Cattaneo_2018_ARE} for a review on general methods for program evaluation.
Most methodological developments in the synthetic control literature have focused on either expanding the causal framework or developing new implementations for prediction/point estimation. Examples of the former include disaggregated data settings \citep*{Abadie-LHour_2021_JASA} and staggered treatment adoption \citep*{BenMichael-Feller-Rothstein_2022_JRSSB}, while examples of the latter include employing different constrained estimation methods (see Table (ref) below for references). Conceptually, implementation of the synthetic control method involves two main steps: first, treated units are “matched” to control units using only their pre-intervention data via (often constrained) regression methods and, second, prediction of the counterfactual outcomes of the treated units are obtained by combining the pre-intervention “matching” weights with the post-intervention data of the control units. As a result, the synthetic control approach offers a prediction or point estimator of the (causal) treatment effect for the treated unit(s) after the intervention was deployed.
Compared to prediction or estimation, considerably less effort has been devoted to develop principled uncertainty quantification for synthetic control methods. The most popular approach in practice is to employ design-based permutation methods taking the potential outcome variables as non-random \citep*{Abadie-Diamond-Hainmueller_2010_JASA}. Other approaches include methods based on large-sample approximations for disaggregated data under correctly specified factor-type models \citep*{Li_2020_JASA}, time-series permutation-based inference \citep*{Chernozhukov-Wuthrich-Zhu_2021_JASA}, large-sample approximations for high-dimensional penalization methods \citep*{Masini-Medeiros_2021_JASA}, and cross-sectional permutation-based inference in semiparametric duration-type settings Shaikh-Toulis_2021_JASA. A conceptually distinct approach to uncertainty quantification is proposed by \citet*{Cattaneo-Feng-Titiunik_2021_JASA} and \citet*{Cattaneo-Feng-Palomba-Titiunik_2022_wp}, who take the potential outcome variables as random and develop prediction intervals for the imputed (counterfactual) outcome of the treated unit(s) in the post-intervention period employing finite-sample probability concentration methods.
This article introduces the software package scpi for prediction and inference using synthetic control methods, implemented in Python, R, and Stata. For prediction or point estimation of treatment effects, the package offers an array of possibly penalized approaches leveraging the latest conic optimization methods \citep*{Domahidi-Chu-Boyd_2013_ECOS, Fu-Narasimhan-Boyd_2020_JSS}; see also Boyd_2004_BOOK for an introduction. For uncertainty quantification, the package focuses on the aforementioned prediction interval methods under random potential outcomes. The rest of the article focuses on the $\mathtt{R}$ implementation of the software, but we briefly illustrate analogous functionalities for Python in Appendix (ref), and for Stata in Appendix (ref).
The R package scpi includes the following six functions:
The objects returned by scest() and scpi() support the methods print() and summary(). In typical applications, the user will first prepare the data using the function scdata() or scdataMulti(), and then produce predictions/point estimators for treatment effects with uncertainty quantification using the function \texttt{scpi()}. The function \texttt{scest()} is useful in cases where only predictions/point estimators are of interest. Numerical illustrations are given in Section (ref).
There are many R, Python, and Stata packages available for prediction/point estimation and inference using synthetic control methods; Table (ref) compares them to the package scpi. As shown in the table, scpi is the first package to offer uncertainty quantification using prediction intervals with random potential outcomes for a wide range of different synthetic control predictors. The package is also one of the first to handle multiple treated units and staggered treatment adoption, offering a wider array of options in terms of predictors and inference methods when compared with the other packages currently available. Furthermore, the package includes misspecification-robust methods, employs the latest optimization packages for conic programs available, and offers automatic parallelization in execution whenever multi-core processors are present, leading to significant improvements in numerical stability and computational speed. Finally, scpi is the only package available in \texttt{Python}, \texttt{R}, and \texttt{Stata}, which gives full portability across multiple statistical software and programming languages, and also the only package employing directly native conic optimization via the \texttt{ECOS} solver (see Table (ref) for details).
The rest of the article is organized as follows. Section (ref) introduces the canonical synthetic control setup, and also briefly discusses extensions to multiple treated units with possibly staggered treatment adoption. Section (ref) gives a brief introduction to the theory and methodology underlying the point estimation/prediction for synthetic control methods, discussing implementation details. Section (ref) gives a brief introduction to the theory and methodology underlying the uncertainty quantification via prediction intervals for synthetic control methods, and also discusses the corresponding issues of implementation. Section (ref) showcases some of the functionalities of the package using a real-world dataset, and Section (ref) concludes. The appendices illustrate the Python (Appendix (ref)) and Stata (Appendix (ref)) implementations of scpi. Detailed instructions for installation, script files to replicate the analyses, links to software repositories, and other companion information can be found in the package's website, \url{https://nppackages.github.io/scpi/}.
We first consider the canonical synthetic control framework with a single treated unit. The researcher observes $J+1$ units for $T_0+T_1$ periods of time. Units are indexed by $i = 1,2,\ldots J, J+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. The single treated unit in our context could be understood as an “aggregate” of multiple treated units; see Section (ref) below for more discussion.
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. Two implicit assumptions are imposed: 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). Then, the observed outcome $Y_{it}$ is \[Y_{it}=
. \]
The causal quantity of interest is the difference between the outcome path taken by the treated unit, and the path it would have taken in the absence of the treatment: \[ \tau_t:=Y_{1t}(1)-Y_{1t}(0), \qquad t>T_0. \] We view the two potential outcomes $Y_{1t}(1)$ and $Y_{1t}(0)$ as random variables, which implies that $\tau_t$ is a random quantity as well, corresponding to the treatment effect on a single treated unit. This contrasts with other analysis that regards the treatment effect as a fixed parameter Abadie_2021_JEL.
The potential outcome $Y_{1t}(1)$ of the treated unit is observed after the treatment. To recover the treatment effect $\tau_t$, it is necessary to have a “good” prediction of the counterfactual outcome $Y_{1t}(0)$ of the treated after the intervention. The idea of the synthetic control method is to find a vector of weights $\mathbf{w} = (w_2,w_3,\dots,w_{J+1})'$ such that a given loss function is minimized under constraints, only using pre-intervention observations. Given the resulting set of constructed weights $\widehat{\mathbf{w}}$, the treated unit's counterfactual (potential) outcome is calculated as $\widehat{Y}_{1t}(0) = \sum_{i=2}^{J+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. In what follows, we briefly describe different approaches for point estimation/prediction leading to $\widehat{Y}_{1t}(0)$, and then summarize the uncertainty quantification methods to complement those predictions.
Building on the canonical synthetic control setup, we can consider other settings involving multiple treated units with possibly staggered treatment adoption. In particular, we briefly discuss three potential extensions of practical interest.
The package scpi allows for estimation/prediction of treatment effects and uncertainty quantification via prediction intervals for the more general synthetic control settings discussed above. However, in order to streamline the exposition, the rest of this article focuses on the case of a single treated unit. See Cattaneo-Feng-Palomba-Titiunik_2022_wp for a formal treatment of more general staggered adoption problems, and its supplemental appendix for further details on how the package scpi can be used in settings with multiple treatment units and staggered treatment adoption. Our companion replication files illustrate both the canonical single treated unit framework and the generalizations discussed above.
We 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,\ldots, 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, during the pre-intervention period, for all $1\leq l\leq M$ and some convex feasibility sets $\mathcal{W}$ and $\mathcal{R}$ that capture the restrictions imposed. Specifically, we consider the following optimization problem:
where \[\underbrace{\mathbf{A}}_{T_0\cdot M\times 1} =
, \quad \underbrace{\mathbf{B}}_{T_0\cdot M \times J} =
,\quad \underbrace{\mathbf{C}}_{T_0\cdot M \times K\cdot M}=
\] and $\mathbf{V}$ is a $T_0\cdot M\times T_0\cdot M$ weighting matrix reflecting the relative importance of different equations and time periods.
From (ref), we can define the pseudo-true residual $\mathbf{u}$ as
where $\mathbf{w}_0$ and $\mathbf{r}_0$ denote the mean squared error population analog of $\widehat{\mathbf{w}}$ and $\widehat{\mathbf{r}}$. As discussed in the next section, the proposed prediction intervals are valid conditional on some information set $\mathscr{H}$. Thus, $\mathbf{w}_0$ and $\mathbf{r}_0$ above are viewed as the (possibly constrained) best linear prediction coefficients conditional on $\mathscr{H}$. We do not attach any structural meaning to $\mathbf{w}_0$ and $\mathbf{r}_0$: they are only (conditional) pseudo-true values whose meaning should be understood in context, and are determined by the assumptions imposed on the data generating process. In other words, we allow for misspecification when constructing the synthetic control weights $\widehat{\mathbf{w}}$, as this is the most likely scenario in practice.
Given the constructed weights $\widehat{\mathbf{w}}$ and coefficients $\widehat{\mathbf{r}}$, the counterfactual outcome at the post-treatment period $T$ for the treated unit, $Y_{1T}(0)$, is predicted by
where $\mathbf{x}_T\in\mathbb{R}^{J}$ 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 in practice it is often the case that $\mathbf{x}_T=(Y_{2T}(0), \cdots, Y_{(J+1)T}(0))'$ and $\mathbf{g}_T$ is excluded when $\mathbf{C}$ is not specified.
The next section discusses implementation details leading to $\widehat{Y}_{1T}(0)$, including the choice of feasibility sets $\mathcal{W}$ and $\mathcal{R}$, weighting matrix $\mathbf{V}$, and additional covariates $\mathbf{C}$.
The function scdata() in scpi prepares the data for point estimation/prediction purposes. This function takes as input an object of class DataFrame and outputs an object of class scpi_data containing the matrices $\mathbf{A},\mathbf{B},\mathbf{C}$ described above, and a matrix of post-treatment predictors $\mathbf{P}=(\mathbf{p}_{T_0+1}, \cdots, \mathbf{p}_{T_0+T_1})'$. The user must provide a variable containing a unit identifier (id.var), a time variable (time.var), an outcome variable (\texttt{outcome.var}), the features to be matched (\texttt{features}), the treated unit (\texttt{unit.tr}), the control units (\texttt{unit.co}), the pre-treatment periods (\texttt{period.pre}), and the post-treatment periods (\texttt{period.post}). These options completely specify $\mathbf{A},\mathbf{B}$, and $\mathbf{P}$. The user can also control the form of $\mathcal{R}$ in (ref) or, equivalently, the form of $\mathbf{C}$, through the options \texttt{cov.adj} and \texttt{constant}. The former option allows the user to flexibly specify covariate adjustment feature by feature, while the latter option introduces a column vector of ones of size $T_0\cdot M$ in $\mathbf{C}$. If $M=1$, this is a simple constant term, but if $M\geq 2$ it corresponds to an intercept which is common across features.
The use of the options cov.adj and constant is best explained through some examples. If the user specifies only one feature ($M=1$), then cov.adj can be an unnamed list:
This particular choice includes a constant term and a linear time trend in $\mathbf{C}$. If instead multiple features ($M \geq 2$) are used to find the synthetic control weights $\widehat\mathbf{w}$, then cov.adj allows for feature-specific covariate adjustment. For example, in a two-feature setting ($M=2$), the code
specifies $\mathbf{C}$ as a block diagonal matrix where the first block $\mathbf{C}_1$ contains a constant term and a trend, while the second block $\mathbf{C}_2$ contains only a trend. If the user wants all features to share the same covariate adjustment, then it is sufficient to input a list with a unique element:
This specification creates a block diagonal matrix $\mathbf{C}$ with identical blocks. In the same example with $M=2$, if constant <- TRUE and cov.adj <- NULL, then $\mathbf{C}$ would not be block diagonal, but rather a column vector of ones of size $2 T_0$.
Finally, if $\mathbf{A}$ and $\mathbf{B}$ form a cointegrated system, by setting the option cointegrated.data to TRUE in scdata(), the matrix $\mathbf{P}$ is prepared in such a way that the function scpi() will properly handle in-sample and out-of-sample uncertainty quantification (see sections (ref) and (ref)).
Once all the design matrices $\mathbf{A},\mathbf{B},\mathbf{C},$ and $\mathbf{P}$ have been created, we can proceed with point estimation/prediction of the counterfactual outcome of interest via the function scest().
The form of the feasibility set $\mathcal{W}$ in (ref) or, equivalently, the constraints imposed on the weights $\mathbf{w}$, can be set using the option w.constr. The package allows for the following family of constraints:
where the inequality constraint on the norm can be made an equality constraint as well. The user can specify the desired form for $\mathcal{W}$ through a list to be passed to the option w.constr:
The four lines above create $\mathcal{W}_1 = \mathbb{R}^J$, $\mathcal{W}_2 = \{\mathbf{w} \in \mathbb{R}_+^{J}: ||\mathbf{w}||_1= 1\}$, $\mathcal{W}_3 = \{\mathbf{w} \in \mathbb{R}^{J}: ||\mathbf{w}||_2\leq 1\}$, and $\mathcal{W}_4 = \{\mathbf{w} \in \mathbb{R}^{J}: ||\mathbf{w}||_1= 1, ||\mathbf{w}||_2\leq 1\}$, respectively. In greater detail,
Popular constraints can be called explicitly using the option name in the list passed to w.constr. Table (ref) gives prototypical examples of such constraints.
In particular, specifying list(name = `simplex', Q = 1) gives the standard constraint used in the canonical synthetic control method, that is, computing weights in (ref) such that they are non-negative and sum up to one, and without including an intercept. This is the default in the function scest() (and scpi()). The following snippet showcases how each of these five constraints can be called automatically through the option name and manually through the options p, Q, \texttt{Q2}, \texttt{lb}, and \texttt{dir}. In the snippet, \texttt{Q} and \texttt{Q2} are set to 1 for ridge and L1-L2 constraints, respectively, for simplicity, but to replicate the results obtained with the option \texttt{name} one should input the proper \texttt{Q} according to the rules of thumb described further below.
{\singlespacing\lstinputlisting{txt/tablecomparison.txt}}
Using the option w.constr in scest() (or scpi()) and the options cov.adj and constant in scdata() appropriately, i.e., setting $\mathcal{W}$ and $\mathcal{R}$ in (ref), many synthetic control estimators proposed in the literature can be implemented. Table (ref) provides a non-exhaustive list of such examples.
We provide rule-of-thumb choices of the tuning parameter $Q$ for Lasso- and Ridge-type constraints.
In case of missing values, we adopt different strategies depending on which units have missing entries and when these occur.
Following \citet*{Cattaneo-Feng-Titiunik_2021_JASA} and \citet*{Cattaneo-Feng-Palomba-Titiunik_2022_wp}, we view the quantity of interest $\tau_{T}$ within the synthetic control framework as a random variable, and hence we refrain from calling it a parameter. Building an analogy with the concept of estimand (or parameter of interest), we refer to $\tau_T$ as a predictand. Consequently, we prefer to call $\widehat{\tau}_{T}=Y_{1T}(1)-\widehat{Y}_{1T}(0)$ based on (ref) a prediction of $\tau_{T}$ rather than an estimator of it, and our goal is to characterize the uncertainty of $\widehat{\tau}_T$ by building prediction intervals rather than confidence intervals. In practice, it is appealing to construct prediction intervals that are valid conditional on a set of observables. We let $\mathscr{H}$ be an information set generated by all features of control units and covariates used in the synthetic control construction, i.e., $\mathbf{B}$, $\mathbf{C}$, $\mathbf{x}_T$, and $\mathbf{g}_T$.
We first decompose the potential outcome of the treated unit based on $\mathbf{w}_0$ and $\mathbf{r}_0$ introduced in (ref):
where $e_T$ is defined by construction. In our analysis, $\mathbf{w}_0$ and $\mathbf{r}_0$ are assumed to be (possibly) random quantities around which $\widehat{\mathbf{w}}$ and $\widehat{\mathbf{r}}$ are concentrating in probability, respectively. Then, the distance between the predicted treatment effect on the treated and the target population one is
where $e_T$ is the out-of-sample error coming from misspecification along with any additional noise occurring at the post-treatment period $T > T_0$, and the term $\mathbf{p}_T'(\widehat\boldsymbol{\beta}-\boldsymbol{\beta}_0)$ is the in-sample error coming from the construction of the synthetic control weights. Our goal is to find probability bounds on the two terms separately to give uncertainty quantification: for some pre-specified levels $\alpha_1, \alpha_2\in(0,1)$, with high probability over $\mathscr{H}$,
It follows that these probability bounds can be combined to construct a prediction interval for $\tau_T$ with conditional coverage at least $1-\alpha_1-\alpha_2$: with high probability over $\mathscr{H}$, \[ \mathbb{P}\big[\widehat{\tau}_T+M_{1,\mathtt{L}}-M_{2,\mathtt{U}} \leq \tau_T\leq \widehat{\tau}_T+M_{1,\mathtt{U}}-M_{2,\mathtt{L}} \big| \mathscr{H} \big] \geq 1-\alpha_1-\alpha_2. \]
\citet*{Cattaneo-Feng-Titiunik_2021_JASA} provide a principled simulation-based method for quantifying the in-sample uncertainty coming from $\mathbf{p}_T'(\widehat\boldsymbol{\beta}-\boldsymbol{\beta}_0)$. Let $\mathbf{Z}=(\mathbf{B}, \mathbf{C})$ and $\mathbf{D}$ be a non-negative diagonal (scaling) matrix of size $d$, possibly depending on the pre-treatment sample size $T_0$. Since $\widehat{\boldsymbol{\beta}}$ solves (ref), $\widehat{\boldsymbol{\delta}}:=\mathbf{D}(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_0)$ is 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{V}\mathbf{Z}\mathbf{D}^{-1}$, $\widehat{\boldsymbol{\gamma}}'=\mathbf{u}'\mathbf{V}\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}\}$. Recall that the information set conditional on which our prediction intervals are constructed contains $\mathbf{B}$ and $\mathbf{C}$. Thus, $\widehat\mathbf{Q}$ can be taken as fixed, and we need to characterize the uncertainty of $\widehat\boldsymbol{\gamma}$.
We construct a simulation-based criterion function accordingly:
where $\widehat{\boldsymbol{\Sigma}}$ is some estimate of $\boldsymbol{\Sigma}=\mathbb{V}[\widehat\boldsymbol{\gamma}|\mathscr{H}]$ and $\mathsf{N}(0, \widehat\boldsymbol{\Sigma})$ represents the normal distribution with mean $\mathbf{0}$ and variance-covariance matrix $\widehat\boldsymbol{\Sigma}$. In practice, the criterion function $\ell^{\star}(\cdot)$ can be simulated by simply drawing normal random vectors $\mathbf{G}^{\star}$.
Since the original constraint set $\Delta$ is infeasible, we need to construct a constraint set $\Delta^{\star}$ used in simulation that is close to $\Delta$. Specifically, define the distance between a point $\mathbf{a}\in\mathbb{R}^d$ and a set $\Lambda\subseteq\mathbb{R}^d$ by \[ \mathrm{dist}(\mathbf{a}, \Lambda)=\inf_{\lambda\in\Lambda}\|\mathbf{a}-\boldsymbol{\lambda}\|, \] where $\|\cdot\|$ is a generic $\ell_p$ vector norm on $\mathbb{R}^d$ with $p\geq 1$ (e.g., Euclidean norm or $\ell_1$ norm). We require
where $\mathcal{B}(\mathbf{0}, \varepsilon)$ is an $\varepsilon$-neighborhood around zero for some $\varepsilon>0$. In words, every point in the infeasible constraint set $\Delta$ has to be sufficiently close to the feasible constraint set used in simulation. We discuss below a principled strategy for constructing $\Delta^\star$, which allows for both linear and non-linear constraints in the feasibility set. Section (ref) provides details on how $\Delta^{\star}$ is constructed and implemented in the scpi package.
Given the feasible criterion function $\ell^\star(\cdot)$ and constraint set $\Delta^\star$, we let
conditional on the data. Under mild regularity conditions, for a large class of synthetic control predictands (ref), with high probability over $\mathscr{H}$, \[ \mathbb{P}\big[M_{1,\mathtt{L}}\leq\mathbf{p}_T'(\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}_0)\leq M_{1,\mathtt{U}} \; \big| \; \mathscr{H} \big]\geq 1-\alpha_1, \] up to some small loss of the (conditional) coverage probability. Importantly, this conclusion holds whether the data are stationary or non-stationary and whether the model is correctly specified (i.e., $\mathbb{E}[\mathbf{u}|\mathscr{H}]=0$) or not. If constraints imposed are non-linear, an additional adjustment to this bound may be needed to ensure the desired coverage.
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}$. Naturally, in order to have a proper bound on $e_T$, it is necessary to determine certain features of its conditional distribution $F_{e_T}(\mathsf{e}) = \mathbb{P}[e_T \leq \mathsf{e} | \mathscr{H}]$. In this section, we outline principled but agnostic approaches to quantify the uncertainty introduced by the post-treatment unobserved shock $e_T$. Since formalizing the validity of our methods usually requires strong assumptions, we also recommend a generic sensitivity analysis to incorporate out-of-sample uncertainty into the prediction intervals. See Section (ref) and Section (ref), in particular Figure (ref) with the corresponding snippet of R code, for further clarifications on how to carry out sensitivity analysis on $e_T$.
Using any of the above methods, we have the following probability bound on $e_T$: \[ \mathbb{P}\big[M_{2,\mathtt{L}}\leq e_T\leq M_{2,\mathtt{U}} \; \big| \; \mathscr{H} \big]\geq 1-\alpha_2. \]
We now discuss the implementation details. The function scpi(), through various options, allows the user to specify different approaches to quantify in-sample and out-of-sample uncertainty based on the methods described above. Most importantly, scpi() permits modelling separately the in-sample error $\mathbf{p}_T'(\widehat\boldsymbol{\beta}-\boldsymbol{\beta}_0)$ and the out-of-sample error $e_T$. In addition, the user can provide bounds on them manually with the options w.bounds and e.bounds, respectively, which can be useful for sensitivity analysis in empirical applications.
In-sample uncertainty stems from the prediction of $\mathbf{p}_T'(\widehat\boldsymbol{\beta}-\boldsymbol{\beta}_0)$, and its quantification reduces to determining $M_{1,\mathtt{L}}$ and $M_{1,\mathtt{U}}$. We first review the methodological proposals for constructing the constraint set $\Delta^\star$ used in simulation discussed in \citet*{Cattaneo-Feng-Palomba-Titiunik_2022_wp}, and then present the main procedure for constructing bounds on the in-sample error.
\uline{Constructing $\Delta^\star$}. Our in-sample uncertainty quantification requires the centered and scaled constraint feasibility set $\Delta$ to be locally identical to (or, at least, well approximated by) the constraint set $\Delta^\star$ used in simulation described in (ref), in the sense of (ref). Suppose that $$ \mathcal{W}\times\mathcal{R}=\Big\{\boldsymbol{\beta}\in\mathbb{R}^d:\bm{m}_{\mathtt{eq}}(\boldsymbol{\beta})=\bm{0},\bm{m}_{\mathtt{in}}(\boldsymbol{\beta})\leq \bm{0}\Big\}, $$ where $\bm{m}_{\mathtt{eq}}(\cdot)\in\mathbb{R}^{d_{\mathtt{eq}}}$ and $\bm{m}_{\mathtt{in}}(\cdot)\in\mathbb{R}^{d_{\mathtt{in}}}$ and denote the $j$th constraint in $\bm{m}_{\mathtt{in}}(\cdot)$ as $m_{\mathtt{in},j}(\cdot)$. Given tuning parameters $\varrho_j>0$, $j=1, \cdots, d_{\mathtt{in}}$, let $\mathcal{B}$ be the set of indices for the inequality constraints such that $m_{\mathtt{in}, j}(\hat{\boldsymbol{\beta}})>-\varrho_j$. Then, we construct $\Delta^\star$ as \[ \Delta^\star=\Big\{ \mathbf{D}(\boldsymbol{\beta}-\widehat{\boldsymbol{\beta}}):\bm{m}_{\mathtt{eq}}(\boldsymbol{\beta})=\bm{0}, m_{\mathtt{in},j}(\boldsymbol{\beta})\leq m_{\mathtt{in},j}(\widehat{\boldsymbol{\beta}}) \text{ for } j\in\mathcal{B}, \text{ and } m_{\mathtt{in},l}(\boldsymbol{\beta})\leq \bm{0} \text{ for } l\notin\mathcal{B} \Big\}. \] In practice, we need to choose possibly heterogeneous parameters $\varrho_j$, $j=1,\ldots, d_\mathtt{in}$, for different inequality constraints. Our proposed choice of $\varrho_j$ is \[ \varrho_j:=\Big\|\frac{\partial}{\partial\boldsymbol{\beta}}m_{\mathtt{in},j}(\widehat{\boldsymbol{\beta}})\Big\|_1\times \varrho, \quad j=1,\ldots,d_\mathtt{in}, \] for some parameter $\varrho$ where $\|\cdot\|_1$ denotes the $\ell_1$-norm. We estimate $\varrho$ according to the following formula if $M=1$: \[ \varrho = \mathcal{C} \frac{\log(T_0)^c}{T_0^{1/2}},\] where $c=1/2$ if the data are i.i.d. or weakly dependent, and $c=1$ if $\mathbf{A}$ and $\mathbf{B}$ form a cointegrated system, while $\mathcal{C}$ is one of the following: \[\mathcal{C}_1 = \frac{\widehat{\sigma}_{u}}{\min_{1\leq j \leq J}\widehat{\sigma}_{b_j}}, \qquad \mathcal{C}_2 = \frac{\max_{1\leq j \leq J}\widehat{\sigma}_{b_j}\widehat{\sigma}_{u}}{\min_{1\leq j \leq J}\widehat{\sigma}^2_{b_j}}, \qquad \mathcal{C}_3 = \frac{\max_{1\leq j \leq J}\widehat{\sigma}_{b_ju}}{\min_{1\leq j \leq J}\widehat{\sigma}^2_{b_j}}, \] with $\mathcal{C}_1$ as the default. $\widehat\sigma_{b_j,u}$ is the estimated (unconditional) covariance between the pseudo-true residual $\mathbf{u}$ and the feature of the $j$th control unit $\mathbf{B}_{j,1}$, and $\widehat\sigma_u$ and $\widehat\sigma_{b_j}$ are the estimated (unconditional) standard deviation of, respectively, $\mathbf{u}$ and $\mathbf{B}_{j,1}$. In the case of multiple features ($M>1$), the package employs the same construction above after stacking the data.
\uline{Degrees-of-Freedom Correction}. Our uncertainty quantification strategy requires an estimator of the conditional variance $\mathbb{V}[\mathbf{u} | \mathscr{H}]$, which may rely on the effective degrees of freedom $\mathsf{df}$ of the synthetic control method. In general, there exists no exact correspondence between the degrees of freedom and the number of parameters in a fitting model Ye_1998_JASA. Therefore, the estimated degrees of freedom $\widehat{\mathsf{df}}$ are defined according to the chosen constraint sets for $\boldsymbol{\beta}$ underlying the estimation procedure in (ref):
\uline{Main procedure}. Given the constraint set $\Delta^\star$, the main procedure for computing the upper and lower bounds on the in-sample error is as follows:
\uline{Execution Speed and Parallelization}. Steps 3 and 4 of the procedure above are the most computationally intensive and we optimize them in two ways. First, to solve the optimization problem in (ref), scpi relies on ECOS, an efficient solver for conic problems \citep*{Domahidi-Chu-Boyd_2013_ECOS, Fu-Narasimhan-Boyd_2020_JSS}. See \citet*{Cattaneo-Feng-Palomba-Titiunik_2022_wp} for more details on how to cast the different constrained SC methods into conic optimization problems. To give the reader a sense of the speed improvement, Table (ref) compares the execution speed of the conic solver we rely on (first column) with other two popular optimizers in R. The first row of the table reports the median computation time of each optimizer, whereas the second row shows the inter-quartile range. On the one hand, using a conic solver in place of a solver for more generic optimization programs (like nloptr) makes our software 4 times faster. On the other hand, our software is tailored to rewrite the SC problem as a conic problem. This gives a 300-fold gain in speed when compared to CVXR, which relies on ECOS but is meant for prototyping generic optimization problems in conic form.
Second, scpi can be sped up further by efficient parallelization of the tasks performed through the package parallel which assigns different simulations to different cores. If $\mathtt{N}_{\text{cores}}$ cores are used, the final execution time would be approximately $\mathtt{T}_{\text{exec}}/\mathtt{N}_{\text{cores}}$, where $\mathtt{T}_{\text{exec}}$ is the execution time when a single core is used.
To quantify the uncertainty coming from $e_T$, we need to impose some probabilistic structure that allows us to model the distribution $\mathbb{P}[e_T\leq \mathtt{e} |\mathscr{H}]$ and, ultimately, estimate $M_{2,\mathtt{L}}$ and $M_{2,\mathtt{U}}$. We discussed three different alternative approaches: (i) non-asymptotic bounds; (ii) location-scale model; and (iii) quantile regression. The user can choose the preferred way of modeling $e_T|\mathscr{H}$ by setting the option e.method to either `gaussian', `ls', or `qreg'.
The user can also choose the information used to estimate (conditional) moments or quantiles of $e_T|\mathscr{H}$. Practically, we allow the user to specify a design matrix $\mathbf{D}_{\mathbf{e}}$ that is then used to run the appropriate regressions depending on the approach requested. By default, we set $\mathbf{D}_{\mathbf{e}} = [\mathbf{B}^\star_1\;\; \mathbf{C}_1]$. Alternatively, the matrix $\mathbf{D}_{\mathbf{e}}$ can be provided directly through the option e.design or by specifying the lags of $\mathbf{B}^\star_1$ (e.lags) and/or the order of the fully interacted polynomial in $\mathbf{B}^\star_1$ (e.order). If the user specifies e.lags = 0 and e.order = 2, then $\mathbf{D}_{\mathbf{e}}$ contains $\mathbf{B}^\star_1$, $\mathbf{C}_1$, and all the unique second-order terms generated by the interaction of the columns of $\mathbf{B}^\star_1$. If instead e.order = 0 and \texttt{e.lags = 0} are set, then $\widehat{\mathbb{E}}[e_T|\mathscr{H}]$ and $\widehat{\mathbb{V}}[e_T|\mathscr{H}]$ are estimated using the sample average and the sample variance of $e_T$ using the pre-intervention data. Recall that if the option \texttt{cointegrated.data} is set to \texttt{TRUE}, $\mathbf{B}^\star_1$ is formed using the first differences of the columns in $\mathbf{B}_1$. Finally, the user can specify $\alpha_2$ with the option \texttt{e.alpha}.
Up to this point, we focused on prediction intervals that possess high coverage for the individual treatment effect in each period. However, it may be desirable to have prediction intervals that have high simultaneous coverage for several periods, usually known as simultaneous prediction intervals in the literature. In other words, our final goal is to construct a sequence of intervals $\{\mathcal{I}_t: T_0+1\leq t\leq T_0+L\}$ for some $1\leq L \leq T_1$ such that with high probability over $\mathscr{H}$, \[ \mathbb{P}\big[\tau_{t}\in\mathcal{I}_t, \text{ for all } T_0+1\leq t\leq T_0+L\,\big|\, \mathscr{H}\big]\geq 1-\alpha_1-\alpha_2. \]
To construct such intervals, we need to generalize the procedures described above to quantify the in-sample error (Section (ref)) and the out-of-sample error (Section (ref)).
With regard to the in-sample uncertainty, we handle two separate cases. On the one hand, if the constraints in $\Delta$ are linear (e.g., simplex or lasso), then
which guarantees that with high probability over $\mathscr{H}$ \[ \mathbb{P}\big[M_{1,\mathtt{L}}\leq \mathbf{p}_{t}'(\boldsymbol{\beta}_0-\widehat{\boldsymbol{\beta}})\leq M_{1,\mathtt{U}}, \text{ for all } T_0+1\leq t\leq T_0+L\,\big|\, \mathscr{H}\big]\geq 1-\alpha_1. \] On the other hand, if $\Delta$ includes non-linear constraints (e.g., constraints involving the $\ell_2$ norm), it is necessary to decrease the lower bound $M_{1,\mathtt{L}}$ and increase the upper bound $M_{1,\mathtt{U}}$ by some quantity $\varepsilon_{\Delta,t}>0$ for each $T_0+1\leq t\leq T_0+L$. To give an example of what $\varepsilon_{\Delta,t}$ looks like, in the case of ridge-type constraints we have \[ \varepsilon_{\Delta,t}=\|\mathbf{p}_t\|_1\times (2\|\widehat{\boldsymbol{\beta}}\|_2)^{-1}\times \varrho^2, \] and see \citet*{Cattaneo-Feng-Palomba-Titiunik_2022_wp} for more general cases. With regard to the out-of-sample uncertainty, our proposed strategy is a generalization of “Approach 1” in Section (ref): find $M_{2,\mathtt{L},t}$ and $M_{2,\mathtt{U},t}$ such that with high probability over $\mathscr{H}$, \[ \mathbb{P}\big[M_{2,\mathtt{L},t}\leq e_t\leq M_{2,\mathtt{U},t}, \text{ for all } T_0+1\leq t\leq T_0+L\,\big|\, \mathscr{H}\big]\geq 1-\alpha_2. \] Suppose that each $e_{t}$, $T_0+1\leq t\leq T_0+L$, is sub-Gaussian conditional on $\mathscr{H}$ (not necessarily independent over $t$) with sub-Gaussian parameters $\sigma_{\mathscr{H},t}\leq \sigma_{\mathscr{H}}$ for some $\sigma_{\mathscr{H}}$. Then, we can take \[ M_{2,\mathtt{L},t}:=\mathbb{E}[e_t|\mathscr{H}]-\sqrt{2\sigma_{\mathscr{H}}^2\log (2L/\alpha_2)} \qquad\text{and}\qquad M_{2,\mathtt{U},t}:=\mathbb{E}[e_t|\mathscr{H}]+\sqrt{2\sigma_{\mathscr{H}}^2\log (2L/\alpha_2)}. \] We can see that, compared to what we had for “Approach 1”, there is an extra term, $\sqrt{\log L}$, which makes the simultaneous prediction intervals longer.
While the three approaches for out-of-sample uncertainty quantification described in Section (ref) are simple and intuitive, their validity requires potentially strong assumptions on the underlying data generating process that links the pre-treatment and post-treatment data. Such assumptions are difficult to avoid because the ultimate goal is to learn about the statistical uncertainty introduced by a single unobserved random variable after the treatment/intervention is deployed, that is, $e_T|\mathscr{H}$ for some $T>T_0$. Without additional data availability, or specific modelling assumptions allowing for transferring information from the pre-treatment period to the post-treatment period, it is difficult to formally construct $M_{2,\mathtt{L}}$ and $M_{2,\mathtt{U}}$ using data-driven methods.
We suggest approaching the out-of-sample uncertainty quantification as a principled sensitivity analysis, using the approaches above as a starting point. Given the formal and detailed in-sample uncertainty quantification described previously, it is natural to progressively enlarge the final prediction intervals by adding additional out-of-sample uncertainty to 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$ statistically insignificant? Using the approaches above, or similar ones, it is possible to construct natural initial benchmarks. For instance, to implement Approach 1, one can use the pre-treatment outcomes or synthetic control residuals to obtain a “reasonable” benchmark estimate of the sub-Gaussian parameter $\sigma_{\mathscr{H}}$ and then progressively enlarge or shrink this parameter to check the robustness of the conclusion. Alternatively, in specific applications, natural levels of uncertainty for the outcomes of interest could be available, and hence used to tabulate the additional out-of-sample uncertainty. We illustrate this approach in Section (ref).
We showcase the features of the package scpi using real data. For comparability purposes, we employ the canonical dataset in the synthetic control literature on the economic consequences of the 1990 German reunification Abadie_2021_JEL, and focus on quantifying the causal impact of the German reunification on GDP per capita in West Germany. Thus, we compare the post-reunification outcome for West Germany with the outcome of a synthetic control unit constructed using 16 OECD countries from 1960 to 1990. Using the notation introduced above, we have $T_0 = 31$ and $J=16$. The only feature we exploit to construct the synthetic control is yearly GDP per capita, and we add a constant term for covariate adjustment. Thus $M=1$ and $K=1$, and $\mathcal{R}=\mathbb{R}$. We explore the effect of the reunification from 1991 to 2003, hence $T_1 = 13$. Finally, we treat the time series for West Germany and those countries in the donor pool as a cointegrating system. Given this information, the command scdata() prepares all the matrices needed in the synthetic control framework described above ($\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$ and $\mathbf{P}$), and returns an object that must be used as input in either scest() to predict $\widehat{Y}_{1T}(0),$ $T>T_0$, or scpi() to conduct inference on $\widehat{\tau}_{T},$ $T>T_0$.
We first call scdata() to transform any data frame into an object of class “scpi_data”.
{\singlespacing\lstinputlisting{txt/1scdata.txt}}
After having prepared the data, the next step involves choosing the desired constraint set $\mathcal{W}$ to construct the vector of weights $\mathbf{w}$. We consider the canonical synthetic control method and thus search for optimal weights in $\mathcal{W} = \{\mathbf{w} \in \mathbb{R}^J_+: ||\mathbf{w}||_1 = 1\}$. Such constraint set is the default in scest() and, consequently, in scpi(), as the latter internally calls the former to construct $\mathbf{w}$. The snippet below illustrates how to call scest() and reports the results displayed in the console with the summary() method.
{\singlespacing\lstinputlisting{txt/2scest.txt}}
The next step is uncertainty quantification using scpi(). In this case, we quantify the in-sample and out-of-sample uncertainty the same way, using $\mathbf{B}$ and $\mathbf{C}$ as the conditioning set in both cases. To do so, it suffices to set the order of the polynomial in $\mathbf{B}$ to 1 (u.order <- 1 and e.order <- 1) and not include lags (u.lags <- 0 and e.lags <- 0). Furthermore, by specifying the option u.miss <- TRUE, we take into account that the conditional mean of $\mathbf{u}$ might differ from 0. This option, together with \texttt{u.sigma <- "HC1"}, specifies the following estimator of $\mathbb{V}[\mathbf{u}|\mathscr{H}]$:
\[\widehat\mathbb{V}[\mathbf{u}|\mathscr{H}] = \operatorname*{diag}\Big( \mathtt{vc}_1^{(1)}(\widehat{\mathbf{u}}_1 -\widehat\mathbb{E}[\mathbf{u}_1 |\mathscr{H}])^2,\cdots, \mathtt{vc}_{T_0}^{(1)}(\widehat{\mathbf{u}}_{T_0} -\widehat\mathbb{E}[\mathbf{u}_{T_0} |\mathscr{H}])^2 \Big). \]
Finally, by selecting e.method <- "gaussian", we perform out-of-sample uncertainty quantification treating $e_T$ as sub-gaussian conditional on $\mathbf{B}$ and $\mathbf{C}$. As a last step, we visualize the constructed synthetic control and compare it with the observed time series for the treated unit, taking advantage of the function scplot().
{\singlespacing\lstinputlisting{txt/3scpi.txt}}
Figure (ref) displays the plot resulting from the scplot call. The vertical bars are 90% prediction intervals, where the non-coverage error rate is halved between the out-of-sample and the in-sample uncertainty quantification, i.e. $\alpha_1=\alpha_2=0.05$.
We also conduct the same exercise using different choices of $\mathcal{W}$ (see Table (ref)). In particular, we construct weights and compute prediction intervals using four other specifications: (i) a lasso-type constraint (Figure (ref)), (ii) a ridge-type constraint (Figure (ref)), (iii) no constraint (Figure (ref)), and (iv) an L1-L2 constraint.
{\singlespacing\lstinputlisting{txt/4scpi.txt}}
Section (ref) clarified the need for some additional sensitivity analysis when it comes to out-of-sample uncertainty quantification. Figure (ref) shows the impact of shrinking and enlarging $\widehat\sigma_{\mathscr{H}}$ on the prediction intervals for $Y_{1t}(0)$, $t=1997$, when we assume that $e_t$ is sub-Gaussian conditional on $\mathscr{H}$. As shown in the figure, the predicted treatment effect $\widehat\tau_{\mathtt{1997}}$ remains different from zero with high probability over $\mathscr{H}$ even doubling $\widehat\sigma_\mathscr{H}$ .
Finally, the package offers the possibility to match the treated unit and the synthetic unit using multiple features and the possibility to compute simultaneous prediction intervals. If we want to match West Germany and the synthetic unit not only on GDP per capita but also on trade openness ($M=2$) and include joint prediction intervals, we can simply modify the object scpi_data as follows.
{\singlespacing\lstinputlisting{txt/5scpimulti.txt}}
Results are reported in Figure (ref), where blue shaded areas depict 90% simultaneous prediction intervals for periods from 1991 to 2004.
This article introduced the R software package scpi, which implements point estimation/prediction and inference/uncertainty quantification procedures for synthetic control methods. The package is also available in the Stata and Python statistical platforms, as described in the appendices. Further information can be found at \url{https://nppackages.github.io/scpi/}.
We thank Alberto Abadie and Bartolomeo Stellato for many insightful discussions. Cattaneo and Titiunik gratefully acknowledge financial support from the National Science Foundation (SES-2019432), and Cattaneo gratefully acknowledges financial support from the National Institute of Health (R01 GM072611-16).