Bayesian Synthetic Control with a Soft Simplex Constraint

Bayesian Synthetic Control with a Soft Simplex Constraint

\def\spacingset#1{ {#1}} \spacingset{1}

\spacingset{1.2}

{\it Keywords:} average treatment effect; constrained linear regression; Metropolis-within-Gibbs sampler; spike-and-slab prior; strong selection consistency.

Introduction

Background on Synthetic Control Methods

In many fields of social science, researchers need to identify the causal effect of an intervention or treatment of interest using only observational data, since conducting case-control experiments is either infeasible or improper. The most commonly used causal inference method in economics is the Difference-in-Differences approach. Its validity, however, relies heavily on the parallel trends assumption: the control unit must be sufficiently similar to the treated unit so that it can mimic the pre-intervention dynamics of the treated. In many empirical settings, this assumption is overly restrictive, especially when the outcome of interest is measured at an aggregate level, such as a country's GDP. The synthetic control method (SCM) was developed to overcome this limitation, which constructs a counterfactual outcome of the treated unit using a weighted average of the control units that closely tracks the pre-treatment trajectory of the treated. For example, in the seminal work of abadie2003economic, SCM was developed to study how terrorist conflicts affected the local economy in the Basque Country between the 1970s and 1990s.

Formally, SCM can be described as follows. Consider a panel with \( T = T_0 + T_1 \) total observations over time for \( N+1 \) units, where only the first unit receives the treatment after period $T_0$. The remaining \( N \) units form the control group, which are unaffected by the intervention. Let \( Z_{i, t} \) denote the outcome of interest for unit \( i \) at time \( t \), where \( i = 1, 2, ..., N+1 \) and \( t = 1, 2, ..., T \). The objective of SCM is to estimate the average treatment effect on the treated (ATT) by estimating the counterfactual outcome \( Z_{1, t}^{(0)} \) for the treated unit at times $t = T_0 + 1, \dots, T$, which represents what the outcome of unit 1 would have been in the absence of the treatment. The standard approach is to define \( Z_{1, t}^{(0)} \) as a weighted combination of the outcomes of the \( N \) untreated units:

where the weights $(w_i)_{i=1}^{N}$ are time-independent, non-negative and sum to one. A larger \( w_i \) indicates that the $i$-th control unit better approximates the treated unit's pre-intervention characteristics. To find the optimal set of weights $(w_i)_{i=1}^{N}$ such that the synthetic control approximates the pre-treatment outcome of the treated unit as closely as possible, a common strategy is to minimize the mean squared error between the treated unit and the synthetic control in the pre-treatment period \(t = 1, \dots, T_0 \):

For ease of presentation, we do not consider time-invariant characteristics in the above formulation, and readers are referred to abadie2021syntheticcontrol for a more general description of SCM. The average difference between the observed outcome \( Z_{1, t} \) and the synthetic control \( Z_{1, t}^{(0)}(\hat{w}) \) serves as an estimate of the ATT:

Since its introduction, SCM has gained widespread adoption in empirical research, and athey2017state described SCM as “arguably the most important innovation in the policy evaluation literature in the last 15 years.” However, the growing application of SCM has also revealed new challenges. The first is the high dimensionality: when the number of the control units is comparable to, or even larger than, the number of pre-treatment periods, conventional least-squares-based methods are no longer applicable. To address this issue, researchers have proposed alternatives for computing the weights in (ref), such as penalized regression methods doudchenko2016balancing, carvalho2018arco, factor-model-based approaches xu_2017, pang2022bayesian, or variable selection kim2020bayesian, shi2023forward. However, none of these methods utilizes the simplex constraint and thus the resulting estimator $\hat{w}$ loses its interpretation as weights. How to implement SCM in a high-dimensional setting while preserving the simplex constraint remains largely unclear.

This naturally leads to another major issue with SCM: whether the simplex constraint should be imposed. Proponents argue that the simplex constraint prevents extrapolation king2006dangers, abadie2021syntheticcontrol. Critics, however, argue that the simplex constraint makes the extension of SCM to high-dimensional settings complicated and is often too restrictive in practice. A well-known example is the study of the economic impact of Hong Kong's reunification with China in 1997 hsiao2012panel. Since no other Asian regions faced a similar “handover” situation as Hong Kong, it is natural to apply SCM to estimate the impact. However, since Hong Kong was long known for its outstanding economic performance at that time (known as the leading region in the “Four Asian Tigers”), researchers found it unreasonable to use a convex combination of other regions to track the economic trajectory of Hong Kong. Allowing the sum of the coefficients to exceed one is in fact consistent with Hong Kong’s historically strong economic performance.

In this work, we develop a new Bayesian method for building the synthetic control in high-dimensional settings. Bayesian SCM has seen growing use in the recent literature for its flexibility in integrating with high-dimensional techniques. In variable selection context, Bayesian approach enables statistical inference by model averaging kass1995bayes, which contrasts with penalized regression methods like Lasso that focus on identifying a single best model. Further, Bayesian SCM provides posterior distributions for both ATT and the counterfactual outcome, which are more informative than point estimates. A potential limitation, however, is that the posterior fitting often requires the use of Markov chain Monte Carlo (MCMC) sampling, which can be very time-consuming compared to optimization-based methods. One recent work by goh2022synthetic considered Bayesian linear regression where the prior distribution of the regression coefficients is uniform over a convex hull, corresponding to the simplex constraint proposed by abadie2003economic. They also discussed the Bayesian formulation of the method of li2020statistical, which involves using a prior over the shifted conical hull. A similar model was studied in martinez2022bayesian, and a Bernstein-von Mises result for the Bayes estimator was proved. kim2020bayesian proposed using either the horseshoe prior or the spike-and-slab prior, two well-known Bayesian approaches to variable selection, to construct synthetic control with selected control units in high-dimensional settings; the simplex constraint was dropped in their method. It was demonstrated in their paper that their method outperforms those frequentist SCM methods based on variable selection. Nevertheless, to our knowledge, there is no existing Bayesian variable selection method for constructing synthetic control that also incorporates the simplex constraint.

Overview of this Work

Motivated by the high-dimensional challenge and the mixed benefits and limitations of the simplex constraint, we propose a new SCM based on a Bayesian model called BVS-SS (Bayesian Variable Selection with Soft Simplex constraint), which captures the sparsity among control units and allows the data to determine whether the simplex constraint is appropriate on a case-by-case basis. BVS-SS aims to leverage the strengths of both simplex-constrained SCM and unconstrained regression-based methods, which serves as an effective solution to the debate surrounding the use of this constraint. There is no similar method in the frequentist SCM literature to our knowledge.

In Section (ref) we formally introduce the BVS-SS model, which utilizes a standard spike-and-slab-type prior for selecting which control units should be used for constructing the synthetic control, a technique that has gained significant success and widespread popularity in the Bayesian community george1993variable, george1997approaches, brown1998multivariate. Different from existing models, BVS-SS assumes that there is an underlying sparse “mean weight vector” $\mu$ that satisfies the simplex constraint, but the actual weights $w_1, \dots, w_N$ involved in (ref) may deviate from $\mu$. Next, we introduce a key variance parameter $\tau$ quantifying the amount of this deviation. By imposing a non-informative prior on $\tau$, we let BVS-SS learn its value from the data, revealing whether the data actually agrees with the simplex assumption. A similar model was utilized in the theoretical analysis of martinez2022bayesian, but they did not consider variable selection and the implementation of the model.

To generate samples from the posterior distribution of BVS-SS, in Section (ref) we propose an original Metropolis-within-Gibbs sampler, where $\mu$ and error variance are updated by Gibbs schemes and $\tau$ is updated by Metropolis-Hastings steps. Our Gibbs-type updating for $\mu$ is different from existing Gibbs schemes for similar problems, as those methods typically only update one coordinate at a time, which is not feasible in our context due to the simplex constraint (i.e., given $\mu_{-j}$, we have $\mu_j = 1 - \sum_{i \neq j} \mu_i$ almost surely). We overcome this challenge by theoretically finding the joint conditional posterior distribution of two coordinates $(\mu_i, \mu_j)$ given all the other parameters, which enables us to update two coordinates simultaneously. This full conditional posterior distribution is complicated due to the spike-and-slab prior for variable selection, but sampling from it can be performed straightforwardly and efficiently. Using these MCMC samples, we can compute the posterior mean estimates and credible intervals of the counterfactual and ATT without extra computational cost.

