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.
96,339 characters · 34 sections · 75 citation commands
Robust Synthetic Control
On November 8, 2016 in the aftermath of several high profile mass-shootings, voters in California passed Proposition 63 in to law prop63. Prop. 63 “outlaw[ed] the possession of ammunition magazines that [held] more than 10 rounds, requir[ed] background checks for people buying bullets,” and was proclaimed as an initiative for “historic progress to reduce gun violence” prop63latimes. Imagine that we wanted to study the impact of Prop. 63 on the rates of violent crime in California. Randomized control trials, such as A/B testings, have been successful in establishing effects of interventions by randomly exposing segments of the population to various types of interventions. Unfortunately, a randomized control trial is not applicable in this scenario since only one California exists. Instead, a statistical comparative study could be conducted where the rates of violent crime in California are compared to a “control” state after November 2016, which we refer to as the post-intervention period. To reach a statistically valid conclusion, however, the control state must be demonstrably similar to California sans the passage of a Prop. 63 style legislation. In general, there may not exist a natural control state for California, and subject-matter experts tend to disagree on the most appropriate state for comparison.
As a suggested remedy to overcome the limitations of a classical comparative study outlined above, Abadie et al. proposed a powerful, data-driven approach to construct a “synthetic” control unit absent of intervention abadie1, abadie3, abadie2. In the example above, the synthetic control method would construct a “synthetic” state of California such that the rates of violent crime of that hypothetical state would best match the rates in California before the passage of Prop. 63. This synthetic California can then serve as a data-driven counterfactual for the period after the passage of Prop. 63. Abadie et al. propose to construct such a synthetic California by choosing a convex combination of other states (donors) in the United States. For instance, synthetic California might be 80% like New York and 20% like Massachusetts. This approach is nearly entirely data-driven and appeals to intuition. For optimal results, however, the method still relies on subjective covariate information, such as employment rates, and the presence of domain “experts” to help identify a useful subset of donors. The approach may also perform poorly in the presence of non-negligible levels of noise and missing data.
As the main result, we propose a simple, two-step robust synthetic control algorithm, wherein the first step de-noises the data and the second step learns a linear relationship between the treated unit and the donor pool under the de-noised setting. The algorithm is robust in two senses: first, it is fully data-driven in that it is able to find a good donor subset even in the absence of helpful domain knowledge or supplementary covariate information; and second, it provides the means to overcome the challenges presented by missing and/or noisy observations. As another important contribution, we establish analytic guarantees (finite sample analysis and asymptotic consistency) -- that are missing from the literature -- for a broader class of models.
{\em Robust algorithm.} A distinguishing feature of our work is that of de-noising the observation data via singular value thresholding. Although this spectral procedure is commonplace in the matrix completion arena, it is novel in the realm of synthetic control. Despite its simplicity, however, thresholding brings a myriad of benefits and resolves points of concern that have not been previously addressed. For instance, while classical methods have not even tackled the obstacle of missing data, our approach is well equipped to impute missing values as a consequence of the thresholding procedure. Additionally, thresholding can help prevent the model from overfitting to the idiosyncrasies of the data, providing a knob for practitioners to tune the “bias-variance” trade-off of their model and, thus, reduce their mean square error (MSE). From empirical studies, we hypothesize that thresholding may possibly render auxiliary covariate information (vital to several existing methods) a luxury as opposed to a necessity. However, as one would expect, the algorithm can only benefit from useful covariate and/or “expert” information and we do not advocate ignoring such helpful information, if available.
In the spirit of combatting overfitting, we extend our algorithm to include regularization techniques such as ridge regression and LASSO. We also move beyond point estimates in establishing a Bayesian framework, which allows one to quantitatively compute the uncertainty of the results through posterior probabilities.
{\em Theoretical performance.} To the best of our knowledge, ours is the first to provide finite sample analysis of the MSE for the synthetic control method, in addition to guarantees in the presence of missing data. Previously, the main theoretical result from the synthetic control literature (cf. abadie1, abadie3, abadie2) pertained to bounding the bias of the synthetic control estimator; however, the proof of the result assumed that the latent parameters, which live in the simplex, have a perfect pre-treatment match in the noisy predictor variables -- our analysis, on the other hand, removes this assumption. We begin by demonstrating that our de-noising procedure produces a consistent estimator of the latent signal matrix (Theorems (ref), (ref)), proving that our thresholding method accurately imputes and filters missing and noisy observations, respectively. We then provide finite sample analysis that not only highlights the value of thresholding in balancing the inherent “bias-variance” trade-off of forecasting, but also proves that the prediction efficacy of our algorithm degrades gracefully with an increasing number of randomly missing data (Theorems (ref), (ref), and Corollary (ref)). Further, we show that a computationally beneficial pre-processing data aggregation step allows us to establish the asymptotic consistency of our estimator in generality (Theorem (ref)).
Additionally, we prove a simple linear algebraic fact that justifies the basic premise of synthetic control, which has not been formally established in literature, i.e. the linear relationship between the treatment and donor units that exists in the pre-intervention continues to hold in post-intervention period (Theorem (ref)). We introduce a latent variable model, which subsumes many of the models previously used in literature (e.g. econometric factor models). Despite this generality, a unifying theme that connects these models is that they all induce (approximately) low rank matrices, which is well suited for our method.
{\em Experimental results.} We conduct two sets of experiments: (a) on existing case studies from real world datasets referenced in abadie1, abadie2, abadie3, and (b) on synthetically generated data. Remarkably, while abadie1, abadie2, abadie3 use numerous covariates and employ expert knowledge in selecting their donor pool, our algorithm achieves similar results without any such assistance; additionally, our algorithm detects subtle effects of the intervention that were overlooked by the original synthetic control approach. { Since it is impossible to simultaneously observe the evolution of a treated unit and its counterfactual}, we employ synthetic data to validate the efficacy of our method. Using the MSE as our evaluation metric, we demonstrate that our algorithm is robust to varying levels of noise and missing data, { reinforcing the importance of de-noising}.
Synthetic control has received widespread attention since its conception by Abadie and Gardeazabal in their pioneering work abadie3, abadie1. It has been employed in numerous case studies, ranging from criminology crime to health policy health to online advertisement to retail; other notable studies include abadie4, economic_liberalization, tax, turkey. In their paper on the state of applied econometrics for causality and policy evaluation, Athey and Imbens assert that synthetic control is “one of the most important development[s] in program evaluation in the past decade” and “arguably the most important innovation in the evaluation literature in the last fifteen years” athey. In a somewhat different direction, Hsiao et al. introduce the panel data method hsiao1, hsiao2, which seems to have a close bearing with some of the approaches of this work. In particular, hsiao1, hsiao2 only uses data for the outcome variable and solves an ordinary least squares problem in learning synthetic control. However, hsiao1, hsiao2 restrict the subset of possible controls to units that are within the geographical or economic proximity of the treated unit. Therefore, there is still some degree of subjectivity in the choice of the donor pool. In addition, hsiao1, hsiao2 do not include a “de-noising” step, which is a key feature of our approach. For an empirical comparison between the synthetic control and panel data methods, see comparison. It should be noted that comparison also adapts the panel data method to automate the donor selection process. doudchenko allows for an additive difference between the treated unit and donor pool, similar to the difference-in-differences (DID) method. Moreover, similar to our exposition, doudchenko relaxes the convexity aspect of synthetic control and proposes an algorithm that allows for unrestricted linearity as well as regularization. In an effort to infer the causal impact of market interventions, bayes_sc introduce yet another evaluation methodology based on a diffusion-regression state-space model that is fully Bayesian; similar to abadie1, abadie3, hsiao1, hsiao2, their model also generalizes the DID procedure. Due to the subjectivity in the choice of covariates and predictor variables, ferman1 provides recommendations for specification-searching opportunities in synthetic control applications. The recent work of xu extends the synthetic control method to allow for multiple treated units and variable treatment periods as well as the treatment being correlated with unobserved units. Similar to our work, xu computes uncertainty estimates; however, while xu obtains these measurements via a parametric bootstrap procedure, we obtain uncertainty estimates under a Bayesian framework.
Matrix completion and factorization approaches are well-studied problems with broad applications (e.g. recommendation systems, graphon estimation, etc.). As shown profusely in the literature, spectral methods, such as singular value decomposition and thresholding, provide a procedure to estimate the entries of a matrix from partial and/or noisy observations recht1. With our eyes set on achieving “robustness”, spectral methods become particularly appealing since they de-noise random effects and impute missing information within the data matrix svdjha1. For a detailed discussion on the topic, see usvt; for algorithmic implementations, see mazumder2010 and references there in. We note that our goal differs from traditional matrix completion applications in that we are using spectral methods to estimate a low-rank matrix, allowing us to determine a linear relationship between the rows of the mean matrix. This relationship is then projected into the future to determine the counterfactual evolution of a row in the matrix (treated unit), which is traditionally not the goal in matrix completion applications. Another line of work within this arena is to impute the missing entries via a nearest neighbor based estimation algorithm under a latent variable model framework lee_song, lee2.
There has been some recent work in using matrix norm methods in relation to causal inference, including for synthetic control. In athey1, the authors use matrix norm regularization techniques to estimate counterfactuals for panel data under settings that rely on the availability of a large number of units relative to the number of factors or characteristics, and under settings that involve limited number of units but plenty of history (synthetic control). This is different from our approach, which increases robustness by “de-noising” using matrix completion methods, and then using linear regression on the de-noised matrix, instead of relying on matrix norm regularizations.
Despite its popularity, there has been less theoretical work in establishing the consistency of the synthetic control method or its variants. abadie1 demonstrates that the bias of the synthetic control estimator can be bounded by a function that is close to zero when the pre-intervention period is large in relation to the scale of the transitory shocks, but under the additional condition that a perfect convex match between the pre-treatment noisy outcome and covariate variables for the treated unit and donor pool exists. ferman2 relaxes the assumption in abadie1, and derives conditions under which the synthetic control estimator is asymptotically unbiased under non-stationarity conditions. To our knowledge, however, no prior work has provided finite-sample analysis, analyzed the performance of these estimators with respect to the mean-squared error (MSE), established asymptotic consistency, or addressed the possibility of missing data, {a common handicap in practice}.
We will denote $\mathbb{R}$ as the field of real numbers. For any positive integer $N$, let $[N] = \{1, \dots, N\}$. For any vector $v \in \mathbb{R}^n$, we denote its Euclidean ($\ell_2$) norm by $\norm{v}_2$, and define $\norm{v}_2^2 = \sum_{i=1}^n v_i^2$. We define its infinity norm as $\norm{v}_{\infty} = \max_{i} \abs{v_i}$. In general, the $\ell_p$ norm for a vector $v$ is defined as $\norm{v}_p = \Big( \sum_{i=1}^n \abs{v_i}^p \Big)^{1/p}$. Similarly, for an $m \times n$ real-valued matrix $\boldsymbol{A} = [A_{ij}]$, its spectral/operator norm, denoted by $\norm{\boldsymbol{A}}_2$, is defined as $\norm{\boldsymbol{A}}_2 = \max_{1 \le i \le k} \abs{\sigma_i}$, where $k = \min\{m,n\}$ and $\sigma_i$ are the singular values of $\boldsymbol{A}$. The Moore-Penrose pseudoinverse $\boldsymbol{A}^{\dagger}$ of $\boldsymbol{A}$ is defined as
with $x_i$ and $y_i$ being the left and right singular vectors of $\boldsymbol{A}$, respectively. We will adopt the shorthand notation of $\norm{\cdot} \equiv \norm{\cdot}_2$. To avoid any confusions between scalars/vectors and matrices, we will represent all matrices in bold, e.g. $\boldsymbol{A}$.
Let $f$ and $g$ be two functions defined on the same space. We say that $f(x) = \mathcal{O}(g(x))$ and $f(x) = \Omega(g(x))$ if and only if there exists a positive real number $M$ and a real number $x_0$ such that for all $x \ge x_0$,
respectively.
The data at hand is a collection of time series with respect to an aggregated metric of interest (e.g. violent crime rates) comprised of both the treated unit and the donor pool outcomes. Suppose we observe $N \geq 2$ units across $T \geq 2$ time periods. We denote $T_0$ as the number of pre-intervention periods with $1 \leq T_0 < T$, rendering $T - T_0$ as the length of the post-intervention stage. Without loss of generality, let the first unit represent the treatment unit -- exposed to the intervention of interest at time $t = T_0 + 1$. The remaining donor units, $2\leq i\leq N$, are unaffected by the intervention for the entire time period $[T]=\{1,\dots, T\}$.
Let $X_{it}$ denote the measured value of metric for unit $i$ at time $t$. We posit
where $M_{it}$ is the deterministic mean while the random variables $\epsilon_{it}$ represent zero-mean noise that are independent across $i,t$. Following the philosophy of latent variable models usvt, lee_song, aldous, hoover1, hoover2, we further posit that for all $2 \le i \le N$, $t \in [T]$
where $\theta_i \in {\mathbb R}^{d_1}$ and $\rho_t \in {\mathbb R}^{d_2}$ are latent feature vectors capturing unit and time specific information, respectively, for some $d_1, d_2 \geq 1$; the latent function $f: {\mathbb R}^{d_1} \times {\mathbb R}^{d_2} \to {\mathbb R}$ captures the model relationship. We note that this formulation subsumes popular econometric factor models, such as the one presented in abadie1, as a special case with (small) constants $d_1 = d_2$ and $f$ as a bilinear function.
The treatment unit obeys the same model relationship during the pre-intervention period. That is, for $t \leq T_0$
where $M_{1t} = f(\theta_1, \rho_t)$ for some latent parameter $\theta_1 \in {\mathbb R}^{d_1}$. If unit one was never exposed to the intervention, then the same relationship as (ref) would continue to hold during the post-intervention period as well. In essence, we are assuming that the outcome random variables for {\em all} unaffected units follow the model relationship defined by (ref) and (ref). Therefore, the “synthetic control” would ideally help estimate the underlying counterfactual means $M_{1t} = f(\theta_1, \rho_t)$ for $T_0 < t \le T$ by using an appropriate combination of the post-intervention observations from the donor pool since the donor units are immune to the treatment.
To render this feasible, we make the key operating assumption (as done in literature cf. abadie1, abadie2, abadie3) that the mean vector of the treatment unit over the pre-intervention period, i.e. the vector $M_{1}^- = [M_{1t}]_{t\leq T_0}$, lies within the span of the mean vectors within the donor pool over the pre-intervention period, i.e. the span of the donor mean vectors $M_{i}^- = [M_{it}]_{2\le i \le N, t \le T_0}$ \footnote{We note that this is a minor departure from the literature on synthetic control starting in abadie3 -- in literature, the pre-intervention {\em noisy} observation (rather than the mean) vector $X_{1}$, is assumed to be a {\em convex} (rather than linear) combination of the noisy donor observations. We believe our setup is more reasonable since we do not want to fit noise.}. More precisely, we assume there exists a set of weights $\beta^* \in {\mathbb R}^{N-1}$ such that for all $t\leq T_0$,
This is a reasonable and intuitive assumption, utilized in literature, hypothesizing that the treatment unit can be modeled as some combination of the donor pool. In fact, the set of weights $\beta^*$ are the very definition of a synthetic control.
In order to distinguish the pre- and post-intervention periods, we use the following notation for all (donor) matrices: $\boldsymbol{A} = [\boldsymbol{A}^{-}, \boldsymbol{A}^{+}]$, where $\boldsymbol{A}^{-} = [A_{ij}]_{2 \le i \le N, j \in [T_0]}$ and $\boldsymbol{A}^{+} = [A_{ij}]_{2 \le i \le N, T_0 < j \le T}$ denote the pre- and post-intervention submatrices, respectively; vectors will be defined in the same manner, i.e. $A_i = [A_i^{-}, A_i^{+}]$, where $A_i^{-} = [A_{it}]_{t \in [T_0]}$ and $A_i^{+} = [A_{it}]_{T_0 < t \le T}$ denote the pre- and post-intervention subvectors, respectively, for the $i$th donor. Moreover, we will denote all vectors related to the treatment unit with the subscript “1”, e.g. $A_1 = [A_1^{-}, A_1^{+}]$.
In contrast with the classical synthetic control work, we allow our model to be robust to incomplete observations. To model randomly missing data, the algorithm observes each data point $X_{it}$ in the donor pool with probability $p \in (0, 1]$, independently of all other entries. While the assumption that $p$ is constant across all rows and columns of our observation matrix is standard in literature, our results remain valid even in situations where the probability of observation is dependent on the row and column latent parameters, i.e. $p_{ij} = g(\theta_i, \rho_j) \in (0, 1]$. In such situations, $p_{ij}$ can be estimated as $\hat{p}_{ij}$ using consistent graphon estimation techniques described in a growing body of literature, e.g. see lee2, usvt, wolfe1,yang1. These estimates can then be used in our analysis presented in Section (ref).
We begin by exploring the intuition behind our proposed two-step algorithm: (1) de-noising the data: since the singular values of our observation matrix, $\boldsymbol{X} = [X_{it}]_{2 \le i \le N, t \in [T]}$, encode both signal and noise, we aim to discover a low rank approximation of $\boldsymbol{X}$ that only incorporates the singular values associated with useful information; simultaneously, this procedure will naturally impute any missing observations. We note that this procedure is similar to the algorithm proposed in usvt. (2) learning $\beta^*$: using the pre-intervention portion of the de-noised matrix, we learn the linear relationship between the treatment unit and the donor pool prior to estimating the post-intervention counterfactual outcomes. Since our objective is to produce accurate predictions, it is not obvious why the synthetic treatment unit should be a convex combination of its donor pool as assumed in abadie1, abadie3, abadie4. In fact, one can reasonably expect that the treatment unit and some of the donor units may exhibit negative correlations with one another. In light of this intuition, we learn the optimal set of weights via linear regression, allowing for both positive and negative elements.
We present the details of our robust method in Algorithm (ref) below. The algorithm utilizes two hyperparameters: (1) a thresholding hyperparameter $\mu \geq 0$, which serves as a knob to effectively trade-off between the bias and variance of the estimator, and (2) a regularization hyperameter $\eta \ge 0$ that controls the model complexity. We discuss the procedure for determining the hyperparameters in Section (ref). To simplify the exposition, we assume the entries of $\boldsymbol{X}$ are bounded by one in absolute value, i.e. $ \abs{X_{it} } \leq 1$.
In order to quantitatively assess the uncertainty of our model, we will transition from a frequentist perspective to a Bayesian viewpoint. As commonly assumed in literature, we consider a zero-mean, isotropic Gaussian noise model (i.e. $\epsilon \mathbin{\overset{}{\kern\z@\sim}} \mathcal{N}(0, \sigma^2 \boldsymbol{I})$) and use the square loss for our cost function. We present the Bayesian method as Algorithm (ref). Note that we perform step one of our robust algorithm exactly as in Algorithm (ref); as a result, we only detail the alterations of step two in the Bayesian version (Algorithm (ref)).
Several of our results, as well as the algorithm we propose, assume that the observation matrix is bounded such that $\abs{X_{it}} \leq 1$. For any data matrix, we can achieve this by using the following pre-processing transformation: suppose the entries of $\boldsymbol{X}$ belong to an interval $[a,b]$. Then, one can first pre-process the matrix $\boldsymbol{X}$ by subtracting $(a+b)/2$ from each entry, and dividing by $(b-a)/2$ to enforce that the entries lie in the range $[-1, 1]$. The reverse transformation, which can be applied at the end of the algorithm description above, returns a matrix with values contained in the original range. Specifically, the reverse transformation equates to multiplying the end result by $(b-a)/2$ and adding by $(a + b)/2$.
For the practitioner who seeks a more interpretable solution, e.g. a convex combination of donors as per the original synthetic control estimator of Abadie et. al, we recommend using an $\ell_1$-regularization penalty in the learning procedure of step 2. Due to the geometry of LASSO, the resulting estimator will be often be a sparse vector. Specifically, for any $\eta > 0$, we define the LASSO estimator to be \[ \hat{\beta}(\eta) = \operatorname*{arg\,min}_{v \in \mathbb{R}^{N-1}} \norm{Y_1^{-} - (\hat{\boldsymbol{M}}^{-})^T v}^2 + \eta \norm{v}_1. \]
Here, we discuss several approaches to choosing the hyperparameter $\mu$ for the singular values. If it is known a priori that the underlying model is low rank with rank at most $k$, then it may make sense to choose $\mu$ such that $\abs{S} = k$. A data driven approach, however, could be implemented based on cross-validation. Precisely, reserve a portion of the pre-intervention period for validation, and use the rest of the pre-intervention data to produce an estimate $\hat{\beta}(\eta)$ for each of the finitely many choices of $\mu$ ($s_1, \dots, s_{N-1}$). Using each estimate $\hat{\beta}(\eta)$, produce its corresponding treatment unit mean vector over the validation period. Then, select the $\mu$ that achieves the minimum MSE with respect to the observed data. Finally, usvt provides a universal approach to picking a threshold; similarly, we also propose another such universal threshold, (ref), in Section (ref). We utilize the data driven approach in our experiments in this work.
The regularization parameter, $\eta$, also plays a crucial role in learning the synthetic control and influences both the training and generalization errors. As is often the case in model selection, a popular strategy in estimating the ideal $\eta$ is to employ cross-validation as described above. However, since time-series data often have a natural temporal ordering with causal effects, we also recommend employing the forward chaining strategy. Although the forward chaining strategy is similar to leave-one-out (LOO) cross-validation, an important distinction is that forward chaining does not break the temporal ordering in the training data. More specifically, for a particular candidate of $\eta$ at every iteration $t \in [T_0]$, the learning process uses $[Y_{11}, \dots, Y_{1,t-1}]$ as the training portion while reserving $Y_{1t}$ as the validation point. As before, the average error is then computed and used to evaluate the model (characterized by the choice of $\eta$). The forward chaining strategy can also be used to learn the optimal $\mu$.
In terms of scalability, the most computationally demanding procedure is that of evaluating the singular value decomposition (SVD) of the observation matrix. Given the ubiquity of SVD methods in the realm of machine learning, there are well-known techniques that enable computational and storage scaling for SVD algorithms. For instance, both Spark (through alternative least squares) and Tensor-Flow come with built-in SVD implementations. As a result, by utilizing the appropriate computational infrastructure, our de-noising procedure, and algorithm in generality, can scale quite well. Also note that for a low rank structure, we typically only need to compute the top few singular values and vectors. Various truncated-SVD algorithms provide resource-efficient implementations to compute the top $k$ singular values and vectors instead of the complete-SVD.
The factor models that are commonly used in the Econometrics literature, cf. abadie1, abadie2, abadie3, often lead to a low rank structure for the underlying mean matrix $\boldsymbol{M}$. When $f$ is nonlinear, $\boldsymbol{M}$ can still be well approximated by a low rank matrix for a large class of functions. For instance, if the latent parameters assumed values from a bounded, compact set, and if $f$ was Lipschitz continuous, then it can be argued that $\boldsymbol{M}$ is well approximated by a low rank matrix, cf. see usvt for a very simple proof. As the reader will notice, while we establish results for low rank matrix, the results of this work are robust to low rank approximations whereby the approximation error can be viewed as “noise”. Lastly, as shown in log_rank, many latent variable models can be well approximated (up to arbitrary accuracy $\epsilon$) by low rank matrices. Specifically, log_rank shows that the corresponding low rank approximation matrices associated with “nice” functions (e.g. linear functions, polynomials, kernels, etc.) are of log-rank.
Although the algorithm does not appear to rely on any helpful covariate information and the experimental results, presented in Section (ref), suggest that it performs on par with that of the original synthetic control algorithm, we want to emphasize that we are not suggesting that practitioners should abandon the use of any additional covariate information or the application of domain knowledge. Rather, we believe that our key algorithmic feature -- the de-noising step -- may render covariates and domain expertise as luxuries as opposed to necessities for many practical applications. If the practitioner has access to supplementary predictor variables, we propose that step one of our algorithm be used as a pre-processing routine for de-noising the data before incorporating additional information. Moreover, other than the obvious benefit of narrowing the donor pool, domain expertise can also come in handy in various settings, such as determining the appropriate method for imputing the missing entries in the data. For instance, if it is known a priori that there is a trend or periodicity in the time series evolution for the units, it may behoove the practitioner to impute the missing entries using “nearest-neighbors” or linear interpolation.
In this section, we derive the finite sample and asymptotic properties of the estimators $\hat{\boldsymbol{M}}$ and $\hat{M}_1$. We begin by defining necessary notations and recalling a few operating assumptions prior to presenting the results, with the corresponding proofs relegated to the Appendix. To that end, we re-write (ref) in matrix form as $\boldsymbol{X} = \boldsymbol{M} + \boldsymbol{E}$, where $\boldsymbol{E} = [\epsilon_{it}]_{2\leq i\leq N, t \in [T]}$ denotes the noise matrix. We shall assume that the noise parameters $\epsilon_{it}$ are independent zero-mean random variables with bounded second moments. Specifically, for all $2\leq i\leq N, t \in [T]$,
We shall also assume that the treatment unit noise in (ref) obeys (ref). Further, we assume the relationship in (ref) holds. To simplify the following exposition, we assume that $\abs{M_{ij}} \le 1$ and $\abs{X_{ij}} \le 1$.
As previously discussed, we evaluate the accuracy of our estimated means for the treatment unit with respect to the deviation between $\hat{M}_1$ and $M_1$ measured in $\ell_2$-norm, and similarly between $\hat{\boldsymbol{M}}$ and $\boldsymbol{M}$. Additionally, we aim to establish the validity of our pre-intervention linear model assumption (cf. (ref)) and investigate how the linear relationship translates over to the post-intervention regime, i.e. if $M_1^{-} = (\boldsymbol{M}^{-})^T \beta^*$ for some $\beta^*$, does $M^{+}_1$ (approximately) equal to $(\boldsymbol{M}^{+})^T \beta^*$? If so, under what conditions? We present our results for the above aspects after a brief motivation of $\ell_2$ regularization.
{\bf Combatting overfitting.} One weapon to combat overfitting is to constrain the learning algorithm to limit the effective model complexity by fitting the data under a simpler hypothesis. This technique is known as regularization, and it has been widely used in practice. To employ regularization, we introduce a complexity penalty term into the objective function (ref). For a general regularizer, the objective function takes the form
for some choice of positive constants $\eta$ and $q$. The first term measures the empirical error of the model on the given dataset, while the second term penalizes models that are too “complex” by controlling the “smoothness” of the model in order to avoid overfitting. In general, the impact/trade-off of regularization can be controlled by the value of the regularization parameter $\eta$. Via the use of Lagrange multipliers, we note that minimizing (ref) is equivalent to minimizing (ref) subject to the constraint that
for some appropriate value of $c$. When $q = 2$, (ref) corresponds to the classical setup known as ridge regression or weight decay. The case of $q =1$ is known as the LASSO in the statistics literature; the $\ell_1$-norm regularization of LASSO is a popular heuristic for finding a sparse solution. In either case, incorporating an additional regularization term encourages the learning algorithm to output a simpler model with respect to some measure of complexity, which helps the algorithm avoid overfitting to the idiosyncrasies within the observed dataset. Although the training error may suffer from the simpler model, empirical studies have demonstrated that the generalization error can be greatly improved under this new setting. Throughout this section, we will primarily focus our attention on the case of $q = 2$, which maintains our learning objective to be (convex) quadratic in the parameter $v$ so that its exact minimizer can be found in closed form:
In this section, we highlight the importance of our de-noising procedure and prescribe a universal threshold (similar to that of usvt) that dexterously distinguishes signal from noise, enabling the algorithm to capture the appropriate amount of useful information (encoded in the singular values of $\boldsymbol{Y}$) while discarding out the randomness. Due to its universality, the threshold naturally adapts to the amount of structure within $\boldsymbol{M}$ in a purely data-driven manner. Specifically, for any choice of $\omega \in (0.1,1)$, we find that choosing
results in an estimator with strong theoretical properties for both interpolation and extrapolation (discussed in Section (ref)). Here, $\hat{p}$ and $\hat{\sigma}^2$ denote the unbiased maximum likelihood estimates of $p$ and $\sigma^2$, respectively, and can be computed via (ref) and (ref).
The following Theorems (adapted from Theorems 2.1 and 2.7 of usvt) demonstrate that Step 1 of our algorithm (detailed in Section (ref)) accurately imputes missing entries within our data matrix $\boldsymbol{X}$ when the signal matrix $\boldsymbol{M}$ is either low rank or generated by an $\mathcal{L}$-Lipschitz function. In particular, Theorems (ref) and (ref) states that Step 1 produces a consistent estimator of the underlying mean matrix $\boldsymbol{M}$ with respect to the (matrix) mean-squared-error, which is defined as
We say that $\hat{\boldsymbol{M}}$ is a consistent estimator of $\boldsymbol{M}$ if the right-hand side of (ref) converges to zero as $N$ and $T$ grow without bound.
The following theorem demonstrates that $\hat{\boldsymbol{M}}$ is a good estimate of $\boldsymbol{M}$ when $\boldsymbol{M}$ is a low rank matrix, particularly when the rank of $\boldsymbol{M}$ is small compared to $(N-1)p$. \\
Suppose that the latent row and column feature vectors, $\{\theta_i\}$ and $\{\rho_j\}$, belong to some bounded, closed sets $K \subset \mathbb{R}^d$, where $d$ is some arbitrary but fixed dimension. If we assume $f: K \times K \rightarrow [-1, 1]$ possesses desirable smoothness properties such as Lipschitzness, then $\hat{\boldsymbol{M}}$ is again a good estimate of $\boldsymbol{M}$. \\
It is important to observe that the models under consideration for both Theorems (ref) and (ref) encompass the mean matrices, $\boldsymbol{M}$, generated as per many of the popular Econometric factor models often considered in literature and assumed in practice. Therefore, de-noising the data serves as an important imputing and filtering procedure for a wide array of applications.
Similar to the setting for interpolation, the prediction performance metric of interest is the average mean-squared-error in estimating $M^-_1$ using $\hat{M}^-_1$. Precisely, we define
If the right-hand side of (ref) approaches zero in the limit as $T_0$ grows without bound, then we say that $\hat{M}_1^-$ is a consistent estimator of $M_1^-$ (note that our analysis here assumes that only $T_0 \rightarrow \infty$).
In what follows, we first state the finite sample bound on the average MSE between $\hat{M}_1^-$ and $M_1^-$ for the most generic setup (Theorem (ref)). As a main Corollary of the result, we specialize the bound in the case where we use our prescribed universal threshold. Finally, we discuss a minor variation of the algorithm where the data is pre-processed, and specialize the above result to establish the consistency of our estimator (Theorem (ref)).
We provide a finite sample error bound for the most generic setting, i.e. for any choice of the threshold, $\mu$, and regularization hyperparameter, $\eta$. \\
{\bf Bias-variance tradeoff.} Let us interpret the result by parsing the terms in the error bound. The last term decays exponentially with $(N-1)T$, as long as the fraction of observed entries is such that, on average, we see a super-constant number of entries, i.e. $p (N-1) T \gg 1$. More interestingly, the first two terms highlight the “bias-variance tradeoff” of the algorithm with respect to the singular value threshold $\mu$. Precisely, the size of the set $S$ increases with a decreasing value of the hyperparameter $\mu$, causing the second error term to increase. Simultaneously, however, this leads to a decrease in $\lambda^*$. Note that $\lambda^*$ denotes the aspect of the “signal” within the matrix $\boldsymbol{M}$ that is not captured due to the thresholding through $S$. On the other hand, the second term, $|S| \sigma^2/ T_0$, represents the amount of “noise” captured by the algorithm, but wrongfully interpreted as a signal, during the thresholding process. In other words, if we use a large threshold, then our model may fail to capture pertinent information encoded in $\boldsymbol{M}$; if we use a small threshold, then the algorithm may overfit the spurious patterns in the data. Thus, the hyperparameter $\mu$ provides a way to trade-off “bias” (first term) and “variance” (second term).
Using the universal threshold defined in (ref), we now highlight the prediction power of our estimator for any choice of $\eta$, the regularization hyperparameter. As described in Section (ref), the prescribed threshold automatically captures the “correct” level of information encoded in the (noisy) singular values of $\boldsymbol{Y}$ in a data-driven manner, dependent on the structure of $\boldsymbol{M}$. However, unlike the statements in Theorems (ref) and (ref), the following bound does not require $\boldsymbol{M}$ to be low rank or $f$ to be Lipschitz. \\
As an implication, if $p = (1 + \vartheta) \sqrt{T_0} / (1 + \sqrt{T_0})$ and $\sigma^2 \le \vartheta$, we have that $\text{MSE}(\hat{M}^-_1) = \mathcal{O}(1/\sqrt{T_0})$. More generally, Corollary (ref) shows that by adroitly capturing the signal, the resulting error bound simply depends on the variance of the noise terms, $\sigma^2$, and the error introduced due to missing data. Ideally, one would hope to overcome the error term when $T_0$ is sufficiently large. This motivates the following setup.
We present a straightforward pre-processing step that leads to the consistency of our algorithm. The pre-processing step simply involves replacing the columns of $\boldsymbol{X}$ by the averages of subsets of its columns. This admits the same setup as before, but with the variance for each noise term reduced. An implicit side benefit of this approach is that required SVD step in the algorithm is now applied to a matrix of smaller dimensions.
To begin, partition the $T_0$ columns of the pre-intervention data matrix $\boldsymbol{X}^{-}$ into $\Delta$ blocks, each of size $\tau = \lfloor T_0 / \Delta \rfloor$ except potentially the last block, which we shall ignore for theoretical purposes; in practice, however, the remaining columns can be placed into the last block. Let $B_j = \{ (j-1) \tau + \ell : 1\leq \ell \le \tau \}$ denote the column indices of $\boldsymbol{X}^{-}$ within partition $j \in [\Delta]$. Next, we replace the $\tau$ columns within each partition by their average, and thus create a new matrix, $\bar{\boldsymbol{X}}^{-}$, with $\Delta$ columns and $N-1$ rows. Precisely, $\bar{\boldsymbol{X}}^{-} = [\bar{X}_{ij}]_{2\leq i\leq N, j \in [\Delta]}$ with
where
For the treatment row, let $\bar{X}_{1j} = \frac{\hat{p}}{\tau} \sum_{t \in B_j} X_{1t}$ for all $j \in [\Delta]$\footnote{Although the statement in Theorem (ref) assumes that an oracle provides the true $p$, we prescribe practitioners to use $\hat{p}$ since $\hat{p}$ converges to $p$ almost surely by the Strong Law of Large Numbers.}. Let $\bar{\boldsymbol{M}}^{-} = [\bar{M}_{ij}]_{2\leq i\leq N, j \in [\Delta]}$ with
We apply the algorithm to $\bar{\boldsymbol{X}}^{-}$ to produce the estimate $\hat{\bar{\boldsymbol{M}}}^{-}$ of $\bar{\boldsymbol{M}}^{-}$, which is sufficient to produce $\hat{\beta}(\eta)$. This $\hat{\beta}(\eta)$ can be used to produce the post-intervention synthetic control means $\hat{M}_{1}^{+} = [\hat{M}_{1t}]_{T_0 < t \leq T}$ in a similar manner as before \footnote{In practice, one can first de-noise $\boldsymbol{X}^+$ via step one of Section (ref), and use the entries of $\hat{\boldsymbol{M}}^+$ in (ref).}: for $T_0 < t \leq T$,
For the pre-intervention period, we produce the estimator $\hat{\bar{M}}_1^{-} = [\hat{\bar{M}}_{1j}]_{j \in [\Delta]}$: for $j \in [\Delta]$,
Our measure of estimation error is defined as
For simplicity, we will analyze the case where each block contains at least one entry such that $\bar{\boldsymbol{X}}^-$ is completely observed. We now state the following result. \\
We note that the method of abadie3 learns the weights (here $\hat{\beta}(0)$) by pre-processing the data. One common pre-processing proposal is to also aggregate the columns, but the aggregation parameters are chosen by solving an optimization problem to minimize the resulting prediction error of the observations. In that sense, the above averaging of column is a simple, data agnostic approach to achieve a similar effect, and potentially more effectively.
For the post-intervention regime, we consider the average root-mean-squared-error in measuring the performance of our algorithm. Precisely, we define
The key assumption of our analysis is that the treatment unit signal can be written as a linear combination of donor pool signals. Specifically, we assume that this relationship holds in the pre-intervention regime, i.e. $M_1^{-} = (\boldsymbol{M}^{-})^T\beta^*$ for some $\beta^* \in {\mathbb R}^{N-1}$ as stated in (ref). However, the question still remains: does the same relationship hold for the post-intervention regime and if so, under what conditions does it hold? We state a simple linear algebraic fact to this effect, justifying the approach of synthetic control. It is worth noting that this important aspect has been amiss in the literature, potentially implicitly believed or assumed starting in the work by abadie3. \\
If we assume that the linear relationship prevails in the post-intervention period, then we arrive at the following error bound. \\
Benefits of regularization. In order to motivate the use of regularization, we analyze the error bounds of Theorems (ref) and (ref) to observe how the pre- and post-intervention errors react to regularization. As seen from Theorem (ref), the pre-intervention error {\em increases} linearly with respect to the choice of $\eta$. Intuitively, this increase in pre-intervention error derives from the fact that regularization reduces the model complexity, which biases the model and handicaps its ability to fit the data. At the same time, by restricting the hypothesis space and controlling the “smoothness” of the model, regularization prevents the model from overfitting to the data, which better equips the model to generalize to unseen data. Therefore, a larger value of $\eta$ {\em reduces} the post-intervention error. This can be seen by observing the second error term of Theorem (ref), which is controlled by the expression $\norm{\hat{\beta}(\eta) - \beta^*}$. In words, this error is a function of the learning algorithm used to estimate $\beta^*$. Interestingly, farebrother demonstrates that there exists an $\eta > 0$ such that
without any assumptions on the rank of $\hat{\boldsymbol{M}}^{-}$. In other words, farebrother demonstrates that regularization can decrease the MSE between $\hat{\beta}(\eta)$ and the true $\beta^*$, thus reducing the overall error. Ultimately, employing ridge regression introduces extraneous bias into our model, yielding a higher pre-intervention error. In exchange, regularization reduces the post-intervention error (due to smaller variance).
We turn our attention to a Bayesian treatment of synthetic control. By operating under a Bayesian framework, we allow practitioners to naturally encode domain knowledge into prior distributions while simultaneously avoiding the problem of overfitting. In addition, rather than making point estimates, we can now quantitatively express the uncertainty in our estimates with posterior probability distributions.
We begin by treating $\beta^*$ as a random variable as opposed to an unknown constant. In this approach, we specify a prior distribution, $p(\beta)$, that expresses our apriori beliefs and preferences about the underlying parameter (synthetic control). Given some new observation for the donor units, our goal is to make predictions for the counterfactual treatment unit on the basis of a set of pre-intervention (training) data. For the moment, let us assume that the noise parameter $\sigma^2$ is a known quantity and that the noise is drawn from a Gaussian distribution with zero-mean; similarly, we temporarily assume $\boldsymbol{M}^{-}$ is also given. Let us denote the vector of donor estimates as $M_{\cdot t} = [M_{it}]_{2 \le i \le N}$; we define $X_{\cdot t}$ similarly. Denoting the pre-intervention data as $D = \{ (Y_{1t}, M_{\cdot t}): t \in [T_0] \}$, the likelihood function $p(Y_1^{-} \given \beta, \boldsymbol{M}^{-})$ is expressed as
an exponential of a quadratic function of $\beta$. The corresponding conjugate prior, $p(\beta)$, is therefore given by a Gaussian distribution, i.e. $\beta \mathbin{\overset{}{\kern\z@\sim}} \mathcal{N}( \beta \given \beta_0, \boldsymbol{\Sigma}_0)$ with mean $\beta_0$ and covariance $\Sigma_0$. By using a conjugate Gaussian prior, the posterior distribution, which is proportional to the product of the likelihood and the prior, will also be Gaussian. Applying Bayes' Theorem (derivation unveiled in the Appendix), we have that the posterior distribution is $p(\beta \given D) = \mathcal{N}( \beta_D, \boldsymbol{\Sigma}_D)$ where
For the remainder of this section, we shall consider a popular form of the Gaussian prior. In particular, we consider a zero-mean isotropic Gaussian with the following parameters: $\beta_0 = 0$ and $\boldsymbol{\Sigma}_0 = \alpha^{-1} \boldsymbol{I}$ for some choice of $\alpha > 0$. Since $\boldsymbol{M}^{-}$ is unobserved by the algorithm, we use the estimated $\hat{\boldsymbol{M}}^{-}$, computed as per step one of Section (ref), as a proxy; therefore, we redefine our data as $D = \{ (Y_{1t}, \hat{M}_{\cdot t}): t \in [T_0] \}$. Putting everything together, we have that $p(\beta \given D) = \mathcal{N}( \beta_D, \boldsymbol{\Sigma}_D)$ whereby
By using the zero-mean, isotropic Gaussian conjugate prior, we can derive a point estimate of $\beta^*$ by maximizing the log posterior distribution, which we will show is equivalent to minimizing the regularized objective function of (ref) for a particular choice of $\eta$. In essence, we are determining the optimal $\hat{\beta}$ by finding the most probable value of $\beta^*$ given the data and under the influence of our prior beliefs. The resulting estimate is known as the maximum a posteriori (MAP) estimate.
We begin by taking the log of the posterior distribution, which gives the form
Maximizing the above log posterior then equates to minimizing the quadratic regularized error (ref) with $\eta = \alpha \sigma^2$. We define the MAP estimate, $\hat{\beta}_{\text{MAP}}$, as
With the MAP estimate at hand, we then make predictions of the counterfactual as
Therefore, we have seen that the MAP estimation is equivalent to ridge regression since the introduction of an appropriate prior naturally induces the additional complexity penalty term.
Although we have treated $\beta^*$ as a random variable attached with a prior distribution, we can venture beyond point estimates to be fully Bayesian. In particular, we will make use of the posterior distribution over $\beta^*$ to marginalize over all possible values of $\beta^*$ in evaluating the predictive distribution over $Y_1^{-}$. We will decompose the regression problem of predicting the counterfactual into two separate stages: the inference stage in which we use the pre-intervention data to learn the predictive distribution (defined shortly), and the subsequent decision stage in which we use the predictive distribution to make estimates. By separating the inference and decision stages, we can readily develop new estimators for different loss functions without having to relearn the predictive distribution, providing practitioners tremendous flexibility with respect to decision making.
Let us begin with a study of the inference stage. We evaluate the predictive distribution over $Y_{1t}$, which is defined as
where
Note that $p(\beta \given D)$ is the posterior distribution over the synthetic control parameter and is governed by (ref) and (ref). With access to the predictive distribution, we move on towards the decision stage, which consists of determining a particular estimate $\hat{M}_{1t}$ given a new observation vector $X_{\cdot t}$ (used to determine $\hat{M}_{\cdot t}$). Consider an arbitrary loss function $L(Y_{1t}, g(\hat{M}_{\cdot t}))$ for some function $g$. The expected loss is then given by
and we choose our estimator $\hat{g}(\cdot)$ as the function that minimizes the average cost, i.e.,
Since $p(\hat{M}_{\cdot t}) \ge 0$, we can minimize (ref) by selecting $\hat{g}(\hat{M}_{\cdot t})$ to minimize the term within the parenthesis for each individual value of $Y_{1t}$, i.e.,
As suggested by (ref), the optimal estimate $\hat{M}_{1t}$ for a particular loss function depends on the model only through the predictive distribution $p(Y_{1t} \given \hat{M}_{\cdot t}, D)$. Therefore, the predictive distribution summarizes all of the necessary information to construct the desired Bayesian estimator for any given loss function $L$.
We analyze the case for the squared loss function (MSE), a common cost criterion for regression problems. In this case, we write the expected loss as
Under the MSE cost criterion, the optimal estimate is the mean of the predictive distribution, also known as the Bayes' least-squares (BLS) estimate:
We begin by exploring two real-world case studies discussed in abadie1, abadie2, abadie3 that demonstrate the ability of the original synthetic control's algorithm to produce a reliable counterfactual reality. We use the same case-studies to showcase the “robustness” property of our proposed algorithm. Specifically, we demonstrate that our algorithm reproduces similar results even in presence of missing data, and without knowledge of the extra covariates utilized by prior works. We find that our approach, surprisingly, also discovers a few subtle effects that seem to have been overlooked in prior studies. In the following empirical studies, we will employ three different learning procedures as described in the robust synthetic control algorithm: (1) linear regression ($\eta = 0$), (2) ridge regression ($\eta > 0$), and (3) LASSO ($\zeta > 0$).
As described in abadie1, abadie2, abadie4, the synthetic control method allows a practitioner to evaluate the reliability of his or her case study results by running placebo tests. One such placebo test is to apply the synthetic control method to a donor unit. Since the control units within the donor pool are assumed to be unaffected by the intervention of interest (or at least much less affected in comparison), one would expect that the estimated effects of intervention for the placebo unit should be less drastic and divergent compared to that of the treated unit. Ideally, the counterfactuals for the placebo units would show negligible effects of intervention. Similarly, one can also perform exact inferential techniques that are similar to permutation tests. This can be done by applying the synthetic control method to every control unit within the donor pool and analyzing the gaps for every simulation, and thus providing a distribution of estimated gaps. In that spirit, we present the resulting placebo tests (for only the case of linear regression) for the Basque Country and California Prop. 99 case studies below to assess the significance of our estimates.
We will also analyze both case studies under a Bayesian setting. From our results, we see that our predictive uncertainty, captured by the standard deviation of the predictive distribution, is influenced by the number of singular values used in the de-noising process. Therefore, we have plotted the eigenspectrum of the two case study datasets below. Clearly, the bulk of the signal contained within the datasets is encoded into the top few singular values -- in particular, the top two singular values. Given that the validation errors computed via forward chaining are nearly identical for low-rank settings (with the exception of a rank-1 approximation), we shall use a rank-2 approximation of the data matrix. In order to exhibit the role of thresholding in the interplay between bias and variance, we also plot the cases where we use threshold values that are too high (bias) or too low (variance). For each figure, the dotted blue line will represent our posterior predictive means while the shaded light blue region spans one standard deviation on both sides of the mean. As we shall see, our predictive uncertainty is smallest in the neighborhood around the pre-intervention period. However, the level of uncertainty increases as we deviate from the the intervention point, which appeals to our intuition.
In order to choose an appropriate choice of the prior parameter $\alpha$, we first use forward-chaining for the ridge regression setting to find the optimal regularization hyperparameter $\eta$. By observing the expressions of (ref) and (ref), we see that $\eta = \alpha \sigma^2$ since ridge regression is closely related to MAP estimation for a zero-mean, isotropic Gaussian prior. Consequently, we choose $\alpha = \eta / \hat{\sigma}^2$ where $\eta$ is the value obtained via forward chaining.
The goal of this case-study is to investigate the effects of terrorism on the economy of Basque Country using the neighboring Spanish regions as the control group. In 1968, the first Basque Country victim of terrorism was claimed; however, it was not until the mid-1970s did the terrorist activity become more rampant abadie3. To study the economic ramifications of terrorism on Basque Country, we only use as data the per-capita GDP (outcome variable) of 17 Spanish regions from 1955-1997. We note that in abadie3, 13 additional predictor variables for each region were used including demographic information pertaining to one's educational status, and average shares for six industrial sectors.
Results. Figure (ref) shows that our method (for all three estimators) produces a very similar qualitative synthetic control to the original method even though we do not utilize additional predictor variables. Specifically, the synthetic control resembles the observed GDP in the pre-treatment period between 1955-1970. However, due to the large-scale terrorist activity in the mid-70s, there is a noticeable economic divergence between the synthetic and observed trajectories beginning around 1975. This deviation suggests that terrorist activity negatively impacted the economic growth of Basque Country.
One subtle difference between our synthetic control -- for the case of linear and ridge regression -- and that of abadie3 is between 1970-75: our approach suggests that there was a small, but noticeable economic impact starting just prior to 1970, potentially due to first terrorist attack in 1968. Notice, however, that the original synthetic control of abadie3 diverges only after 1975. Our LASSO estimator's trajectory also agrees with that of the original synthetic control method's, which is intuitive since both estimators seek sparse solutions.
To study the robustness of our approach with respect to missing entries, we discard each data point uniformly at random with probability $1-p$. The resulting control for different values of $p$ is presented in Figure (ref) suggesting the robustness of our (linear) algorithm. Finally, we produce Figure (ref) by applying our algorithm without the de-noising step. As evident from the Figure, the resulting predictions suffer drastically, reinforcing the value of de-noising. Intuitively, using an appropriate threshold $\mu$ equates to selecting the correct model complexity, which helps safeguard the algorithm from potentially overfitting to the training data.
Placebo tests. We begin by applying our robust algorithm to the Spanish region of Cataluna, a control unit that is not only similar to Basque Country, but also exposed to a much lower level of terrorism abadie2. Observing both the synthetic and observed economic evolutions of Cataluna in Figure (ref), we see that there is no identifiable treatment effect, especially compared to the divergence between the synthetic and observed Basque trajectories. We provide the results for the regions of Aragon and Castilla Y Leon in Figures (ref) and (ref).
Finally, similar to abadie2, we plot the differences between our estimates and the observations for Basque Country and all other regionals, individually, as placebos. Note that abadie2 excluded five regions that had a poor pre-intervention fit but we keep all regions. Figure (ref) shows the resulting plot for all regions with the solid black line being Basque Country. This plot helps visualize the extreme post-intervention divergence between the predicted means and the observed values for Basque. Up until about 1990, the divergence for Basque Country is the most extreme compared to all other regions (placebo studies) lending credence to the belief that the effects of terrorism on Basque Country were indeed significant. Refer to Figure (ref) for the same test but with Madrid and Balearic Islands excluded, as per abadie2. The conclusions drawn should remain the same, pointing to the robustness of our approach.
Bayesian approach. We plot the resulting Bayesian estimates in the figures below under varying thresholding conditions. It is interesting to note that our uncertainty grows dramatically once we include more than two singular values in the thresholding process. This confirms what our theoretical results indicated earlier: choosing a smaller threshold, $\mu$, would lead to a greater number of singular values retained which results in higher variance. On the other hand, notice that just selecting 1 singular value results in an apparently biased estimate which is overestimating the synthetic control. It appears that selecting the top two singular values balance the bias-variance tradeoff the best and is also agrees with our earlier finding that the data matrix appears to be of rank 2 or 3. Note that in this setting, we would find it hard to reject the null-hypothesis because the observations for the treated unit lie within the uncertainty band of the estimated synthetic control.
We study the impact of California's anti-tobacco legislation, Proposition 99, on the per-capita cigarette consumption of California. In 1988, California introduced the first modern-time large-scale anti-tobacco legislation in the United States abadie1. To analyze the effect of California's anti-tobacco legislation, we use the annual per-capita cigarette consumption at the state-level for all 50 states in the United States, as well as the District of Columbia, from 1970-2015. Similar to the previous case study, abadie3 uses 6 additional observable covariates per state, e.g. retail price, beer consumption per capita, and percentage of individuals between ages of 15-24, to predict their synthetic California. Furthermore, abadie3 discarded 12 states from the donor pool since some of these states also adopted anti-tobacco legislation programs or raised their state cigarette taxes, and discarded data after the year 2000 since many of the control units had implemented anti-tobacco measures by this point in time.
Results. As shown in Figure (ref), in the pre-intervention period of 1970-88, our control matches the observed trajectory. Post 1988, however, there is a significant divergence suggesting that the passage of Prop. 99 helped reduce cigarette consumption. Similar to the Basque case-study, our estimated effect is similar to that of abadie3. As seen in Figure (ref), our algorithm is again robust to randomly missing data.
Placebo tests. We now proceed to apply the same placebo tests to the California Prop 99 dataset. Figures (ref), (ref), and (ref) are three examples of the applied placebo tests on the remaining states (including District of Columbia) within the United States.
Finally, similar to abadie1, we plot the differences between our estimates and the actual observations for California and all other states, individually, as placebos. Note that abadie1 excluded twelve states but we keep all states. Figure (ref) shows the resulting plot for all states with the solid black line being California. This plot helps visualize the extreme post-intervention divergence between the predicted means and the observed values for California. Up until about 1995, the divergence for California was clearly the most significant and consistent outlier compared to all other regions (placebo studies) lending credence to the belief that the effects of Proposition 99 were indeed significant. Refer to Figure (ref) for the same test but with the same twelve states excludes as in abadie1. Just like the Basque Country case study, the exclusion of states should not affect the conclusions drawn.
Bayesian approach. Similar to the Basque Country case study, our predictive uncertainty increases as the number of singular values used in the learning process exceeds two. In order to gain some new insight, however, we will focus our attention to the resulting figure associated with three singular values, which is particularly interesting. Specifically, we observe that our predictive means closely match the counterfactual trajectory produced by the classical synthetic control method in both the pre- and post-intervention periods (up to year 2000), and yet our uncertainty for this estimate is significantly greater than our uncertainty associated with the estimate produced using two singular values. As a result, it may be possible that the classical synthetic control method overestimated the effect of Prop. 99, even though the legislation did probably discourage the consumption of cigarettes -- a conclusion reached by both our robust approach and the classical approach.
We conduct synthetic simulations to establish the various properties of the estimates in both the pre- and post-intervention stages.
Experimental setup. For each unit $i \in [N]$, we assign latent feature $\theta_i$ by drawing a number uniformly at random in $[0,1]$. For each time $t \in [T]$, we assign latent variable $\rho_t = t$. The mean value $m_{it} = f(\theta_i, \rho_t)$. In the experiments described in this section, we use the following:
where $f_1, f_2, f_3, f_4$ define the periodicities: $f_1 = \rho_t \bmod(360), f_2 = \rho_t \bmod(180), f_3 = 2\cdot \rho_t \bmod(360), f_4 = 2.0 \cdot \rho_t \bmod(180)$. The observed value $X_{it}$ is produced by adding i.i.d. Gaussian noise to mean with zero mean and variance $\sigma^2$. For this set of experiments, we use $N = 100, T = 2000$, while assuming the treatment was performed at $t = 1600$.
Training error approximates generalization error. For the first experimental study, we analyze the relationship between the pre-intervention MSE (training error) and the post-intervention MSE (generalization error). As seen in Table (ref), the post-intervention MSE closely matches that of the pre-intervention MSE for varying noise levels, $\sigma^2$. Thus suggesting efficacy of our algorithm. Figures (ref) and (ref) plot the estimates of algorithm with no missing data (Figure (ref)) and with 40% randomly missing data (Figure (ref)) on the same underlying dataset. All entries in the plots were normalized to lie within $[-1, 1]$. These plots confirm the robustness of our algorithm. Our algorithm outperforms the algorithm with no singular value thresholding under all proportions of missing data. The estimates from the algorithm which performs no singular value thresholding (green) degrade significantly with missing data while our algorithm remains robust.
Benefits of de-noising. We now analyze the benefit of de-noising the data matrix, which is the main contribution of this work compared to the prior work. Specifically, we study the generalization error of method using de-noising via thresholding and without thresholding as in prior work. The results summarized in Table (ref) show that for range of parameters the generalization error with de-noising is consistency better than that without de-noising.
Bayesian approach. From the synthetic simulations (figures below), we see that the number of singular values included in the thresholding process plays a crucial role in the model's prediction capabilities. If not enough singular values are used, then there is a significant loss of information (high bias) resulting in a higher MSE. On the other hand, if we include too many singular values, then the model begins to overfit to the dataset by misinterpreting noise for signal (high variance). As emphasized before, the goal is to find the simplest model that both fits the data and is also plausible, which is achieved when four singular values are employed.
The classical synthetic control method is recognized as a powerful and effective technique for causal inference for comparative case studies. In this work, we motivate a robust synthetic control algorithm, which attempts to improve on the classical method in the following regimes: (a) randomly missing data and (b) large levels of noise. We also demonstrate that the algorithm performs well even in the absence of covariate or expert information, but do not propose ignoring information which may eliminate “bad” donors. Our data-driven algorithm, and its Bayesian counterpart, uses singular value thresholding to impute missing data and “de-noise” the observations. Once “de-noised”, we use regularized linear regression to determine the synthetic control. Motivating our algorithm is a modeling framework, specifically the Latent Variable Model, which is a generalization of the various factor models used in related work. We establish finite-sample bounds on the MSE between the estimated “synthetic” control and the latent true means of the treated unit of interest. In situations with plentiful data, we show that a simple data aggregation method can lead to an asymptotically consistent estimator. Experiments on synthetically generated data (where the truth is known) and on real-world case-studies allow us to demonstrate the promise of our approach, which is an improvement over the classical method.