scpi: Uncertainty Quantification for Synthetic Control Methods

The synthetic control method was introduced by Abadie and Gardeazabal (2003), 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 Abadie (2021) for a review on synthetic control methods, and Abadie and Cattaneo (2018) 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 (Abadie and L’Hour, 2021) and staggered treatment adoption (Ben-Michael, Feller and Rothstein, 2022), while examples of the latter include employing different constrained estimation methods (see Table 3 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 (Abadie, Diamond and Hainmueller, 2010). Other approaches include methods based on large-sample approximations for disaggregated data under correctly specified factor-type models (Li, 2020), time-series permutation-based inference (Chernozhukov, Wüthrich and Zhu, 2021), large-sample approximations for high-dimensional penalization methods (Masini and Medeiros, 2021), and cross-sectional permutation-based inference in semiparametric duration-type settings (Shaikh and Toulis, 2021). A conceptually distinct approach to uncertainty quantification is proposed by Cattaneo, Feng and Titiunik (2021) and Cattaneo, Feng, Palomba and Titiunik (2022), 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 (Domahidi, Chu and Boyd, 2013; Fu, Narasimhan and Boyd, 2020); see also Boyd and Vandenberghe (2004) 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 A, and for Stata in Appendix B.

The R package scpi includes the following six functions:

• •

  scdata() and scdataMulti(). These functions take as input a DataFrame object and process it to prepare the data matrices used for point estimation/prediction and inference/uncertainty quantification. The function scdata() is specific to the single treated unit case, whereas scdataMulti() can be used with multiple treated units and/or when treatment is adopted in a staggered fashion. Both functions allow the user to specify multiple features of the treated unit(s) to be matched by the synthetic unit(s), as well as feature-specific covariate adjustment, and can handle both independent and identically distributed (i.i.d.) and non-stationary (cointegrated) data.

• •

  scest(). This function handles “scpi_data” objects produced with scdata() or “scpi_data_multi” objects produced with scdataMulti(), and then implements a class of synthetic control predictions/point estimators for quantification of treatment effects. The implementation allows for multiple features, with and without additional covariate adjustment, and for both stationary and non-stationary data. The allowed prediction procedures include unconstrained weighted least squares as well as constrained weighted least squares with simplex, lasso-type, ridge-type parameter space restrictions and combinations thereof (see Table 2 below).

• •

  scpi(). This function takes as input an “scpi_data” object produced with scdata() or an “scpi_data_multi” object produced with scdataMulti(), and then computes prediction intervals for a class of synthetic control predictions/point estimators for quantification of treatment effects. It relies on scest() for point estimation/prediction of treatment effects, and thus inherits the same functionalities of that function. In particular, scpi() is designed to be the main function in applications, offering both predictions/point estimators for treatment effects as well as inference/uncertainty quantification (i.e., prediction intervals) for synthetic control methods. The function also allows the user to model separately in-sample and out-of-sample uncertainty, offering a broad range of options for practice.

• •

  scplot() and scplotMulti(). These functions process objects whose class is either “scest” or “scpi”. These objects contain the results of the point estimation/prediction or uncertainty quantification methods, respectively. The commands build on the ggplot2 package in R (Wickham, 2016) to compare the time series for the outcome of the treated unit(s) with the outcome time series of the synthetic control unit, along with the associated uncertainty. The functions return a ggplot object that can be further modified by the user.

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 scpi(). The function scest() is useful in cases where only predictions/point estimators are of interest. Numerical illustrations are given in Section 5.

There are many R, Python, and Stata packages available for prediction/point estimation and inference using synthetic control methods; Table 1 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 Python, R, and Stata, which gives full portability across multiple statistical software and programming languages, and also the only package employing directly native conic optimization via the ECOS solver (see Table 4 for details).

The rest of the article is organized as follows. Section 2 introduces the canonical synthetic control setup, and also briefly discusses extensions to multiple treated units with possibly staggered treatment adoption. Section 3 gives a brief introduction to the theory and methodology underlying the point estimation/prediction for synthetic control methods, discussing implementation details. Section 4 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 5 showcases some of the functionalities of the package using a real-world dataset, and Section 6 concludes. The appendices illustrate the Python (Appendix A) and Stata (Appendix B) 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, 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 2.1 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

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:

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 (see Abadie, 2021, for references).

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})^{\prime}$ 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.

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})^{\prime}\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})^{\prime}$ 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})^{\prime}$ 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

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 (3.1), 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))^{\prime}$ 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}$.

Following Cattaneo, Feng and Titiunik (2021) and Cattaneo, Feng, Palomba and Titiunik (2022), 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 (3.3) 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 (3.2):

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}^{\prime}(\widehat{\bm{\beta}}-\bm{\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}$,