In Section (ref), we theoretically investigate the effect of the simplex constraint and establish high-dimensional consistency results for BVS-SS with $\tau$ fixed to zero (i.e., regression coefficients are forced to satisfy the simplex constraint). In Theorem (ref), we show that the posterior inclusion probability of the true model converges to $1$ in probability under mild assumptions on the true data-generating process (DGP). While similar results have been obtained for various Bayesian variable selection models yang2016computational, the presence of the simplex constraint introduces a substantial new challenge for the analysis: the marginal likelihood for each model in BVS-SS does not admit a closed-form expression and instead involves an intractable integral over the simplex. To overcome this difficulty, we leverage techniques from constrained least-squares estimation and Gaussian concentration inequalities to derive closed-form upper and lower bounds for the marginal likelihood. The signal size required for detection has roughly the same order as that in the unconstrained case yang2016computational; see Assumption (ref). To our knowledge, this is the first consistency result for high-dimensional Bayesian variable selection under geometric constraints on the regression coefficients. Moreover, building on this result, we establish the order of the posterior expected predictive loss in Theorem (ref), providing theoretical guarantees on the ATT estimation. Our results complement the existing theoretical developments for frequentist SCM abadie2010synthetic, firpo2018synthetic, carvalho2018arco, li2020statistical, chernozhukov2021exact.

Section (ref) presents simulation studies which show that BVS-SS delivers robust and better performance than conventional methods, regardless of whether the simplex constraint is satisfied by the true DGP. Moreover, the posterior estimate of $\tau$ successfully captures the extent to which the true signals deviate from the simplex constraint. Section (ref) presents two real-data applications of ATT estimation with BVS-SS, one from carvalho2018arco on the effectiveness of an anti-tax evasion program and the other from shi2023forward on the impact of an anti-corruption policy. Section (ref) concludes the paper with a brief discussion. Details of the sampling algorithm, all proofs and additional numerical results are deferred to the supplementary materials.

The BVS-SS Model

Model and Prior

For generality and ease of notation, henceforth we use $Y \in \mathbb{R}^M$ to denote a response vector of $M$ observations and $X \in \mathbb{R}^{M \times N}$ to denote an arbitrary design matrix with $N$ explanatory variables. Assume that

where $w \in \mathbb{R}^N$ is the vector of unknown regression coefficients, $\epsilon$ is the vector of i.i.d. Gaussian noise with variance $\phi^{-1}$, and $I$ denotes the identity matrix. For the synthetic control problem described in Section (ref), we have $M = T_0$, $Y_t = Z_{1, t}$ and $X_{t j} = Z_{j + 1, t}$ for $t = 1, \dots, T_0$ and $j = 1, \dots, N$.

We propose a hierarchical prior distribution on $w$ which induces sparsity and incorporates the simplex constraint. To describe it, introduce an indicator vector $\gamma \in \{0, 1\}^N$ such that $\gamma_i = 1$ if and only if $w_i \neq 0$. Let $|\gamma| = \sum_{j = 1}^N \gamma_j$ denote the number of selected predictors. Given a vector $w$, we use $w_\gamma$ to denote the subvector of $w$ with entries indexed by $\{i \colon \gamma_i = 1\}$. We consider the following prior on $(w, \phi)$:

and when $\gamma_i = 0$, we let $\mu_i = w_i = 0$ with probability one. In (ref), $\kappa_1, \kappa_2, a_1, a_2, \alpha > 0$ and $\theta \in (0, 1)$ are fixed hyperparameters, $\text{sym-Dirichlet}(\alpha)$ denotes the symmetric Dirichlet distribution on the simplex with concentration parameter $\alpha$ (the dimension of the simplex is clear from context), and the Gamma distribution is in shape-rate parameterization. When $\alpha = 1$ and $|\gamma| = \ell$, the prior of $\mu_\gamma$ is a uniform distribution on $ \Delta^{|\gamma| - 1}$ with density $p(\mu_\gamma) = \Gamma(\ell) $ for each $\mu_\gamma \in \Delta^{|\gamma| - 1}$, where $\Gamma(\ell) = (\ell - 1)!$ is the inverse volume of the simplex.

In plain words, we assume that given $\gamma$ and $\phi$, $w_\gamma$ follows a normal distribution with mean vector $\mu_\gamma \in \Delta^{|\gamma| - 1}$ and prior variance depending on the parameter $\tau$. This conditional prior for $w$ can be viewed as an interpolation between the simplex constraint and the unconstrained setting. If the data suggests that $w$ is likely to satisfy the simplex constraint, then the posterior distribution of $\tau$ should concentrate around zero, while if the simplex constraint is significantly violated, $\tau$ should stay away from zero in the posterior. As will be demonstrated in our simulation study, $\tau$ can be reasonably estimated even with a relatively small sample size, which tells the user if the simplex constraint is appropriate for the data set being analyzed. This approach effectively reconciles the different views on the use of the simplex constraint, and thus we call our model BVS-SS (Bayesian Variable Selection with Soft Simplex constraint). We note that our prior for $w$ is similar to that studied in martinez2022bayesian; however, they did not consider variable selection, and in their implementation they assumed that $w$ satisfies the hard simplex constraint (though the soft constraint was used in their theory). Regarding the hyperparameters, $a_1, a_2$ should be chosen to reflect one's prior belief on whether simplex constraint is likely to hold, and a larger value of $\alpha$ encourages a more even distribution of the weights $\{ \mu_i \colon\gamma_i = 1\}$. In practice, one often has limited information about the true weight vector $w$, but the effect of these hyperparameters quickly becomes negligible as sample size increases.

The other elements of our model are standard. The prior for $\phi$ is conjugate, and one can choose small values for $\kappa_1, \kappa_2$ as a noninformative prior. The hyperparameter $\theta$ is the marginal prior probability of an explanatory variable being included in the regression model, and a smaller value of $\theta$ represents a heavier penalty on the model complexity.

Posterior Distribution and ATT Estimation

Let $X_\gamma$ denote the submatrix of $X$ with columns indexed by $\{i \colon \gamma_i = 1\}$. Under the prior distribution of BVS-SS, we have

A routine calculation using Woodbury identity yields the marginal likelihood of $(\mu, \tau, \phi)$:

Further, using the conjugacy of normal-gamma prior, we find that the conditional posterior distributions of $w_\gamma$ and $\phi$ are given by

Using the prior independence between $\mu, \tau, \phi$, we can express the joint posterior distribution by $p(\gamma, \mu_\gamma, \tau, \phi \mid y) \propto p(y \mid \gamma, \mu_\gamma, \tau, \phi)p(\mu_\gamma \mid \gamma) p(\gamma) p(\tau) p(\phi).$

Prediction with BVS-SS can be conducted as follows. Let $\tilde{X} \in \mathbb{R}^{\tilde{M} \times N}$ denote the design matrix of another $\tilde{M}$ observations, and suppose we are interested in estimating $\tilde{X} w$. Letting $\mathbb{E}_y$ denote the posterior distribution given $Y = y$ and using (ref), we obtain that

Given $n$ MCMC samples $(\gamma^{(k)}, \mu_\gamma^{(k)}, \tau^{(k)}, \phi^{(k)})_{k = 1}^n$, we can approximate the above expectation using the sample average, which yields the posterior predictive mean of $\tilde{X} w$. The covariance matrix of the posterior predictive distribution of $\tilde{X} w$ can be calculated similarly.

For synthetic control, the goal is to estimate ATT among the $\tilde{M}$ observations. Denote their response under the treatment by $\tilde{Y}^{(1)}$, which is assumed observed, and denote the counterfactual response without treatment by $\tilde{Y}^{(0)}$. Let $\delta = \tilde{Y}^{(1)} - \tilde{Y}^{(0)}$, and define ATT by

which is the main parameter of interest. By replacing $\tilde{Y}^{(0)}$ using the synthetic control $\mathbb{E}_y[ \tilde{X} w ]$, we obtain an estimator for $\bar{\delta}$. If we have samples $(w^{(k)})_{k=1}^n$ approximating the marginal posterior distribution of $w$, we can also obtain the posterior predictive distribution of ATT by plugging $w^{(k)}$ into the estimator $ \tilde{M}^{-1}\sum_{i = 1}^{\tilde{M}} ( \tilde{Y}^{(1)}_i - w^\top \tilde{X}_{(i)} )$, where $\tilde{X}_{(i)}$ denotes the $i$-th row of $\tilde{X}$ (treated as a column vector).

Metropolis-within-Gibbs Sampling for BVS-SS

