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.
85,077 characters · 35 sections · 37 citation commands
Recovering Counterfactual Distributions via Wasserstein GANs
Since its inception by abadie2003economic and subsequent formalization by abadie2010synthetic, the Synthetic Control Method has emerged as the canonical framework for comparative case studies in empirical economics. By constructing a counterfactual unit as a convex combination of untreated donors, the method provides a transparent, data-driven mechanism to estimate treatment effects. While the standard framework focuses exclusively on mean outcomes, recent policy evaluations increasingly demand an understanding of distributional impacts. Whether assessing the effect of minimum wage laws on income inequality or the impact of trade liberalization on firm productivity dispersion, structural interest often lies in the tails or the shape of the distribution rather than its first moment. Addressing this imperative, recent scholarship has advanced the concept of DSC to recover the full counterfactual probability law of the treated unit gunsilius2023distributional.
Despite these conceptual advances, existing implementations remain vulnerable to the limitations of the Euclidean metrics they employ. Standard approaches typically estimate counterfactuals by minimizing the integrated squared error between the quantile functions of the treated and synthetic units. This strategy relies on pointwise vertical comparisons and implicitly assumes that the counterfactual distribution is a location-scale transformation of the donors. Such a restriction is non-trivial; it imposes a form of rank invariance that prevents the estimator from recovering complex morphological changes, such as a unimodal distribution evolving into a bimodal one. Furthermore, in settings characterized by support mismatch—where the outcomes of the treated unit lie outside the support of the donor distributions—the squared error objective becomes locally constant. This lack of informative gradients renders gradient-based optimization ineffective and leads to severe estimator instability.
Our methodological departure from Euclidean metrics to OT is primarily motivated by a pervasive challenge in financial econometrics: characterizing the geometry of risk under structural breaks. Financial time series—unlike smooth macroeconomic aggregates—are frequently characterized by heavy-tailed contamination, disjoint supports due to liquidity gaps, and complex multimodal distributions during regime switches. In such regimes, standard distributional estimators that rely on quantile averaging or probability density matching often yield vanishing gradients or spurious unimodal artifacts, failing to capture the true latent factor structure. By reformulating the Synthetic Control problem as a Wasserstein minimization via a WGAN, we provide a robust estimator capable of recovering counterfactual probability laws even when the post-shock distribution shares no common support with the pre-shock history. While we demonstrate this capability through a canonical policy experiment—the Kansas tax cuts, which serves as an ideal stress test involving a massive, exogenous distributional shift—the proposed framework is broadly applicable to asset pricing and risk management, particularly for recovering counterfactual return distributions and tail dependencies in the presence of idiosyncratic shocks.
To resolve these impediments, we propose a robust estimator grounded in the theory of OT. We identify synthetic control weights by minimizing the Wasserstein-1 distance between the treated and synthetic distributions. The departure from Euclidean metrics to the Wasserstein distance represents a fundamental shift in how discrepancies are quantified. While Euclidean metrics penalize deviations at fixed coordinate points, the Wasserstein distance measures the minimum cost required to transport probability mass from one distribution to another. This metric respects the underlying geometry of the outcome space and induces a topology that preserves gradient information even when distributions have disjoint supports.
The economic intuition behind this shift is straightforward. Consider a scenario where a policy shifts the entire outcome distribution of a treated unit by a constant amount relative to a control unit. A Euclidean metric comparing densities or quantile functions at fixed points would register this merely as a maximum error, failing to capture the proximity of the two distributions. In contrast, the Wasserstein distance linearly reflects the magnitude of this shift, correctly identifying the structural similarity between the units. This property allows our estimator to effectively search for donors that are structurally closest to the treated unit based on the geometry of the data, rather than failing due to a lack of local density overlap. Consequently, our approach relaxes the rigidity of quantile averaging and permits the construction of synthetic counterfactuals that exhibit complex shapes, such as mixtures, which are not present in any single donor unit.
Our work contributes to three distinct strands of literature. First, we extend the methodological frontier of the SCM. Following the seminal work of abadie2010synthetic, numerous extensions have been proposed to improve weight selection, including penalized estimators and matrix completion methods athey2021matrix, ben2021augmented. While these innovations enhance robustness for mean outcomes, they do not address the estimation of full distributions. Second, we connect with the growing literature on OT in economics. galichon2016optimal demonstrated the power of OT for structural estimation in matching models, yet its application to causal inference and time-series counterfactuals remains nascent. We bridge this gap by showing that the geometry induced by OT is uniquely suited for the Synthetic Control framework. Third, we add to the integration of machine learning into econometrics. Unlike approaches that utilize machine learning purely for prediction chernozhukov2018double, we employ deep learning architectures to solve a structural minimization problem defined by a theoretically grounded statistical metric.
The primary theoretical contributions of this paper are the asymptotic theory and a formal identification result for synthetic controls in the space of probability measures. We demonstrate that the validity of our estimator does not require the restrictive rank-invariance assumptions inherent in quantile-based methods. Instead, we posit a latent factor model where the evolution of the outcome distribution is governed by unobserved time-varying common factors and unit-specific factor loadings. We show that under a scaled isometry condition, a synthetic control that successfully matches the pre-treatment distribution of the treated unit in terms of the Wasserstein distance effectively recovers the underlying factor loadings. This result generalizes the classic parallel trends assumption to the Wasserstein space, providing a rigorous justification for using linear mixtures of donors to approximate the counterfactual probability law of the treated unit.
Computationally, estimating the Wasserstein distance in high-dimensional or continuous spaces presents significant challenges, often requiring linear programming solutions that are cubically expensive in the sample size. To render this estimation feasible, we introduce a computational framework based on the Kantorovich-Rubinstein duality. Rather than solving the primal transport problem directly, we estimate the dual potentials using a WGAN. In this architecture, a critic network approximates the Lipschitz-continuous function that maximizes the transport cost, while a generator network—representing the synthetic weights—minimizes it. This adversarial formulation bypasses the sorting bottlenecks of traditional quantile methods and allows for the scalable, differentiable estimation of weights even when the outcome variables are multivariate. By solving the dual problem, our estimator robustly handles continuous supports and provides stable gradients, ensuring convergence even in regimes where standard density-matching techniques fail.
We corroborate our theoretical and computational arguments through a series of Monte Carlo simulations designed to stress-test the estimator under pathological conditions often encountered in economic data. Our results demonstrate that the WGAN-based solver significantly outperforms standard $L_2$ benchmarks in terms of both stability and accuracy. Specifically, in scenarios involving heavy-tailed data contamination and severe support mismatch, the standard quantile-based estimator exhibits substantial variance and bias, whereas our Wasserstein-based method remains robust. We further distinguish our approach by its ability to handle structural complexity. In a simulation where the true counterfactual involves a transition from a unimodal to a bimodal distribution—mimicking phenomena such as labor market polarization—we show that while quantile-averaging methods structurally fail and yield a unimodal artifact, our mixture-based approach successfully recovers the distinct modes. Additionally, we provide a proof-of-concept demonstration in a multivariate setting, confirming that our adversarial framework scales effectively to higher-dimensional outcomes where traditional sorting-based quantile methods are computationally intractable.
The remainder of the paper is organized as follows. Section 2 reviews the relevant literature. Section 3 outlines the econometric framework and defines the notation. Section 4 details the estimation strategy, introducing the WGAN architecture and the optimization algorithm based on the dual formulation of OT. Section 5 establishes the asymptotic theory and the formal identification result. Section 6 reports the results of our Monte Carlo simulations to evaluate finite-sample performance. Section 7 presents an empirical application. Section 8 concludes. All mathematical proofs and implementation details are provided in the Appendix.
Since its formalization by abadie2010synthetic, the Synthetic Control Method has become a cornerstone of comparative case studies in empirical economics. The framework's ability to construct data-driven counterfactuals has catalyzed a rich methodological literature aimed at enhancing the robustness of average treatment effect estimates. Significant extensions include the introduction of penalized estimators to mitigate overfitting in high-dimensional donor pools abadie2021penalized and the augmented synthetic control method, which corrects for regularization bias ben2021augmented. Further integrating the method with modern causal inference tools, arkhangelsky2021synthetic developed synthetic difference-in-differences to leverage both unit and time fixed effects, while athey2021matrix reformulated the counterfactual imputation problem via matrix completion. Collectively, these contributions have substantially refined the estimation of conditional means $\mathbb{E}[Y]$, establishing a rigorous standard for aggregate policy evaluation.
Complementing these aggregate approaches, a growing literature emphasizes the necessity of evaluating distributional impacts. As articulated in the program evaluation literature by bitler2006mean and firpo2009unconditional, structural parameters of interest—such as inequality measures or tail risks—often depend on the full probability law rather than the first moment. Addressing this imperative, gunsilius2023distributional pioneered the extension of SCM to the space of probability measures. This framework constructs the counterfactual distribution by minimizing the $L_2$ distance between the quantile function of the treated unit and a weighted average of donor quantile functions. Subsequent work has further solidified this quantile-based paradigm: zhang2024asymptotic established the asymptotic optimality of such estimators under general regularity conditions, while kato2023asymptotically introduced distribution-matching variations designed to enhance finite-sample performance. These contributions have successfully generalized the logic of synthetic controls to distributional objects under the assumption that the counterfactual preserves the rank structure of the donors.
Despite these theoretical strides, existing implementations remain vulnerable to data imperfections common in empirical settings. A primary concern arises from the reliance on squared-error loss functions. As noted in the robust statistics literature huber1992robust, such metrics are highly sensitive to contamination, implying that tail outliers can disproportionately distort weight estimation. A second, more fundamental challenge pertains to the geometry of the support. The validity of standard quantile averaging is predicated on the assumption that the target distribution lies within the convex hull of the donors. In scenarios characterized by support mismatch—where the target distribution is disjoint from the donors—metrics based on pointwise vertical discrepancies fail to capture the underlying proximity, often resulting in uninformative gradients. These geometric limitations suggest the need for a metric that respects the topology of the outcome space, motivating our departure from Euclidean distances.
Our methodology aligns with the geometric turn in recent econometric theory, characterized by the increasing integration of OT principles into structural modeling. While classically rooted in pure mathematics, OT was systematically introduced to economics by galichon2016optimal, primarily as a rigorous tool for analyzing matching markets and hedonic models. The central object of this theory—the Wasserstein distance—offers a metric for probability measures that explicitly respects the underlying geometry of the outcome space. This geometric sensitivity provides a decisive advantage in settings where standard information-theoretic divergences, such as the Kullback-Leibler divergence or the Total Variation distance, fail to provide meaningful comparisons for distributions with disjoint supports. Consequently, the Wasserstein metric has proven instrumental in facilitating advances across a spectrum of econometric challenges, ranging from distributionally robust optimization blanchet2019quantifying to partial identification in incomplete models.
More pertinent to our framework is the subsequent shift from employing OT solely for identification to utilizing it for statistical inference. A growing literature on Wasserstein Minimum Distance Estimation (W-MDE) demonstrates that matching model-implied distributions to data via transport costs yields consistent estimates, offering a robust alternative in settings where the likelihood function is intractable or singular bernton2019parameter. This inference paradigm shares a fundamental conceptual link with Generative Adversarial Networks (GANs). In the econometric domain, this connection has been formalized by kaji2023adversarial, who establish that adversarial estimators can be interpreted as minimizing an Integral Probability Metric (IPM) defined by a discriminator class. From a statistical perspective, liang2021how further provides minimax rates of convergence for distribution estimation under this framework, offering theoretical guarantees for non-parametric estimation. Collectively, these results provide the necessary foundation for utilizing neural networks in causal inference: rather than functioning as opaque “black box” generators, the adversarial architecture serves as a flexible, data-driven device for moment selection, ensuring that the synthetic controls satisfy the geometric constraints of the outcome space.
Despite these significant strides in cross-sectional distribution estimation, the application of OT-based adversarial inference to panel data and dynamic causal frameworks remains underexplored. While gunsilius2023distributional explicitly acknowledges the theoretical appeal of Wasserstein metrics for synthetic controls—noting their ability to capture geometric proximity beyond rank invariance—no feasible computational strategy has yet been proposed to operationalize this concept for trajectory reweighting. Existing extensions of SCM that incorporate machine learning, such as athey2021matrix, predominantly focus on low-rank approximations for mean outcomes rather than recovering full distributional counterfactuals via geometric transport. This paper bridges this specific gap. By adapting the Wasserstein GAN architecture to the unique structure of the synthetic control problem, we leverage the geometric robustness of OT to resolve the support mismatch failure that limits existing quantile-based estimators. This contribution effectively extends the practical reach of adversarial inference from static density estimation to the dynamic, counterfactual setting of comparative case studies.
Finally, our work contributes to the rapidly expanding interface between deep learning and causal inference. The integration of machine learning into econometrics has historically proceeded through two distinct waves. The first wave leveraged regularization techniques to manage high-dimensional controls and improve predictive accuracy belloni2014high, mullainathan2017does. This foundation paved the way for a second wave focused on valid statistical inference, best exemplified by the Double Machine Learning framework of chernozhukov2018double and the Generalized Random Forests of wager2018estimation. These methods utilize flexible function approximators to handle nuisance parameters, thereby enabling the unbiased estimation of average and heterogeneous treatment effects. However, these canonical approaches operate primarily within the paradigm of discriminative learning, targeting conditional expectations $\mathbb{E}[Y|X]$ while leaving the modeling of complex, high-dimensional counterfactual distributions largely unaddressed.
To capture distributional heterogeneity beyond the conditional mean, recent scholarship has increasingly turned to generative modeling. Unlike discriminative approaches, generative frameworks aim to learn the underlying probability density of the potential outcomes. Notable contributions in this vein include the GANITE framework of yoon2018ganite, which employs GANs to impute missing counterfactuals in cross-sectional data, and hartford2017deep, who introduced Deep IV to address endogeneity via deep neural networks. While these studies demonstrate that adversarial training can effectively recover complex data-generating processes, the application of deep generative models to the Synthetic Control framework remains nascent. Current machine learning extensions of SCM predominantly focus on low-rank approximations; for instance, athey2021matrix formulate the counterfactual estimation as a matrix completion problem using nuclear norm regularization. While powerful, such methods rely heavily on the assumption that the matrix of potential outcomes possesses a low-rank structure. They do not explicitly model the geometric transport of probability mass between distributions. Our work bridges this divide by proposing a WGAN-based solver that moves beyond the low-rank assumption, offering a methodology capable of preserving the distributional properties of the trajectory even when the data exhibits complex, non-linear structures.
This section formalizes the problem of recovering counterfactual distributions. We depart from Euclidean metrics on quantile functions gunsilius2023distributional by grounding our framework in the geometry of OT. This shift enables causal identification even when the support of the target distribution is disjoint from the donor pool, while accommodating granular, micro-level data structures.
Consider a panel of $J+1$ units observed over $T$ periods. Let unit $i=1$ denote the treated unit exposed to an intervention after period $T_0$, where $1 \le T_0 < T$. The set $\mathcal{D} = \{2, \dots, J+1\}$ serves as the donor pool of units that remain untreated throughout the observation window.
Unlike standard synthetic control methods that focus on point estimates $\mathbb{E}[Y_{it}]$, our object of interest is the full probability law of the outcome. We utilize the full granular information available at the micro-level. For each unit $i$ and period $t$, we observe a sample of $M_{it}$ independent observations, denoted by $\{y_{it}^{(k)}\}_{k=1}^{M_{it}}$. We define the observed outcome $P_{it} \in \mathcal{P}_1(\mathcal{X})$ as the empirical measure:
where $\delta$ denotes the Dirac mass. Crucially, we do not require the sample size $M_{it}$ to be balanced. The number of micro-observations can vary arbitrarily across units and time periods, allowing for the flexible integration of unbalanced survey data.
For the post-intervention period $t > T_0$, we observe the treated distribution $P_{1t} = \mathcal{L}(Y_{1t}(1))$, but we seek to estimate the counterfactual distribution $P_{1t}^N = \mathcal{L}(Y_{1t}(0))$. We construct a synthetic control as a convex combination of the donor distributions. Let $\Lambda = \Delta^{J}$ denote the simplex defined as follows:
For any weight vector $\lambda \in \Lambda$, we define the synthetic distribution $P_{\lambda, t}$ as the linear mixture of donor measures:
This formulation imposes no parametric constraints on the shape of the counterfactual, allowing for the recovery of complex, multimodal distributions.
The core innovation of our framework lies in the choice of the discrepancy metric used to determine the optimal weights. While standard methods typically minimize the integrated squared error between quantile functions, we adopt the Wasserstein-1 distance to ensure robustness against support mismatch.
Formally, utilizing the Kantorovich-Rubinstein duality, we define the transport cost between the treated unit and a synthetic mixture at any specific time period $t$ as:
where $Lip_{1}(\mathcal{X})$ denotes the set of 1-Lipschitz continuous functions. Unlike Euclidean metrics, this distance provides valid gradient signals even when the supports of the distributions are disjoint, as the gradient is governed by the underlying transport cost rather than pointwise density overlaps.
To construct the final synthetic control, our estimator aggregates information across the pre-intervention window. We adopt a period-wise estimation strategy to capture the structural fit at each time step. For each pre-intervention period $t \le T_0$, we define the period-specific optimal weights $\lambda_{t}^{*}$ as the minimizer of the regularized Wasserstein loss:
where $H(\lambda) = - \sum_{j=2}^{J+1} \lambda_j \log \lambda_j$ is the entropic regularization term. The parameter $\eta > 0$ controls the trade-off between fitting the transport cost and maximizing the entropy of the weights, ensuring solution uniqueness and numerical stability. Note that while the objective function varies with $t$ through the observed distributions, the regularization term $H(\lambda)$ depends solely on the weight vector and is time-invariant.
The final synthetic control weights $\lambda^{*}$ are obtained by aggregating these period-specific estimates using temporal weights $\{w_{t}\}_{t=1}^{T_{0}}$ (typically $w_{t} = 1/T_{0}$ for a simple average):
This two-step aggregation procedure stabilizes the estimation by reducing the variance associated with any single period's idiosyncratic shocks, while ensuring that the final weights reflect a consistent structural similarity across the entire pre-intervention history.\footnote{We employ a linear aggregation of period-specific weights rather than a single global optimization to accommodate potential structural breaks or heteroskedasticity in the pre-treatment period. However, alternative aggregation schemes, such as weighting periods by their inverse Wasserstein loss, are straightforward extensions of this framework.}
We now establish the conditions under which minimizing the pre-intervention discrepancy identified in Eq. (ref) recovers the true counterfactual distribution. Following the literature on synthetic controls, we assume the outcomes are generated by a latent factor structure, which we extend to the space of probability measures.
Remark. Assumption 2 implies that while aggregate shocks may alter the location, scale, or shape of outcomes (scaling distances by $\kappa_t$), they preserve the relative structural geometry between units. This condition serves as a distributional generalization of the parallel trends assumption.
Based on this structural geometry, Theorem 1 establishes that a synthetic control constructed to match the pre-treatment observables identifies the latent counterfactual.
\paragraph{Theorem 1 (Counterfactual Recovery).} Suppose Assumptions 1 and 2 hold. Let $\lambda^{*} \in \Lambda$ be a weight vector satisfying the following conditions:
Then, $\lambda^{*}$ perfectly recovers the counterfactual distribution in the post-intervention period:
\paragraph{Proof.} See Appendix A.
A critical theoretical advantage of our framework over standard distributional synthetic controls lies in its robustness to support mismatch and model misspecification. We analyze two distinct geometric failure modes of existing estimators and demonstrate how the Wasserstein metric resolves them.
\paragraph{Gradient Existence under Support Mismatch.} A key computational advantage of the Wasserstein metric is the preservation of informative gradients even when distributions have disjoint supports. This property explains the stability of our estimator in Simulation 2, where $L_2$-based methods fail. We formalize this property in the following proposition.
\paragraph{Proof.} See Appendix B.
Beyond support mismatch, our framework also addresses structural rigidity in the shape of the counterfactual. Existing approaches gunsilius2023distributional typically define the synthetic control as a weighted average of quantile functions: $Q_{synth}(\tau) = \sum \lambda_j Q_{j}(\tau)$. This formulation implicitly assumes that the counterfactual distribution belongs to a location-scale family generated by the donors. Such an assumption imposes a rigid unimodal constraint: an average of unimodal quantile functions typically results in a unimodal distribution. Our framework, by defining the synthetic control as a linear mixture of measures ($P_{synth} = \sum \lambda_j P_{j}$), relaxes this restriction. As demonstrated in our bimodal simulation (Simulation 3), a mixture of unimodal donors can naturally approximate a multimodal target. Thus, our estimator is robust to misspecification regarding the shape of the counterfactual distribution.
While the identification results in the previous section rely on population distributions, in practice we observe only finite samples. This section details the estimation procedure, replacing the intractable supremum over Lipschitz functions with a neural network approximation, and outlines the adversarial optimization framework used to solve for the synthetic weights.
Let $\{y_{it}^{(k)}\}_{k=1}^{M_{it}}$ denote the observed micro-samples for unit $i$ at time $t$. We approximate the underlying population distribution $P_{it}$ using the empirical measure $\hat{P}_{it}$:
where $\delta$ is the Dirac mass.
Directly computing the Wasserstein distance $W_1(\hat{P}_{\lambda, t}, \hat{P}_{1t})$ involves solving a high-dimensional linear programming problem, which becomes computationally prohibitive when the sample size $M_{it}$ is large. Moreover, we require differentiable gradients with respect to the weights $\lambda$ to optimize the synthetic control. To address these challenges, we adopt the WGAN framework arjovsky2017wasserstein.
We parameterize the Kantorovich potential, which is the critic function, using a neural network $f_{\theta}: \mathcal{X} \to \mathbb{R}$ with parameters $\theta \in \Theta$, where $\Theta$ is a compact parameter space. The space of 1-Lipschitz functions is thus approximated by the parameterized family $\mathcal{F}_{\Theta} = \{f_{\theta} : \theta \in \Theta\}$. Consequently, the estimation of the period-specific weights $\lambda_{t}^{*}$ is recast as the solution to the following empirical minimax problem:
where the regularized empirical objective function $\mathcal{L}_{t}$ is given by:
In this formulation, the expectations over empirical measures are computed as sample means (e.g., $\mathbb{E}_{x \sim \hat{P}_{jt}}[f_{\theta}(x)] = \frac{1}{M_{jt}} \sum_{k} f_{\theta}(y_{jt}^{(k)})$). The inner maximization over $\theta$ approximates the Wasserstein distance, which is finding the tightest lower bound on the transport cost, while the outer minimization over $\lambda$ finds the optimal weights that minimize this estimated distance plus the entropic penalty.
Solving the minimax problem in Eq. (ref) requires the critic network $f_{\theta}$ to belong to the set of 1-Lipschitz functions, a necessary condition for the Kantorovich-Rubinstein duality to hold. Standard weight clipping techniques often lead to optimization pathologies. To strictly enforce this constraint, we employ the Gradient Penalty method proposed by gulrajani2017improved.
We augment the objective function of the critic with a regularization term that penalizes deviations of the gradient norm from unity. The modified loss for the critic becomes:
where $\zeta > 0$ is the penalty coefficient, and $\hat{P}_{int}$ denotes the distribution of random interpolates sampled uniformly along straight lines connecting points from the treated and synthetic distributions. This soft constraint ensures that $f_{\theta}$ approximates the Wasserstein distance with valid gradients almost everywhere.
We solve the saddle-point problem using an alternating Gradient Descent-Ascent (GDA) procedure. The optimization involves two distinct update rules tailored to the geometry of the parameters.
For a fixed weight vector $\lambda$, we update the network parameters $\theta$ to maximize the separation between the treated and synthetic distributions. We apply standard gradient ascent or Adam on the penalized objective $\mathcal{L}_{critic}$.
For a fixed critic $f_{\theta}$, the minimization with respect to $\lambda$ is strictly convex due to the entropic regularization. However, standard gradient descent is inefficient and may violate the simplex constraint $\sum \lambda_{j}=1$. Instead, we employ Exponentiated Gradient Descent, which naturally preserves the simplex geometry and ensures strictly positive weights. The multiplicative update rule for each donor $j \in \{2, \dots, J+1\}$ is derived from the first-order optimality condition:
where $\alpha_{\lambda}$ is the learning rate and $Z^{(k+1)}$ is the normalization factor ensuring the weights sum to one.
The complete estimation procedure, incorporating the period-wise strategy and temporal aggregation, is detailed in Algorithm 1.
The estimation of the optimal weights $\lambda^*$ involves solving a minimax optimization problem. Based on the identification results in Equation ((ref)), the estimator is the solution to:
where $\mathcal{L}(\lambda, f) = \mathbb{E}_{P_1}[f] - \sum_{j=2}^{J+1} \lambda_j \mathbb{E}_{P_j}[f]$.
We propose a stochastic approximation algorithm that alternates between optimizing the critic function $f$ and the synthetic weights $\lambda$.
We parameterize the Kantorovich potential $f$ using a neural network $f_\theta: \mathcal{X} \to \mathbb{R}$ with parameters $\theta \in \Theta$. The network architecture consists of a Multi-Layer Perceptron (MLP) with Leaky ReLU activation functions to ensure efficient gradient flow.
To enforce the 1-Lipschitz constraint required by the Wasserstein duality, we adopt the Gradient Penalty (GP) method proposed by gulrajani2017improved. Unlike weight clipping, which biases the critic towards simple functions, GP enforces the Lipschitz constraint softly by penalizing the norm of the gradient. The empirical objective for the critic becomes:
where $P_{\hat{x}}$ is defined by sampling uniformly along straight lines between pairs of points from the target and donor distributions, and $\gamma_{GP} > 0$ is the penalty coefficient\footnote{We conduct a comprehensive sensitivity analysis in Appendix D, demonstrating that the estimator's performance remains robust to variations in key hyperparameters, including the gradient penalty coefficient and network width.}.
For the outer minimization problem with respect to $\lambda$, we must handle the simplex constraint $\lambda \in \Delta^{J-1}$ and the entropic regularization term. Standard Projected Gradient Descent is suboptimal here as it ignores the geometry induced by the entropy.
Instead, we employ Entropic Mirror Descent (also known as the Exponentiated Gradient algorithm), which is the natural first-order method for entropy-regularized problems. The update rule is multiplicative, ensuring that $\lambda$ strictly remains in the positive simplex without explicit projection steps.
At iteration $t$, given the current critic $f_{\theta^{(t)}}$, the gradient of the transport cost with respect to $\lambda_j$ is simply the negative expected value of the critic for donor $j$:
The weight update follows the mirror descent step with learning rate $\alpha_\lambda$:
In practice, since the regularization term $\eta \lambda \log \lambda$ is strictly convex, we can also decouple the update: update $\lambda$ based on the transport gradient $g^{(t)}$ and then apply the proximal operator for the entropy, which simplifies to the softmax-like scaling.
The complete training procedure is detailed in Algorithm (ref). Since the sample size $N$ in synthetic control applications is typically moderate (e.g., hundreds of time periods), we can often use full-batch gradient descent rather than mini-batches to reduce variance.
In this section, we establish the statistical identification of the synthetic control weights. While Section 3 establishes the causal identification, for example, the recovery of the unobserved counterfactual distribution under structural assumptions, this section focuses on the statistical identification of the weight vector $\lambda^*$ itself. Specifically, we establish conditions under which the regularized minimization problem admits a unique solution, ensuring that the estimator is well-defined.
Let $(\mathcal{X}, d)$ be a complete separable metric space. Let $\mathcal{M}(\mathcal{X})$ denote the vector space of finite signed Borel measures on $\mathcal{X}$, and let $\mathcal{P}_1(\mathcal{X}) \subset \mathcal{M}(\mathcal{X})$ denote the space of Borel probability measures with finite first moments. We endow $\mathcal{P}_1(\mathcal{X})$ with the 1-Wasserstein distance $W_1$.
We observe the target distribution $P_1 \in \mathcal{P}_1(\mathcal{X})$ and a set of donor distributions $\mathcal{D} = \{P_2, \dots, P_{J+1}\} \subset \mathcal{P}_1(\mathcal{X})$. The synthetic control is defined as a linear mixture $P_\lambda = \sum_{j=2}^{J+1} \lambda_j P_j$, where the weights $\lambda$ reside in the simplex:
The estimator $\lambda^*$ is defined as the minimizer of the regularized transport cost:
To ensure the mapping from weights to distributions is injective, we impose the following condition:
We examine the asymptotic properties of the estimator in the regime of Large Micro-Sample, Fixed Donors. Specifically, we consider the case where the number of micro-level observations per period, denoted by $N$ (or $M_{it}$), tends to infinity, while the number of donor units $J$ remains fixed. This regime is relevant for distributional synthetic controls where researchers have access to granular data (e.g., census or survey micro-data) for each aggregate unit.
Let $\{Y_{1}^{(k)}\}_{k=1}^N$ be $N$ i.i.d. observations from the target $P_1$, and $\{Y_{j}^{(k)}\}_{k=1}^N$ be $N$ i.i.d. observations from donor $P_j$ for $j \in \{2, \dots, J+1\}$. The empirical measures are $\hat{P}_{1,N} = \frac{1}{N} \sum_{k=1}^N \delta_{Y_{1}^{(k)}}$ and $\hat{P}_{j,N} = \frac{1}{N} \sum_{k=1}^N \delta_{Y_{j}^{(k)}}$.
The estimator minimizes the empirical regularized objective:
To derive the limiting distribution of $\hat{\lambda}_N$, we utilize the theory of Sieve M-estimation (Chen, 2007). We treat the critic network $f_\theta$ as a sieve approximator for the dual potential space $Lip_1(\mathcal{X})$.
Let $\mathcal{F}_N = \{f_\theta : \theta \in \Theta_N\}$ be the class of functions parameterized by the neural network, where the complexity of the network (e.g., width and depth) may grow with $N$. Let $\lambda^*$ be the population minimizer and $f^*(\cdot, \lambda)$ be the optimal Kantorovich potential for a given $\lambda$.
We impose the following regularity conditions:
While Theorem (ref) establishes that the WGAN-based estimator achieves the parametric $\sqrt{N}$ convergence rate asymptotically, relying on the Gaussian approximation for hypothesis testing presents practical challenges in the context of Synthetic Controls.
First, the number of donor units $J$ is typically small, rendering the estimation of the asymptotic variance matrix $\Sigma = (H^*)^{-1} V (H^*)^{-1}$ unstable. Second, and more critically, the optimal weights $\lambda^*$ often lie on the boundary of the simplex $\Delta^{J-1}$ (i.e., sparsity). In such cases, the finite-sample distribution of the estimator is often non-standard and bounded, deviating significantly from the unrestricted Gaussian limit derived under interior solutions.
Therefore, for valid finite-sample inference, we adopt the conformal-style permutation framework abadie2010synthetic, firpo2009unconditional , which remains exact regardless of the sample size or the boundary nature of the solution, provided the exchangeability assumption holds.
We test the null hypothesis $H_0$ that the intervention has no effect on the distribution of the outcome. Under $H_0$, the treated unit and the donor units are exchangeable regarding the assignment of the intervention. We construct a reference distribution of the treatment effect statistic by iteratively assigning the treatment status to each donor unit $j \in \{2, \dots, J+1\}$.
Let $\hat{\tau}_1$ be the estimated treatment effect for the actual treated unit, defined as the Wasserstein distance between the observed post-treatment distribution and its synthetic counterpart:
Similarly, for each donor $j$, we estimate a placebo synthetic control using the remaining units and compute the placebo effect $\hat{\tau}_j$.
The p-value for the null hypothesis of no effect is calculated as the proportion of units whose estimated effect is at least as large as that of the treated unit:
where $\mathbb{I}(\cdot)$ is the indicator function.
The validity of this p-value relies on the assumption of random assignment of the treatment among the units.
This subsection evaluates the finite-sample stability of the proposed estimator under heavy-tailed data contamination. In empirical economics, particularly within finance and income distribution studies, observed outcomes frequently exhibit fat tails or measurement errors that violate sub-Gaussian assumptions. We test the ability of the Wasserstein-based solver to recover the true structural weights $\bm{\lambda}^{true}$ in a regime where the treated unit's post-treatment distribution follows a Huber contamination model.
We simulate a panel comprising $J=4$ donor units and one treated unit. To align with the asymptotic regime defined in Section 6.2 ("Large Micro-Sample"), we set the number of micro-observations $N=300$ for each unit $j$ and period $t$. The pre-intervention window is set to $T_0=3$ to stress-test the estimator under limited temporal information.
The outcomes are generated via a latent factor structure consistent with Assumption (ref):
where $\alpha_t$ is a common time trend, $F_t \sim \mathcal{N}(0,1)$ represents a time-varying common factor, and $\mu_j \sim \mathcal{U}[0,1]$ is the unit-specific factor loading (corresponding to the latent factors $\mu_i$ in Section 3). The clean counterfactual outcome is constructed as a linear mixture of the donors using the ground-truth weight vector $\bm{\lambda}^{true} = (0.15, 0.25, 0.35, 0.25)^\top$, such that $\lambda_{1} = \sum_{j=2}^{J+1} \lambda_j^{true} \lambda_j$.
To induce contamination, we model the observed treated sample $\{Y_{t,i}^{Treated}\}$ as a mixture distribution:
where $\epsilon \in \{0.01, \dots, 0.04\}$ represents the contamination rate. The contaminating distribution $P_{outlier}$ is drawn from a high-variance normal distribution $\mathcal{N}(\mu_{out}, 10^2)$, with the location parameter $\mu_{out}$ shifted significantly into the right tail of the clean distribution.
We assess the performance of two estimators across $N_{SIM}=100$ Monte Carlo replications. The first is the Benchmark (CDF-$L_2$), serving as a proxy for standard distributional synthetic controls (e.g., Gunsilius, 2023), which minimizes the integrated squared error between empirical CDFs. The second is our Proposed (WGAN-GP) estimator, which minimizes the Wasserstein-1 distance via the Kantorovich-Rubinstein dual, parameterized by a 1-Lipschitz neural network.
Figure (ref) illustrates the geometry of the contamination for period $t=2$. The outliers are strictly localized in the right tail. Under the Wasserstein metric, this shift results in a linear increase in transport cost, whereas under Euclidean metrics, it induces a global distortion of the cumulative probability mass.
We evaluate estimator performance using the Average Root Mean Squared Error (RMSE) of the estimated weights relative to $\bm{\lambda}^{true}$. Figure (ref) presents the results.
The results reveal a fundamental difference in the robustness of the two metrics. The Benchmark $CDF-L_2$ estimator exhibits an immediate breakdown, with the RMSE stabilizing at a high level ($\approx 0.39$) even under minimal contamination ($1\%$). This fragility is intrinsic to the geometry of the $L_2$ norm on function spaces. An outlier in the far tail shifts the empirical CDF ($F_N(y)$) by a constant amount $1/N$ over the entire domain between the bulk and the outlier. The squared error integral accumulates this deviation over the potentially infinite gap, rendering the loss function highly sensitive to the location of extreme values.
In contrast, the proposed WGAN estimator maintains a remarkably stable RMSE ($\approx 0.06$). By constraining the critic function $f$ to be 1-Lipschitz , the Wasserstein objective ensures that the gradient contribution of any data point is bounded. Consequently, the cost of fitting distant outliers grows only linearly with their distance, effectively "capping" their influence. The optimal transport plan thus prioritizes matching the structural core of the distribution rather than distorting the weights to chase tail events.
This subsection evaluates estimator performance under support mismatch, a scenario where the support of the target distribution is disjoint from the convex hull of the donor supports. This violation of the classical overlap assumption renders the synthetic control problem theoretically ill-posed. In such regimes, the objective is no longer to achieve perfect recovery , but to test whether the estimator degrades predictability or signals uncertainty through stable, maximum-entropy weights.
We modify the latent factor structure to induce a controllable geometric gap. Let $Z \sim \mathcal{N}(0,1)$ be a latent generator. The treated unit's factors are generated via a bounded transformation $U_t^{Treated} = 0.5 \tanh(Z)$, confining the target support to the open interval $\mathcal{S}_{1} = (-0.5, 0.5)$.
The donor factors are generated by projecting the same latent $Z$ onto a disjoint domain controlled by a gap parameter $\gamma \in [0, 0.9]$:
The effective support of the donor pool is $\mathcal{S}_{donors} \approx \mathbb{R} \setminus (-\gamma, \gamma)$. As $\gamma$ increases, the distance between the closure of the treated support and the donor support, $d(\bar{\mathcal{S}}_1, \bar{\mathcal{S}}_{donors})$, strictly increases. At $\gamma=0.9$, the distributions are mutually singular with a significant gap.
To handle the potential non-uniqueness of weights in disjoint regimes, we employ the Entropy-Regularized WGAN (WGAN-E) estimator. The regularization term serves as a strict convexifier, selecting the maximum entropy solution among minimizing sequences:
We set $\eta=0.01$ to prioritize transport fidelity while ensuring numerical stability when the transport gradient becomes flat. The benchmark remains the standard CDF-$L_2$ estimator.
We analyze the variance of the estimated weights as a function of the mismatch degree $\gamma$. Figure (ref) plots the average variance of $\hat{\bm{\lambda}}$ ($N_{SIM}=100$) on a logarithmic scale.
The results demonstrate a topological divergence. The Benchmark $L_2$ estimator exhibits instability as the support gap widens, with variance increasing by orders of magnitude at $\gamma=0.9$. This failure is geometric: when $\text{supp}(P_{synth}) \cap \text{supp}(P_{target}) = \emptyset$, the $L_2$ distance between densities (or CDFs) becomes locally constant with respect to infinitesimal weight changes. Consequently, $\nabla_{\bm{\lambda}} \mathcal{L}_{L_2} \approx \mathbf{0}$, causing the numerical optimizer to oscillate between corner solutions driven by machine precision noise.
In contrast, the WGAN-E estimator remains stable. To elucidate the mechanism, Figure (ref) plots the mean estimated weights $E[\hat{\lambda}_j]$.
Figure (ref) confirms that the stability is driven by the regularized transport geometry. Under disjoint supports, the Wasserstein distance strictly penalizes the geometric distance between probability masses. As $\gamma$ increases, the transport cost from any donor to the target becomes approximately equal (dominated by the gap $\gamma$). In this regime, the data term $W_1(\cdot)$ provides no discriminatory power, and the objective function is dominated by the entropy term. Consequently, the estimator converges to the barycenter of the simplex (uniform weights), properly signaling that the data is uninformative. The $L_2$ estimator, lacking both informative gradients and regularization, produces spurious sparsity ($\lambda_j \in \{0,1\}$), creating a false sense of identification.
The third simulation investigates a fundamental theoretical limitation of standard distributional synthetic controls: model misspecification arising from multimodal outcomes. We construct a scenario where the target distribution is bimodal, representing a mixture of sub-populations, while the donor pool consists exclusively of unimodal distributions. This design exposes the structural inability of quantile-averaging methods to recover complex densities that do not satisfy location-scale preservation constraints.
The target distribution $P_t^{Treated}$ is generated as an equiprobable mixture of two distinct Poisson measures :
This specification yields a discrete Probability Mass Function with distinct modes centered at 5 and 20. The donor pool consists of $J=4$ unimodal Poisson distributions $P_{j,t} = \mathcal{P}(\lambda_j)$, where the intensity parameters are given by $\bm{\lambda}_{donors} = \{2, 12, 18, 25\}$.
We benchmark the proposed WGAN-E solver against the standard $W_2$-Quantile estimator. The benchmark constructs the synthetic control by minimizing the $L_2$ distance between the quantile functions of the treated unit and a convex combination of donor quantiles :
This formulation relies on the implicit assumption that the quantile function of the mixture equals the mixture of the quantile functions. However, the mapping from a probability measure to its quantile function is non-linear with respect to convex combinations. Formally, for a weight vector $\bm{\lambda}$, the following inequality generally holds :
Consequently, the quantile averaging approach is theoretically misspecified for multimodal targets constructed from unimodal donors. In contrast, the WGAN-E estimator optimizes weights directly over the space of probability measures, treating the synthetic control as a statistical mixture $\sum \lambda_j P_j$ rather than a mixture of inverses.
Figure (ref) compares the capacity of both estimators to reconstruct the target Probability Mass Function.
The results demonstrate a structural divergence between the two methods. The $W_2$-Quantile solver displayed in the middle panel fails to replicate the bimodal structure. By averaging the quantiles of unimodal donors, the estimator produces a single, smoothed unimodal distribution centered between the two true modes. This outcome is a statistical artifact that represents neither sub-population, confirming that quantile-based synthesis requires the restrictive assumption that the counterfactual belongs to the same location-scale family as the donors .
Conversely, the proposed WGAN-E solver shown in the bottom panel successfully recovers the bimodal topology of the target. By optimizing the Wasserstein distance over the convex hull of donor measures, the estimator correctly identifies a weight vector that mixes donors with low intensity parameters around 5 and high intensity parameters around 20 to reconstruct the two distinct clusters. This result establishes that the measure-based approach remains consistent in scenarios involving complex or polarized distributions where quantile aggregation fails.
Finally, we demonstrate a decisive theoretical and computational advantage of the WGAN-based framework: its natural extension to multivariate outcome distributions. Standard distributional synthetic control methods relying on the inverse probability transform encounter fundamental definition challenges in dimensions $d \ge 2$. This limitation arises because Euclidean space lacks a canonical total ordering, rendering the definition of a multivariate quantile function theoretically ambiguous and computationally intractable . Consequently, existing quantile-based estimators cannot be straightforwardly applied to joint distributions of multiple outcomes, such as inflation and unemployment, or inequality and growth.
In contrast, the Wasserstein distance requires only a metric space and scales to high-dimensional distributions without structural modification. We validate this property through a bivariate proof-of-concept where the outcome lies in $\mathbb{R}^2$ .
We design a scenario where the donor pool consists of $J=4$ bivariate Gaussian distributions $P_j = \mathcal{N}(\bm{\mu}_j, \bm{\Sigma})$. The mean vectors $\bm{\mu}_j$ are positioned symmetrically at the four corners of a square grid, specifically $\{-3, 3\} \times \{-3, 3\}$. All units share a covariance matrix $\bm{\Sigma}$ with a correlation coefficient $\rho=0.5$. The target distribution $P_t^{Treated}$ is constructed as a known mixture of these four donor distributions using the ground-truth weight vector $\bm{\lambda}^{true} = (0.15, 0.25, 0.35, 0.25)^\top$ .
To solve this problem, we adapt the WGAN-E architecture by modifying the critic network to accept two-dimensional inputs $D: \mathbb{R}^2 \to \mathbb{R}$. The generator minimizes the $W_1$ distance as defined in Equation (ref). We note that the $W_2$-Quantile benchmark is excluded from this simulation as the estimator is undefined for multivariate data .
Figure (ref) provides a qualitative validation of the synthesis. The side-by-side comparison reveals a compelling match. The Ground Truth panel displays the multimodal structure of the target characterized by four distinct clusters with varying probability masses. The WGAN-E synthesis accurately replicates this topology. The solver correctly identifies the location of all four modes and assigns the correct relative density to each cluster, reflecting the underlying weights $\bm{\lambda}^{true}$.
This result serves as a proof-of-concept that the adversarial framework effectively circumvents the dimensionality curse of quantile-based methods. By leveraging neural networks to approximate the Kantorovich potential, our method decouples the computational complexity from the geometric complexity of the outcome space. This enables robust distributional synthesis in multivariate settings where traditional sorting-based methods are inapplicable .
To demonstrate the finite-sample robustness of the proposed WGAN-DSC estimator, we apply it to a canonical stress test scenario: the Kansas tax experiment of 2012. While originating from fiscal policy, this episode serves as an ideal testing ground for financial econometrics because it exhibits the stylized facts often encountered in asset pricing and risk management: a massive, exogenous structural shock, heavy-tailed distributional dynamics, and a potential regime switch where the post-treatment support deviates from the historical convex hull.
The 2012 reform represents a distinct regime switch intended to alter the incentive structure for capital accumulation. From an econometric perspective, this intervention tests whether the estimator can distinguish between a shift in the first moment and a fundamental deformation of the distribution's geometry. While theoretical proponents anticipated a first-order stochastic dominance shift, the empirical question is whether the shock primarily induced idiosyncratic volatility and tail risk without a corresponding mean shift. Standard difference-in-differences estimators, which focus on conditional means, fail to capture the geometry of such shifts for the same reason they limit the analysis of asymmetric financial returns.
We construct empirical measures of log real household income using micro-level data from the IPUMS CPS March Supplement (2005--2016). Crucially, this application tests the estimator's ability to handle support mismatch. Similar to liquidity gaps in high-frequency trading, the post-treatment distribution may drift into regions sparsely covered by the pre-treatment donor pool. As established in Section (ref), while Euclidean metrics would face vanishing gradients in such disjoint regimes, the Wasserstein metric leverages the underlying metric space to provide a meaningful signal of distance regardless of overlap.
Figure (ref) verifies the geometric identification condition. The income distribution of Kansas is enveloped by the donor pool, indicating that a convex combination of donors can theoretically reconstruct the target distribution. Unlike standard methods that match only mean moments, our approach requires the donor pool to span the support of the treated unit's probability measure to ensure valid transport.
We apply the WGAN-E estimator to construct a synthetic counterfactual by minimizing the regularized Wasserstein-1 distance between the empirical measure of Kansas and a weighted barycenter of donor states for the pre-treatment years 2005--2011.
Figure (ref) illustrates the goodness-of-fit for 2010, ensuring the validity of the learned weights. The WGAN-based synthetic control reproduces the observed distribution with remarkable precision, capturing not just the location but also the higher-order moments (skewness and kurtosis). This validates that the estimated donor weights $\bm{\lambda}^*$ successfully recover the latent factor structure driving the regional economy prior to the shock.
We estimate the counterfactual distributions for the post-treatment years 2013--2015 to quantify the distributional impact. The results summarized in Figure (ref) reveal a distinct structural decoupling that mean-variance models would miss.
In 2013, the first full year of the new regime, we observe an immediate distributional shock with the Wasserstein distance spiking to $0.088$. Visually, the synthetic counterfactual exhibits a probability mass concentration distinct from the realized outcome. By 2014 and 2015, while the magnitude of the Wasserstein distance slightly attenuates, the structural gap persists.
The shaded distributional gap highlights that the policy did not simply shift the distribution uniformly; rather, the actual distribution lags behind the counterfactual in specific quantiles. This pattern suggests a heterogeneous response analogous to an idiosyncratic volatility shock in asset markets. The successful recovery of the counterfactual distribution in this setting suggests the WGAN-DSC is well-suited for financial applications requiring tail risk attribution and illiquid asset valuation, where capturing the morphological shift of the stochastic process is paramount.
This paper proposes a robust framework for Distributional Synthetic Controls (DSC) grounded in the geometry of Optimal Transport (OT). While the extension of the synthetic control method from mean outcomes to entire probability distributions represents a significant leap in causal inference, existing approaches relying on $L_2$-quantile minimization remain vulnerable to the geometric limitations of Euclidean metrics. We demonstrate that such methods implicitly assume overlapping supports and unimodal structures, rendering them fragile in the presence of outliers, support mismatch, or complex multimodal outcomes.
Our methodological contribution is to reformulate the synthetic control problem as the minimization of the Wasserstein-1 distance between probability measures. By leveraging the Kantorovich-Rubinstein duality within a Generative Adversarial Network (WGAN) architecture, we provide a feasible computational strategy that scales to continuous and high-dimensional outcomes.
The theoretical and simulation results presented herein establish three decisive advantages of this approach. First, the Wasserstein metric provides valid gradient signals even when the supports of the treated and donor distributions are disjoint, resolving the support mismatch failure that plagues standard density-matching estimators. Second, by penalizing transport cost rather than vertical discrepancies, our estimator exhibits superior robustness to heavy-tailed contamination. Third, identifying the synthetic control as a linear mixture of measures—rather than a mixture of quantile functions—allows for the accurate recovery of complex, multimodal counterfactuals where traditional quantile averaging fails.
Practically, this framework expands the scope of the synthetic control method to a broader class of empirical questions. Researchers can now rigorously evaluate the impact of interventions on inequality, polarization, and multidimensional welfare metrics without being constrained by the rigid distributional assumptions of prior methods. As causal inference increasingly moves towards characterizing heterogeneity beyond the mean, the integration of OT principles offers a principled and robust path forward.
\addappheadtotoc