The standard solution to posterior computation for the spike-and-slab variable selection is to run a Metropolis--Hasting sampler targeting the posterior distribution of $\gamma$ by add/delete/swap proposals; see george1997approaches, guan2011bayesian, chang2024dimension, among many others. However, this approach is not applicable in our setting since the likelihood $p(y \mid \gamma, \phi)$ in BVS-SS does not have a closed-form expression.

Instead of directly sampling from the marginal posterior of $\gamma$, we devise a Metropolis-within-Gibbs scheme targeting the full posterior distribution $p(\gamma, \mu, \tau, \phi \mid y)$. The Metropolis--Hastings update for $\tau$ and the Gibbs update for $\phi$ are standard, which we detail in Appendix (ref). Since $\gamma_i = \mathbbm{1}_{ \{\mu_i \neq 0\}}$, we update $(\gamma_i, \mu_i)$ simultaneously, which can be implemented straightforwardly as a Gibbs step in the unconstrained spike-and-slab variable selection. However, this is again infeasible for BVS-SS due to the simplex constraint: given $\mu_{-j}$, we have $\mu_j = 1 - \sum_{i \neq j} \mu_i$ almost surely. Therefore, we propose to a novel Gibbs scheme where we update two coordinates of $(\gamma, \mu)$ simultaneously.

Let $\mu_{-(i, j)}$ denote the subvector of $\mu$ with $(\mu_i, \mu_j)$ removed and define $\gamma_{-(i, j)}$ analogously. To implement the Gibbs update of $(\gamma_i, \gamma_j, \mu_i, \mu_j)$ for $i \neq j$, we need to find the full conditional posterior $p(\gamma_i, \gamma_j, \mu_i, \mu_j \mid y, \gamma_{-(i, j)}, \mu_{-(i, j)}, \tau, \phi)$. If $\sum_{k \neq i, j} \mu_k = 1$, neither $i$ nor $j$ can be selected due to the simplex constraint, and thus $\gamma_i = \gamma_j = \mu_i = \mu_j = 0$. Now assume $\sum_{k \neq i, j} \mu_k < 1$, in which case we must have $(\gamma_i, \gamma_j) = (1, 0), (0, 1)$ or $(1, 1)$. The conditional posterior probability $p( \gamma_i = 1, \gamma_j = 0 \mid y, \gamma_{-(i, j)}, \mu_{-(i, j)}, \tau, \phi)$ is easy to evaluate (up to a normalizing constant) since the simplex constraint implies $\mu_i = 1 - \sum_{k \neq i, j} \mu_k$. The case $\gamma_i = 0, \gamma_j = 1$ can be analyzed similarly. When $\gamma_i = \gamma_j = 1$, we first find the full conditional posterior $p( \mu_i, \mu_j \mid y, \gamma_{-(i, j)}, \mu_{-(i, j)}, \tau, \phi, \gamma_i = \gamma_j = 1)$, which is a degenerate truncated normal distribution due to the constraint $\mu_i > 0, \mu_j > 0$ and $\mu_i + \mu_j = 1 - \sum_{k \neq i, j}$. This is derived in Lemma (ref). Integrating out $(\mu_i, \mu_j)$, we get the conditional posterior probability $p( \gamma_i = \gamma_j = 1 \mid y, \gamma_{-(i, j)}, \mu_{-(i, j)}, \tau, \phi)$ up to a normalizing constant. Normalizing over the three cases yields the posterior distribution $p( \gamma_i, \gamma_j \mid y, \gamma_{-(i, j)}, \mu_{-(i, j)}, \tau, \phi)$, which has a closed-form expression in terms of the cumulative distribution function of the normal distribution; see Lemma (ref). The resulting sampling algorithm is summarized in Algorithm (ref), where $\beta_{i,j}, \Lambda_{i,j}$ are defined in Lemma (ref). More details are provided in Appendix (ref).

Theoretical Results for BVS-SS

In this section, we study the theoretical properties of BVS-SS. For ease of analysis, we fix the hyperparameter $\tau$ and consider two complementary scenarios: $\tau = 0$ and $\tau \rightarrow \infty$. Following the convention in the variable selection literature yang2016computational, we interpret a model $\gamma$ as both a binary vector and a subset of $[N]:=\{1,2,\dots,N \}$; that is, to indicate that the $i$-th covariate is selected, we can write either $\gamma_i = 1$ or $i \in \gamma$. We use $|\cdot|$ to denote set cardinality, and thus $|\gamma|$ denotes the model size.

High-dimensional Selection Consistency with \texorpdfstring{$\tau = 0$}{tau = 0}

For the first scenario, we fix the prior parameter $\tau = 0 $ and assume that the true DGP satisfies the simplex constraint. Our goal is to establish, in high-dimensional regimes, the “strong selection consistency” property of BVS-SS which means that the posterior probability of the true model, denoted by $\gamma^*$, converges to one in probability under the true DGP; here, “strong” means that this property is stronger than other consistency criteria such as pairwise selection consistency casella2009consistency. Since $\tau = 0$ implies $w = \mu$ almost surely, we can rewrite the model as

For technical convenience, we slightly modify the prior specification as follows:

where $\mathbb{S}_L = \left\{ \gamma \subset [N] \colon 1 \leq |\gamma| \leq L \right\}$ denotes the sparse model space with size bounded by $L$, and Dir$(\mathsf{1})$ denotes the Dirichlet distribution with parameter vector $\mathsf{1} = (1, \dots, 1)$ (i.e., uniform distribution over the simplex $\Delta^{|\gamma| - 1}$). Compared to the original prior specification given by (ref), the changes are minor. First, rather than letting $\gamma_i \overset{\mathrm{i.i.d.}}{\sim} \mathrm{Bernoulli}(\theta)$, we reweight the prior probability of each $\gamma$ by three factors, where $\Gamma(|\gamma|)^{-1}$ is the volume of the simplex $\Delta^{|\gamma| - 1}$, and $\mathsf{1}^\top (X_\gamma^\top X_\gamma)^{-1} \mathsf{1}$ is the variance of $\mathsf{1}^\top \hat{\mu}_\gamma$ with $\hat{\mu}_\gamma$ denoting the OLS estimator for model $\gamma$. This adjustment serves only to simplify the presentation, and we prove in Lemma (ref) below that the effect of this reweighting is typically negligible when $\theta$ decays much faster than $\sqrt{M}$. The sparsity condition $|\gamma| \leq L$ is standard and necessary for high-dimensional consistency analysis. Second, we let the prior of $\phi$ depend on $\gamma$ so that for a larger model, the prior mean of $\phi$ is larger (i.e., the error variance is smaller). This modification enables a simpler comparison of the marginal likelihoods of different models after integrating out $\phi$. When $\kappa_1$ is large (e.g. of order $M$), the difference between the two prior specifications for $\phi$ becomes negligible. Third, for the conditional prior for $\mu_\gamma$ given $\gamma$, we set $\alpha = 1$ in (ref) to ensure that the prior density $p(\mu_{\gamma} \mid \gamma)$ only depends on $|\gamma|$.

By integrating out $\mu_\gamma$, we can express the posterior density of $(\gamma, \phi)$ as

Unlike the unconstrained Bayesian linear regression, $p(\gamma, \phi \mid y) $ is not available in closed form due to the intractable integral involved in $p(y \mid \gamma, \phi)$. Consequently, the marginal posterior probability of $\gamma$, $p(\gamma \mid y) \propto p(\gamma) \int p(\phi \mid \gamma) p(y \mid \gamma, \phi) \mathrm{d} \phi$, has no closed-form expression either. To establish the high-dimensional consistency result, we make the following assumptions on the design matrix, prior parameters and true DGP, where all parameters involved may vary with the sample size $M$, except the universal constants denoted by $c_{\mu}, c_{\theta}, c_M$ and $c_j$ for integer $j$. In particular, the number of variables $N$, the maximum model size $L$, and the true model size $\ell^* = |\gamma^*|$ are all allowed to go to infinity with $M$.

Assumption (ref) requires that each column of $X$ has the same magnitude and that the minimum eigenvalue of the Gram matrix is well controlled for any model of size at most $L$. It is a mild condition, often known as “restricted eigenvalue,” commonly used in the literature shang2011consistency,narisetty2014bayesian,yang2016computational. Assumption (ref) provides guidance for choosing the prior: the prior inclusion probability $\theta$ should be small enough in order to penalize large models and avoid overfitting. Similarly, Assumption (ref) requires that the choices of $\kappa_1, \kappa_2$ are not too extreme. Assumption (ref) assumes that the true DGP satisfies the simplex constraint, and (ref) imposes a lower bound on the smallest nonzero coefficient of the true DGP. This is also a standard assumption and often known as the “$\beta$-min condition.” Intuitively, this assumption is needed since coefficients close to zero cannot be consistently identified through variable selection. Note that the order of the detection threshold required by (ref) is essentially the same as that in the unconstrained Bayesian variable selection yang2016computational. Finally, Assumption (ref) specifies the sample size needed for consistent variable selection. We use $\mathbb{P}^*$ to denote the probability measure associated with the true DGP specified in Assumption (ref).

Under the above assumptions, we prove that $p(\gamma^* \mid y)$ converges to one in probability with respect to $\mathbb{P}^*$ given sufficiently large $c_{\mu}, c_{\theta}$ and $c_M$. A key step in our proof is to establish the following posterior probability ratio bound.

By Theorem (ref), we can bound $p(\gamma \mid y)/p(\gamma^* \mid y)$ by analyzing the behavior of unconstrained least-squares estimators, which is well understood in the variable selection literature. This analysis then yields the desired strong selection consistency.

The strong selection consistency ensures that the true model $\gamma^*$ can be identified with high probability. Next, we consider the estimation of the regression coefficient vector $\mu$ and the task of prediction. The posterior expected $L^2$ loss for estimating $\mu$ can be expressed as

where we recall $\mathbb{E}_y$ denotes the expectation over the posterior distribution given $Y = y$. Similarly, given a new design matrix $\tilde{X}$, we can define the posterior expected predictive loss by $\mathbb{E}_y [ \| \tilde{X} \mu - \tilde{X} \mu^* \|_2^2 ]$. In Theorem (ref), we characterize the rate at which the two loss functions go to zero.

The consistency of our ATT estimator follows from the the prediction loss bound in Theorem (ref). We expect that this result can be extended to more general scenarios, including cases where $[X^\top \; \tilde{X}^\top]$ is a stationary process li2020statistical.

Limiting Behavior when \texorpdfstring{$\tau \rightarrow \infty$}{tau rightarrow infty}

We also investigate the performance of BVS-SS when $\tau \to \infty$. In this case, the prior places increasing mass on values of $w$ that do not satisfy the simplex constraint. We compare the behavior of the posterior distribution $p(\gamma \mid y)$ with the conventional unconstrained version, which we denote by $\tilde{p}(\gamma \mid y)$ (see Appendix (ref) for details). The following theorem shows that the two posterior distributions become essentially the same as $\tau \rightarrow \infty$.

In Theorem (ref), we treat the data $(X, y)$ as fixed and let $\tau \rightarrow \infty$, but we expect that our argument can be extended to a fully non-asymptotic analysis. Such an extension may be combined with existing variable selection consistency results to show that our model also achieves selection consistency when the true DGP significantly violates the simplex constraint but $\tau$ is sufficiently large.

Simulation Studies

Simulation Settings

We simulate the observed training data $X \in \mathbb{R}^{M \times N}$ and $Y \in \mathbb{R}^{M}$ by $Y = X w^* + \epsilon,$ where $\epsilon_i \overset{\mathrm{i.i.d.}}{\sim} N(0, 1/ \phi^* )$, and $w^* \in \mathbb{R}^N$ is the true regression coefficient vector such that $w^*_j = 0$ for $j = J+1, \dots, N$ (that is, only the first $J$ predictors have nonzero effects on $Y$). Similarly, we simulate the data under treatment $\tilde{X} \in \mathbb{R}^{\tilde{M} \times N}$ and $\tilde{Y}^{(1)} \in \mathbb{R}^{\tilde{M}}$ by

The vector $\delta \in \mathbb{R}^{\tilde{M}}$ represents the treatment effects, and $\delta^* \in \mathbb{R}$ denotes the true ATT. In our simulation, we use $M = \tilde{M} \in \{25, 50, 100, 200\}$ and $N \in \{20, 50\}$. When $N = 20$, we set $J = 5$, and when $N = 50$, we set $J = 10$. Given $J$, we define a vector $\mu^* \in \Delta^{N - 1}$ by

Then, we define $w^*$ by $w^* = \lambda \mu^*$ where $\lambda > 0$ is a scaling factor. When $\lambda = 1$, the true DGP satisfies the simplex constraint. Assuming $\lambda > 0$, we have $\lambda = \|w^*\|_1$. We consider $\lambda \in \{1, 2, 3\}$ in our simulation studies. Given $w^*$, we set the variance parameters by $v^* = \phi^* = 4 / \| w^* \|_2^2$. For each choice of $(M, N, \lambda)$, we generate 100 replicates where entries of $X, \tilde{X}$ are also sampled independently from the standard normal distribution. See Appendix (ref) for simulation results where the true DGP is driven by a factor model. We fix true ATT $\delta^* = 0.5$ throughout our simulation studies.

We set the hyperparameters of BVS-SS by $\kappa_1 = \kappa_2 = 1, a_1 = 0.01, a_2 = 0.1, \alpha = 1, \theta = 0.2.$ The choice of $(a_1, a_2)$ guarantees that $\tau$ has prior mean $0.1$ and prior variance $1$. Further, for numerical stability, we truncate the prior distribution of $\tau$ by assuming that $\tau \geq 10^{-6}$. For every simulated data set, we initialize Algorithm (ref) at $\tau^{(0)} = \phi^{(0)} = 1$ and $\mu^{(0)}$ sampled from its prior distribution, and we run Algorithm (ref) for $1,000$ iterations and discard the first $500$ iterations as burn-in.

For comparison, we have also implemented the following methods for ATT estimation. First, we use the R package glmnet to perform variable selection and ATT estimation via Lasso; the parameter of Lasso is chosen by cross-validation. Second, we compute the OLS estimator for $w^*$ using the entire data set $(X, Y)$ without variable selection, denoted by $\hat{w}_{\mathrm{LS}}$. Then ATT is estimated by

where $\tilde{X}_{(i)}$ is the $i$-th row of $\tilde{X}$ treated as a column vector. We also compute the simplex-constrained estimator $\hat{w}_{\mathrm{QP}}$ (QP stands for “quadratic programming”) defined by $\hat{w}_{\mathrm{QP}} = \operatorname*{argmin}_{w \in \Delta^{N-1}} \| Y - X w \|_2^2,$ which was introduced in abadie2003economic and can be obtained by quadratic programming goldfarb2006dual, abadie2011synth. The corresponding ATT estimator $\widehat{\mathrm{ATT}}_{\mathrm{QP}}$ is obtained similarly by replacing $\hat{w}_{\mathrm{LS}}$ with $\hat{w}_{\mathrm{QP}}$ in (ref). Finally, we consider two oracle estimators that are not implementable in practice. Denote by $\gamma^* = \gamma(w^*)$ the true indicator vector corresponding to $w^*$. We compute the OLS and QP estimators using $(X_{\gamma^*}, Y)$. That is, these two estimators have access to $\gamma^*$ (which is not possible in reality) and will be used for comparison in our simulation studies. Note that since the regression coefficient estimator of Lasso and OLS are scale-equivariant, the ATT estimators of Lasso and two OLS methods does not depend on the signal size factor $\lambda$.

Estimation of Average Treatment Effect

For each method considered, we compute the mean squared error of ATT estimation over 100 replicates; denote it by $\mathrm{MSE}(\widehat{\mathrm{ATT}})$. To measure the efficiency of an ATT estimator, we consider the oracle OLS estimator OLS$(\gamma^*)$ as the reference estimator, and define the relative efficiency of an estimator $\widehat{\mathrm{ATT}}$ by

If $ \mathrm{RE}( \widehat{\mathrm{ATT}} ) = 100\%$, then this estimator is as efficient as the oracle OLS estimator.

We show the relative efficiency of the other five estimators with $\lambda = \|w^* \|_1 = 1, 2, 3$ in Figure (ref). When $\lambda = 1$, our method outperforms Lasso substantially when $M = 25$ or $50$, which shows that some prior knowledge about the simplex constraint can be very helpful when the sample size is small or moderate. By comparing BVS-SS with the two QP methods, we see that performing variable selection also significantly improves ATT estimation when the sample size is small. Indeed, when $N = 50, M = 25$, QP with the entire data set is not feasible since the underlying optimization problem does not have a unique solution. We also note that as sample size increases, the performance of BVS-SS quickly approaches that of the two oracle estimators, OLS$(\gamma^*)$ and QP$(\gamma^*)$. OLS with the entire data set performs poorly unless the sample size is very large, which is also expected since it does not perform variable selection.

When $\lambda = 2$ or $3$, the simplex constraint is violated by the true DGP. Though BVS-SS utilizes the simplex constraint, Figure (ref) shows that its performance remains robust. BVS-SS outperforms Lasso in all scenarios, and when the ratio $N/M$ is large, the advantage of BVS-SS is very significant. Indeed, we observe that when $N = 50, M = 25$, Lasso almost fails to recover any signal contained in $w^*$ (see also Table (ref)). When the sample size becomes sufficiently large (say, $M \geq 100$), BVS-SS again achieves close-to-optimal performance, and it outperforms QP$(\gamma^*)$ by a wide margin, which knows the true model but assumes the simplex constraint. The robust and superior performance of BVS-SS across the three settings is largely due to the use of the soft simplex constraint, which enables BVS-SS to adapt to the unknown level of $\|w^*\|_1$ in the given data.

Variable Selection and Estimation of Variance Parameters

We now examine the performance of BVS-SS in terms of variable selection and the estimation of the two variance parameters $\tau, \phi$. First, we compare the variable selection result of BVS-SS with that of Lasso in Table (ref), where “$\ell^1$-loss” denotes the $\ell^1$-distance between the selected model $\hat{\gamma}$ and the true model $\gamma^*$, and “model size” refers to the cardinality of $\hat{\gamma}$ (i.e., the number of nonzero entries of $\hat{w}$). For BVS-SS, the two statistics are averaged over the collected MCMC samples. It is evident that BVS-SS has a much higher accuracy than Lasso across all settings. Especially when $M$ is large, BVS-SS can correctly identify most entries of $\gamma$, while the performance of Lasso remains similar to that observed with small $M$.

Next, we visualize the distribution of the posterior mean estimate of $\phi$ (error precision) across 100 replicates in Figure (ref). This statistic reflects how well the BVS-SS model fits the data. We only show the result for $N = 50$, as the distribution for $N = 20$ is similar. Observe that regardless of $\|w^*\|_1$, the posterior mean of $\phi$ increases with $M$, which is expected since a larger sample size enables BVS-SS to detect more signals in the data, resulting in a smaller estimate of the error variance. Notably, when $\|w^*\|_1 = 2$ or $3$, the posterior estimation of $\phi$ appears highly accurate for $M \geq 100$, indicating that our model achieves an optimal balance between model complexity and fitting.

Finally, we examine the posterior mean estimate of $\tau$ when $N = 50$, for which the result is presented in Figure (ref). As discussed in Section (ref), $\tau$ can be considered as an indicator of whether the simplex constraint is satisfied by the given data, and we now explain why Figure (ref) provides compelling empirical evidence supporting this. When $\|w^*\|_1 = 1$, we see that the posterior mean of $\tau$ decreases as $M$ increases, since a larger sample size offers stronger evidence that the simplex constraint is satisfied, which draws $\tau$ towards zero. In contrast, when $\|w^*\|_1 = 2$ or $3$, the trend reverses, since a larger sample size enables BVS-SS to detect deviations from the simplex constraint. When $M = 200$, we find that the ratio of the posterior mean estimate of $\tau$ and that of $\phi$ equals $0.0003$ for $\|w^*\|_1 = 1$, $0.024$ for $\|w^*\|_1 = 2$, and $0.081$ for $\|w^*\|_1 = 3$. This actually aligns well with the true data-generating mechanism: since we use $J = 10$ for $N = 50$, we may estimate the “true” mean squared deviation of $w^*$ from the simplex constraint by $0.1 \times \sum_{j=1}^{10} (w^*_j - 0.1)^2$, which equals $0.021$ for $\|w^*\|_1 = 2$ and $0.065$ for $\|w^*\|_1 = 3$.

Empirical Examples

In this section, we revisit two empirical examples from the SCM literature. We will use the notation introduced in Section (ref): $T_0$ is the number of pre-treatment time periods, $T_1$ is the number of post-treatment time periods, the total number of units is $N + 1$ where the first one is used as the response and the other $N$ units as the explanatory variables. (For BVS-SS, $M = T_0$ and $\tilde{M} = T_1$.) The first example has $N < T_0$, while the second has $N \gg T_0$, and it will be shown that our method works well in both settings.

Nota Fiscal Paulista: Anti-tax Evasion

An anti-tax evasion program, Nota Fiscal Paulista (NFP), was implemented in S\ {a}o Paulo, Brazil in October 2007 with the goal of reducing tax evasion by incentivizing consumers to request electronic receipts in exchange for potential tax rebates and participation in monthly lotteries promoted by the government. However, those restaurant owners and retailers may pass on part of the tax costs to consumers by raising product prices. Although the NFP was implemented sequentially across various sectors, including restaurants, bakeries, bars, and food service retailers, the literature has suggested that the food away from home (FAH) index could be used as a suitable indicator for price levels of these sectors. carvalho2018arco proposed an artificial counterfactual (ArCo) estimation procedure, by connecting Lasso to SCM without simplex constraint, to explore how the NFP affects the inflation on FAH. They also included monthly GDP growth, retail sales growth and monthly credit growth. Their ArCo estimation indicated that the ATT of NFP is 0.4478%.

We obtain the data analyzed in carvalho2018arco from the $\texttt{R}$ package $\texttt{ArCo}$. The pre-treatment periods span $T_0 = 33$ months, and the post-treatment periods span $T_1 = 23$ months. Slightly different from the original study, the data includes only $N = 8$ control units, which is one fewer than that mentioned in their paper, and there is no data related to credit and retail sales. We run our Algorithm (ref) for 1000 iterations to estimate the impact of NFP, with the first half of the samples dropped for burn-in. We set the hyperparameters by $\kappa_1 = \kappa_2 = 1, a_1 = 0.01, a_2 = 0.1, \alpha = 1, \theta = 0.25.$ Table (ref) summarizes the posterior mean and 95% credible intervals of ATT, $\tau$, $\phi$ and the model size $|\gamma|$. The ATT estimation is close to that of carvalho2018arco, but our credible interval is much smaller compared to their confidence set. While the analysis of carvalho2018arco selected all the predictors to construct the synthetic control, our method only chooses about two predictors on average. This is consistent with the findings of our simulation study presented in Section (ref), which shows that Lasso has the tendency of over-estimating the model size compared to BVS-SS. In Appendix (ref), we provide the trace plots for $\phi$ and $\log(\tau)$ and visualize the posterior distributions of ATT, $\log(\tau)$, $\phi$, model size $|\gamma|$, and the counterfactual estimate are provided.

China's Anti-corruption Campaign

shi2023forward investigated how China's unprecedented anti-corruption campaign, launched in November 2012, affected the importation of luxury watches. This was motivated by the fact that high-end imported watches are often used as discreet bribery gifts in place of blatant cash payments, and empirical studies have shown that the imports of luxury watches in China move with government leadership changes lan2018swiss. shi2023forward used data from the United Nations Comtrade Database, taking the monthly growth rate of the category “watches with case of, or clad with, precious metal” as the outcome and constructing the synthetic counterfactual from $N = 87$ other commodity categories. January 2013 is considered the time of treatment, which was the month immediately following the announcement of the Eight-Point anti-corruption policy in China. The pre-treatment period spans from February 2010 to December 2012, yielding $T_0 = 35$ observations, and the post-treatment period covers January 2013 to December 2015, yielding $T_1 = 36$ observations. Their forward selection algorithm identified three categories as control units: “knitted or crocheted fabric”, “cork and articles of cork”, and “salt, sulfur, earth, stone, plaster, lime, and cement”. Their estimated counterfactual predicted that without the anti-corruption campaign, luxury watch imports would have continued to increase by 1.7% in January 2013, while in reality imports dropped by 42%. Their treatment effect estimation indicated a reduction of 3.09% in luxury watch imports per month over the post-treatment period.

We obtain the data of shi2023forward from the $\texttt{R}$ package $\texttt{fdPDA}$, which contains all the 88 commodity categories (including the response and $N=87$ control categories), and run our Algorithm (ref) for 1000 iterations with the same hyperparameters as in the last example. Table (ref) indicates that our ATT estimation (-2.1%) is smaller than that of shi2023forward in absolute value, but the 95% credible interval indicates that the anti-corruption has a significant effect on the high-end watch imports. On average, our model size $|\gamma|$ is larger than that of shi2023forward (who only selected 3 control categories), but Figure (ref) shows that the posterior distribution of $|\gamma|$ still concentrates on small models with size $\leq 5$. Additional simulation results are provided in Appendix (ref).

It is worth noting that in our posterior samples, the categories “cork” and “salt” appear with negligible probabilities. Instead, the most frequently selected categories are those related to textiles---for example, “knitted or crocheted fabrics” (51.4%), “special woven fabrics such as lace, tapestries, trimmings, and embroidery” (16.6%), and “wool, fine or coarse animal hair” (12%). In addition, the category “optical, photographic, cinematographic, measuring, checking, medical or surgical instruments and apparatus; parts and accessories” also appears with a relatively high probability (10.2%). Compared with “cork” and “salt”, these categories seem more plausible to serve as control units. The textile products share decorative features similar to those of wristwatches, while the optical and photographic instruments also rely on precision mechanical engineering.

Conclusion and Future Directions

We propose a novel Bayesian synthetic control method, BVS-SS, that integrates a relaxed simplex constraint with Bayesian spike-and-slab variable selection. Our approach introduces a hierarchical prior to determine how closely the data should adhere to the simplex constraint, a feature that is missing in existing frequentist or Bayesian synthetic control frameworks and helps mitigate the ongoing debate surrounding the use of simplex constraints. To efficiently compute the posterior distribution of BVS-SS, we develop a novel Metropolis-within-Gibbs sampler, which overcomes the hindrance caused by the simplex constraint by updating two regression coefficients simultaneously from the full conditional posterior. We also show that in high-dimensional settings, when the true DGP satisfies the simplex constraint, BVS-SS with $\tau = 0$ consistently selects the true control units and achieves consistency in terms of both regression coefficient estimation and prediction. Simulation studies and real-data examples illustrate the advantages and usefulness of our proposed algorithm.

There are several interesting directions that future researchers may dive in. First, it is known that Gibbs sampling methods are susceptible to high collinearity among the explanatory variables, which is a common situation in SCM applications where factor models are often used to model the DGP. Though it is unclear if such issues affect the methodology proposed in this work (as our algorithm updates two coordinates at a time), it would be interesting to investigate the use of other MCMC techniques, such as pseudo-marginal sampling andrieu2009pseudo. Second, it will be interesting to further generalize BVS-SS to account for time-invariant characteristics, which are often used together with time-variant predictors in the SCM literature for constructing counterfactual outcomes. Third, more numerical studies or model extensions can be performed to study the performance of BVS-SS when the model is mis-specified (e.g., the data is non-stationary or contains heterogeneous errors), which would be of great interest to economists. We hope our method will serve as an inspiration for future research and facilitate the widespread use of Bayesian synthetic control methods.

Acknowledgements

Y.X. is grateful to Qi Li, Yonghong An and Angda Li for their helpful comments and suggestions. This work was supported in part by NSF via grants DMS-2311307 and DMS-2245591.

\setcounter{section}{0} \setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{lemma}{0} \setcounter{proposition}{0} \setcounter{corollary}{0} \setcounter{assumption}{0} \setcounter{figure}{0} \setcounter{table}{0}

\phantomsection

Details of the Metropolis-within-Gibbs Sampler

For Bayesian variable selection involving a prior on the indicator vector $\gamma$ as we have considered, the standard approach to posterior sampling is to integrate out the regression coefficient vector and use add-delete-swap proposals to construct a Metropolis--Hastings algorithm targeting the marginal posterior distribution of $\gamma$ (or the joint distribution of $\gamma$ and some variance parameter); see, e.g., george1997approaches, chipman2001practical, guan2011bayesian, zhou2022dimension. However, this approach is not applicable to our model, since $\mu_\gamma$ cannot be integrated out in closed form from the joint posterior distribution of $(\gamma, \mu_\gamma, \tau, \phi)$. To tackle this very unique challenge, we propose a novel Metropolis-within-Gibbs sampler targeting the joint posterior distribution of $(\mu, \tau, \phi)$. Note that since $\gamma$ can be viewed as a function of $\mu$ with $\gamma_i = \mathbbm{1}_{ \{\mu_i \neq 0\}}$, henceforth we will denote the conditional distribution in (ref) by $p(y \mid \mu, \tau, \phi)$, and $\gamma$ should be understood as $\gamma = \gamma(\mu)$.

Updating of \texorpdfstring{$\tau$ and $\phi$}{tau and phi}

The updates for $\tau$ and $\phi$ are standard. We update $\phi$ by drawing it from the full conditional distribution given in (ref), and we use a Metropolis--Hastings update for $\tau$ that is invariant with respect to $p(\tau \mid y, \mu, \phi)$. Explicitly, we propose a new value for $\tau^*$ by a Gaussian random walk on log-scale with variance $\eta > 0$ and accept $\tau^*$ with probability

where the likelihood is given in (ref) and $p(\tau)$ denotes the prior density. In Algorithm (ref), we repeat this step $n_\tau$ times so that $\tau$ is more likely to be updated to some value with higher conditional posterior density given $\mu$ and $\phi$.

Updating of \texorpdfstring{$\mu$}{mu}

For $\mu$, we propose a Gibbs-type updating scheme, which is quite different from any existing Gibbs scheme for similar problems. Let $\mu_{-j}$ denote the subvector of $\mu$ with $\mu_j$ removed. Since the prior of BVS-SS assigns probability one to $\mu \in \Delta^{N - 1}$, given $\mu_{-j}$, we have $\mu_j = 1 - \sum_{i \neq j} \mu_i$ almost surely, and thus the full conditional distribution $p(\mu_j \mid y, \mu_{-j}, \tau, \phi)$ is also degenerate. Hence, a Gibbs scheme that updates each $\mu_j$ from its full conditional posterior will not move at all. This motivates us to devise a more complicated updating scheme by updating $\mu_i, \mu_j$ (with $i\neq j$) simultaneously from the full conditional posterior $p(\mu_i, \mu_j \mid y, \mu_{-(i, j)}, \tau, \phi)$, where $\mu_{-(i, j)}$ denotes the subvector of $\mu$ with the $\mu_i, \mu_j$ removed.

Let $\mu_{-(i, j)}, \tau, \phi$ be given and define $s = 1 - \sum_{k \neq i, j} \mu_k$. If $s = 0$, the simplex constraint implies that $\mu_i = \mu_j = 0$ almost surely. If $s \in (0, 1]$, the conditional posterior distribution of $\mu_i, \mu_j$ falls into one of the following three cases:

Therefore, to find a closed-form expression for this conditional distribution, it suffices to find $p(\gamma_i, \gamma_j \mid y, \mu_{-(i, j)}, \tau, \phi)$ and the conditional distribution $p(\mu_i \mid y, \mu_{-(i, j)}, \tau, \phi, \gamma_i = \gamma_j = 1)$.

Let $\gamma_{-(i, j)} = \{ k \notin \{i, j\} \colon \mu_k \neq 0\} \subset [N]$ denote the model corresponding to $\mu_{- (i, j)}$. We can define the four possible models under consideration by

We will also interpret $\gamma^0, \gamma^i, \gamma^j, \gamma^{ij}$ as indicator vectors in $\{0, 1\}^N$. We first find the full conditional distribution of $(\mu_i, \mu_j)$ given other parameters and $\gamma_i = \gamma_j = 1$ in the following lemma. For simplicity, we assume $\alpha = 1$ in this section, and in Remark (ref) we explain how to extend the results to other integer values of $\alpha$.

Next, we find the conditional posterior probabilities of $\gamma^0, \gamma^i, \gamma^j$, and $\gamma^{ij}$.

Algorithm

The pseudocode for our proposed Metropolis-within-Gibbs sampler targeting $p(\mu, \tau, \phi \mid y)$ is given in Algorithm (ref) in the main text. The most computationally intensive part of Algorithm (ref) is the evaluation of expressions of the form $V_{\gamma, \tau}^{-1} z$ for a given vector $z$, which may seem highly time-consuming if $|\gamma|$ is large. But this can be implemented efficiently by computing and updating the Cholesky decomposition of $V_{\gamma, \tau}$, a technique commonly employed in Bayesian spike-and-slab variable selection smith1996nonparametric, george1997approaches. More specifically, at the beginning of the $t$-th iteration, we compute the Cholesky decomposition of $V_{\gamma, \tau^{(t-1)}}$, where $\gamma$ is the model selected in the last iteration. Then, each time we draw $(\mu_i, \mu_j)$ from its full conditional posterior, we need to change one or two coordinates of $\gamma$, and the resulting Cholesky decomposition of $V_{\gamma, \tau^{(t-1)}}$ can be obtained by standard updating algorithms golub2013matrix.

Certain variations of Algorithm (ref) are also straightforward to implement. For example, one can replace the fixed-scan Gibbs updating of $\mu$ in Algorithm (ref) by a random-scan update: in the $t$-th iteration, we uniformly sample a pair $i < j$ such that $\mu^{(t-1)}_i + \mu^{(t-1)}_j > 0$ and update $(\mu^{(t)}_i, \mu^{(t)}_j)$ from its conditional posterior distribution (an acceptance-rejection step is needed to correct for the proposal bias), and then draw $\phi^{(t)}$ from the full conditional posterior and update $\tau^{(t)}$ by one Metropolis--Hastings step. In this case, the iterative complex factorization algorithm proposed by zhou2019fast can be used to efficiently solve the linear system $V_{\gamma, \tau}^{-1} z$.

Constrained Least-squares Estimators

In this section, we review some known results about the least-squares estimator under a linear constraint amemiya1985advanced. We present the results in the general form although for our analysis, we only need to consider the constraint $\mathsf{1}^\top \mu = 1$, which is part of the simplex constraint.

Next, we prove some linear algebra results that will be useful for analyzing the constrained least-squares estimator.

Proofs for Section (ref)

As in Section (ref), we use $|\gamma|$ and $\ell$ interchangeably. All corresponding indices and symbols (e.g., star, tilde, etc.) are defined equivalently.

Proof of Lemma (ref)

Proof of Theorem (ref)

We divide the proof into multiple parts. First, in Section (ref), we construct the events that happen with probability at least $1 - c_1 N^{-c_2}$ under the true data-generating probability measure $\mathbb{P}^*$, where $c_1, c_2 > 0$ are some universal constants. In Section (ref), we derive an explicit expression for $p(\gamma, \phi \mid y)$ by using the property of the constrained least-squares estimators. In Section (ref), we find a lower bound on $p(\gamma^*, \phi \mid y)$ by analyzing both the unconstrained and constrained least-squares estimators under $\gamma^*$, In Section (ref), we prove Theorem (ref) by integrating over $\phi$. Throughout our proof, we use

to denote the projection matrix onto the column space of $X_\gamma$.

High-probability Events

We begin by constructing the events that are essential to our non-asymptotic analysis of the posterior distribution. They enable us to control the behavior of the error vector $\epsilon$. The proof is based on standard concentration inequalities for normal distributions.

Posterior Density of \texorpdfstring{$ (\gamma, \phi)$}{(gamma, phi)}

The presence of the simplex constraint makes the technical analysis substantially different from that for unconstrained variable selection models. The main challenge is to bound the marginal likelihood

To address this, we begin by deriving an equivalent expression for $p(y \mid \gamma, \mu_\gamma, \phi)$ using Theorem (ref). For each $\gamma \in \mathbb{S}_L$, define the OLS estimator for the model $\gamma$ by

Let $\check{\mu}_\gamma$ denote the least-squares estimator under the constraint $\mathsf{1}^\top \mu_\gamma = 1$. By Theorem (ref),

Further, for any $u \in \Delta^{|\gamma| - 1}$, we have

where

Equation (ref) decomposes the residual sum of squares for any $u$ satisfying the simplex constraint into three parts. The first term, $\| y - X_\gamma \hat{\mu}_\gamma \|_2^2$, is the residual sum of squares for the least-squares estimator in the unconstrained setting. The second term, $\| X_\gamma (u - \check{\mu}_\gamma) \|_2^2 $, measures the deviation of $u$ from the constrained least-squares estimator $\check{\mu}_\gamma$. Note that $X_\gamma \check{\mu}_\gamma$ is the projection of $X_\gamma \hat{\mu}_\gamma$ onto the affine space $\{ X_\gamma u \colon \mathsf{1}^\top u = 1 \}$; in particular, we always have $\mathsf{1}^\top \check{\mu}_\gamma = 1$. The last quantity, $b_\gamma \geq 0$, measures the deviation of $\check{\mu}_\gamma$ from $\hat{\mu}_\gamma$, which captures how well the model satisfies the simplex constraint. When $\hat{\mu}_\gamma \in \Delta^{|\gamma| - 1}$, we have $b_\gamma = 0$.

For each $\gamma$, we also define two matrices $D_\gamma \in \mathbb{R}^{|\gamma| \times (|\gamma| - 1) }$ and $\Pi_\gamma \in \mathbb{R}^{ (|\gamma| - 1) \times |\gamma| }$ by

Clearly $D_\gamma$ is a full-rank matrix whose column space is orthogonal to the one vector $\mathsf{1}$, and $\Pi_\gamma u$ yields the subvector of $u$ containing the first $|\gamma| - 1$ entries. Since the dimension of $D_\gamma$ or $\Pi_\gamma$ is always clear from context, we will omit the subscript and simply write $D$ and $\Pi$.

In Proposition (ref), we show that the decomposition given in (ref) enables us to compare the marginal likelihood $p(y \mid \gamma, \phi)$ with that in the unconstrained setting. The expectation with respect to $u\sim \text{Dir}(\mathsf{1})$ is converted to the probability mass (up to a normalizing constant) of a $(|\gamma| - 1)$-dimensional Gaussian distribution centered at $\Pi \check{\mu}_\gamma$.

Lower Bound on the Posterior Density of \texorpdfstring{$(\gamma^*, \phi)$}{(gamma*, phi)}

Next, we consider the true model $\gamma^*$. By Proposition (ref), to lower bound $p(\gamma^*, \phi \mid y)$, we need to find an upper bound on $b_\gamma^*$ and a lower bound on $\mathbb{P} ( \tilde{U}_{\gamma^*, \phi} \in {\tilde{\Delta}^{\ell^* - 1}} )$ where $ \tilde{U}_{\gamma^*, \phi} \sim \mathcal{N}( \Pi \, \check{\mu}_{\gamma^*}, (\phi D^\top X_{\gamma^*}^\top X_{\gamma^*} D)^{-1})$. Consider $b_\gamma^*$ first.

Bounding $\mathbb{P} ( \tilde{U}_{\gamma^*, \phi} \in {\tilde{\Delta}^{\ell^* - 1}} )$ is more difficult, and we divide the argument into two lemmas.

Proof of Theorem (ref)

Proof of Theorem (ref)

To prove Theorem (ref), we first use Theorem (ref) to identify the rate at which $p(\gamma \mid y)/p(\gamma^* \mid y)$ goes to zero for any single $\gamma \neq \gamma^*$. We consider two different cases separately: $ \gamma^* \subset \gamma$ (overfitted) and $\gamma^* \not\subset \gamma$ (underfitted).

Bounds for Overfitted Models

Bounds for Underfitted Models

Proof of Theorem (ref)

Proof of Theorem (ref)

Proof of Theorem (ref)

Consider the model $ Y = X w + \epsilon$ with $\epsilon \sim N(0, \phi^{-1} I )$. We use $\tilde{p}( \cdot \mid y)$ to denote the posterior distribution under the following prior:

This is the classical setup for the spike-and-slab Bayesian variable selection, except that we place an additional prior distribution on $\tau$ so that it aligns with our BVS-SS model.

Additional Simulation Studies

We present in the main text the simulation results for a sparse regression model, where each entry of $X$ is generated independently from a standard normal distribution. In this section, we provide additional simulations to evaluate the performance of BVS-SS in estimating the ATT, selecting relevant variables, and estimating variance components under two alternative data-generating processes.

First, in Section (ref), we consider a non-sparse setting where the data $(Y, X)$ are generated from a factor model. Second, in Section (ref), we examine a sparse model in which the covariates in $X$ are correlated. These simulation designs follow the factor-model setups commonly used in the literature hsiao2012panel, shi2023forward.

Non-sparse Factor Models

We first assume that the training data $X \in \mathbb{R}^{M \times N}$ and $Y \in \mathbb{R}^{M}$ are generated from a non-sparse factor model,

where $(Y, X) \in \mathbb{R}^{ M \times (N + 1) }$ denotes the block matrix obtained by concatenating $Y$ and $X$, $e$ is an $M\times (N+1)$ error matrix with each entry $e_{ij} \stackrel{\text{i.i.d.}}{\sim} N(0,0.5^2)$, $F \in \mathbb{R}^{4 \times M}$ is a factor matrix with 4 factors, and $\Lambda \in \mathbb{R}^{(N+1) \times 4}$ is the corresponding loading matrix. For the matrix $F$, we let the four common factors evolve according to

in which $u_{1,i}, u_{2,i}, u_{3,i}$ are i.i.d. standard normal. For $\Lambda$, we specify each loading $\lambda_{j,l}$ as

where $J$ is a parameter whose value will be specified in the simulation.

Then we simulate the data under treatment $\tilde{X} \in \mathbb{R}^{\tilde{M} \times N}$ and $\tilde{Y}^{(1)} \in \mathbb{R}^{\tilde{M}}$ by

where $\delta_i$ is defined as in Section (ref), and $\tilde{e} \in \mathbb{R}^{\tilde{M}\times (N+1)}$ are defined similarly. $\tilde{F} \in \mathbb{R}^{4 \times \tilde{M}}$ follows the same factor structure as $F$; that is, we generate $(F, \tilde{F})$ by (ref) for $i = 1, \dots, M+\tilde{M}$. All control units are correlated with the treated unit through the common factors. However, the dependence persists primarily for the first $J$ units, while the correlation for the remaining units diminish as the sample size grows.

As in the main text, we use OLS$(\gamma^*)$ as the reference, and we report the relative efficiencies of the other five estimators in Figure (ref). Under a factor structure, the simplex constraint is typically violated, and thus the performances of two QP estimators deteriorate due to the violation of the simplex constraint. BVS-SS outperforms Lasso especially when the ratio $N/M$ is large and quickly approaches the performance of OLS$(\gamma^*)$. Note that OLS also approaches the benchmark OLS($\gamma^*$) as the sample size increases when $N=20$, consistent with the finding that OLS is a reasonable tool for the first stage estimation in SCM hsiao2012panel. However, OLS fails to provide a unique estimation when the sample size is small (e.g., $M=25, N=50$). These results reinforce the key message in Section (ref) in the main text: the robust and superior performance of BVS-SS is largely driven by the use of the soft simplex constraint, which provides substantial gains even when the true DGP follows a non-sparse and unknown factor structure.

Meanwhile, the posterior mean of $\phi$ does not exhibit a clear increasing trend. This is due to our simulation setting, in which the true DGP is a factor model while our Bayesian estimation is based on a sparse regression model. The trend of $\tau$ remains relatively stable, exhibiting only a very slight decrease as $M$ increases. Moreover, the magnitude of the posterior mean of $\tau$ is very close to that in the sparse regression model with $||w^*||_1=2$ presented in Section (ref). This further supports the notion that the factor model is likely to violate the simplex constraint.

Sparse Factor Models

Next, we consider a sparse model in which $X\in \mathbb{R}^{M \times N}$ is generated from a factor model with factor matrix $F \in \mathbb{R}^{4 \times M}$ and loading matrix $\Lambda \in \mathbb{R}^{N \times 4}$, and the untreated potential outcome for $Y \in \mathbb{R}^{M \times 1}$ is a linear combination of $X$,

The treatment and sparse coefficients $w^*$ are specified in the same way as in Section (ref). The sparsity still remains through $w^*$, but all units in the control group are linked by some factors that do not directly affect the treated unit. Analogously, we generate the following post-treatment data,

where $\tilde{Y}^{(1)} \in \mathbb{R}^{\tilde{M} \times 1}$, $\tilde{X} \in \mathbb{R}^{\tilde{M} \times N}$, $\tilde{F} \in \mathbb{R}^{4 \times \tilde{M}}$. $\delta \in \mathbb{R}^{\tilde{M}}$ and $\tilde{\epsilon} \in \mathbb{R}^{\tilde{M}}$ are defined as in Section (ref). The specifications of $(F,\tilde{F})$ and $\Lambda$ are the same as in Section (ref), except that $\Lambda$ contains $N$ rows rather than $N+1$.

We report the relative efficiencies of the five estimators in Figure (ref) with OLS$(\gamma^*)$ fixed at 100%. The QP estimators behave similarly to that in the main text, outperforming BVS-SS when the simplex constraint is satisfied ($\| w^* \|_1=1$), and performing poorly when $\| w^* \|_1=2$ or $3$. BVS-SS still outperforms Lasso when the ratio $N/M$ is large in all scenarios, though Lasso and OLS estimators approach the orcale OLS($\gamma^*$) slightly better when $M=200$. These results show that our method can still provide robust and superior performance for ATT estimation when the control units are correlated.

Regarding variable selection, BVS-SS generally attains higher accuracy than Lasso, particularly when the sample size $M$ is large. In smaller samples, however, its model selection performance deteriorates. For example, when $N=50, J=10, M\leq 50$, BVS-SS fails to select any true control unit on average, although Lasso performs nearly as poorly under these small-sample scenarios. This pattern is understandable: because our sampler is a Gibbs sampler that updates coefficients pairwise, strong correlations among covariates can lead to highly correlated variables being selected in place of the true signals. Importantly, despite the low accuracy of variable selection, BVS-SS delivers substantially better ATT estimation than Lasso when $M \leq 50$, as shown in Figure (ref), highlighting its robustness even when exact model recovery is difficult. This observation also underscores the need for more advanced MCMC algorithms which might improve selection accuracy in such correlated designs.

Finally, the posterior means of $\tau$ and $\phi$ exhibit patterns very similar to those reported in the main text, and we visualize their distributions in Figures (ref) and (ref).

Additional Results for Empirical Examples

This section presents supplementary details for the empirical examples analyzed in the main text.

Additional Results for NFP Anti-tax Evasion

We present more results for the NFP anti-tax evasion study in this section. Figure (ref) shows the trajectories of the posterior samples of $\phi$ and $\log(\tau)$ (with burn-in included). Probably due to the limited number of observations ($T_0=33$), the MCMC samples of $\tau$ exhibit large variability, indicating that the simplex assumption is hard to confirm or reject, though the ratio of the posterior means of $\tau$ and $\phi$ given in Table (ref) is relatively small. The trace plot suggests that our sampler has converged fast.

The posterior distributions of ATT, $\log(\tau)$, $\phi$ and model size $|\gamma|$ are visualized in the histograms in Figure (ref). The posterior mean of the ATT exhibits a clear unimodal and bell-shaped distribution, concentrated in the interval $[0.26,0.34]$. The distribution of $\phi$ is symmetric with a mode between $[3.5,4.0]$, consistent with the posterior mean in Table (ref), while the distribution of $\log(\tau)$ is right-skewed, indicating that a large portion of the posterior draws for $\tau$ are not concentrated near zero. This pattern is also reflected in the posterior mean of $\tau$ reported in Table (ref). Taken together, the evidence suggests that the data in this empirical example do not favor the simplex constraint. The posterior distribution of $|\gamma|$ is left-skewed, with more than half of the posterior mass supporting a relatively small model that includes only two control units.

Finally, Figure (ref) plots the observed FHA inflation index together with the estimated counterfactual path. Before the intervention (marked by the vertical green line), the estimated pre-treatment counterfactual (blue dashed line) tracks the observed FHA index well, indicating a good in-sample fit. Following the intervention, the post-treatment counterfactual (red dashed line) diverges from the observed outcomes, capturing the estimated treatment effect over time.

Additional Results for China's Anti-corruption Campaign

The section provides more details on the China's anti-corruption example. From Figure (ref), one can observe that the posterior distribution of $\tau$ is still dispersed, but that of $\phi$ is more concentrated and much larger on average, compared to the previous example. The estimation of $\phi$ appears to be more accurate, suggesting that the error variance is around 0.05 (1/20.86). The ratio $\tau/\phi$ becomes much smaller, indicating that the simplex constraint is more likely to be satisfied (at least approximately) in this data set.

Figure (ref) presents the posterior distributions. The ATT distribution is unimodal but a little right-skewed. The distribution of $\log(\tau)$ is in line with its trace plot, with a large probability mass for $\log(\tau) < -2$, but the distribution appears relatively diffuse since the bars in the histogram are fairly evenly distributed. More than 80% samples of $\phi$ concentrate in the interval of $[15,30]$, providing evidence that our BVS-SS has detected abundant signals in the data. BVS-SS still suggests a small model size for building the counterfactual, as around 70% MCMC samples have model size smaller than six.

Finally, we provide a counterfactual plot for the growth rate of the luxury watches import in Figure (ref). Compared with the previous example, the in-sample fit is noticeably better, likely due to the availability of a richer set of control units. A pronounced initial decline in the import growth rate is detected in the post-intervention periods, which implies a substantial immediate policy impact. However, the effect attenuates over time, with the gap between the observed and counterfactual paths gradually narrowing. Overall, the figure illustrates that our BVS-SS reproduces the pre-treatment dynamics effectively and provides a coherent counterfactual trajectory to assess the policy impact.