Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
90,010 characters · 22 sections · 18 citation commands
Distributional synthetic controls
The method of synthetic controls, introduced in \citet*{abadie2003economic} and \citet*{abadie2010synthetic}, has become a main tool for estimating causal effects in comparative case studies with aggregate interventions and a limited number of large units. It is designed for settings where some units are subject to a policy intervention and others are not; the respective outcomes of interest are measured in each population before and after the policy intervention, potentially for many periods. The control groups are used to account for unobserved trends in the outcome over time that are unrelated to the effect of the policy intervention. The insight is that an optimally weighted average of the available potential control groups, the synthetic control group, often provides a more appropriate comparison than a single control group alone \citep*{abadie2019using}. The method has been applied in many settings, such as the analysis of the decriminalization of indoor prostitution \citep*{cunningham2017decriminalizing}, the minimum wage debate (\citeauthor*{allegretto2017credible} allegretto2017credible, \citeauthor*{jardim2017minimum} jardim2017minimum, \citeauthor*{neumark2017reply} neumark2017reply), and immigration (\citeauthor*{borjas2017wage} borjas2017wage, \citeauthor*{peri2019labor} peri2019labor); see \citet*{abadie2019using} for further references.
The original method of synthetic controls is designed for settings with aggregate scalar- or vector-valued quantities, in which linear regression approaches are not applicable because of data limitations \citep*{abadie2019using}. Researchers and policy makers are frequently interested in estimating the causal impacts of policy interventions for different quantiles of the population, however. An example is assessing the effects of minimum wage policies on the overall income distribution. An excellent illustration for the relevance of estimating distribution functions is \citet*{ropponen2011reconciling} who applies the changes-in-changes estimator by \citet*{athey2006identification} to estimate the effect of minimum wage changes on employment levels in order to reconcile the classical conflicting results about the minimum wage debate in \citet*{card2000minimum} and \citet*{neumark2000minimum}.
In this article we introduce the distributional synthetic controls method, an extension of the classical method of synthetic controls from aggregate values to quantile functions on the real line. In this case the classical synthetic controls estimator reduces to a constrained linear regression (\citeauthor*{abadie2015comparative} abadie2015comparative, \citeauthor*{doudchenko2016balancing} doudchenko2016balancing) and our method provides a natural extension of this method. The idea of constructing a weighted average of control units that replicates the geometric properties of the target distribution as closely as possible is especially worthwhile in these settings. For instance, in the classical example of a policy change at the state-level \citep*{abadie2010synthetic} our method allows to obtain the entire counterfactual distribution of the treated state after the policy change instead of merely aggregate values. What is more, it allows for causal inference on quantiles of the distribution or any function of the counterfactual distribution such as Lorenz curves \citep*{gastwirth1971general}. This could be important for the research on inequality.
The key concept for our method is a generalization of the linear average of the aggregate outcomes from the classical synthetic control method to appropriate concepts of averages between distribution functions. We argue that the appropriate generalization is to barycenters\footnote{Barycenters, or Fr\'echet means \citep*[chapter 9]{kendall2009shape}, are the natural generalization of linear averages to general metric spaces. Below we introduce these notions formally.} defined on the space of distribution functions with finite second moments, equipped with the $2$-Wasserstein metric\footnote{In different areas the Wasserstein distance is known under the name Mallows distance or earth mover's distance \citep*{levina2001earth}.} \citep*[chapter 6]{villani2008optimal}, or “Wasserstein barycenters” for short \citep*{agueh2011barycenters}. The reason for choosing the Wasserstein barycenter instead of classical linear mixtures is that we want a method that is “geometrically faithful” in replicating the target distribution. This means that it attempts to replicate all geometric properties of the target distribution, such as the support, the moments, the quantiles, and so on. In contrast, standard linear mixture models do not replicate the support or quantiles of the target and are therefore not geometrically faithful as depicted in Figure (ref).
The goal of a synthetic controls method is to obtain a synthetic unit that replicates all quantiles of the target as closely as possible. Wasserstein barycenters provide the framework for this. An example is depicted in Figure (ref), where one wants to find an optimal combination of $4$ Gaussian distributions $\left\{F_{j}\right\}_{j=2,\ldots,5}$ (gray dashed) with the same small variance and different means that replicate as closely as possible the target distribution $F_1$ (black dashed) which is another Gaussian with the same variance but a different mean. The Wasserstein barycenter with equal weights for each donor distribution is geometrically faithful in that it replicates the target distribution perfectly. In contrast, a classical linear mixture approach $\bar{F} = \frac{1}{4}\sum_{j=2}^5 F_j$ is not geometrically faithful: it is a multimodal distribution that does not replicate any part of the target distribution. This is true in any setting, not just for Gaussians: Wasserstein barycenters replicate all moments and quantiles of the target (see e.g. Figure 1 in \citet*{cuturi2014fast} which shows this in a bivariate setting), while standard linear mixtures do not have this property.
The proposed method constitutes a novel concept and is the most natural extension of a linear regression of outcomes to probability measures in Wasserstein space. Note that this differs from the well-known quantile regression \citep*{koenker1978regression}. In classical quantile regression one regresses a specific quantile on covariates. Our method regresses quantile functions on whole quantile functions. The difference is that the explanatory variables in our case are also functional quantities, while the explanatory variables in classical quantile regression are finite-dimensional. In this respect, note that \citet*{bonneel2016wasserstein} introduce a similar version to ours, but for color transport on images, i.e. $2$-dimensional histograms. In this case the simple approach becomes a rather complicated to solve bi-level optimization problem which has many local optima. The authors provide a fixed points approach which only finds a stationary point; a local solution is enough for color transport as the weights are not important. In contrast to color transport, obtaining consistent weights is imperative for a causal inference approach as they are the basis for computing the counterfactual distribution. The proposed method is efficient to implement, works for all type of distributions, and always obtains a unique set of optimal weights.
The proposed method based on Wasserstein barycenters provides a synthetic control unit at the level of interest: the state-level, not the individual level. That is, it matches distributions globally, i.e. it finds an optimal weight for every distribution of the control units that replicates the target distribution as closely as possible. In contrast, linear point-wise approaches like \citet*{chen2020distributional} which is based on classical quantile regression or more na\"ive approaches that (i) decompose distributions into bins and (ii) match these bins between the different distributions are local: they obtain different weights for each point in the distribution and hence obtain a “weight-function”. As a result, these ad-hoc synthetic controls methods are unstable with the choice of points or bins; in particular, they need to make strong assumptions that the optimal weights for each point or bin are “coherent” in the sense that they are the same or at least similar within a given state (e.g. Assumption 1 (ii) in \citet*{chen2020distributional}), something that is rarely satisfied in practice. Such local approaches also cannot provide uniform inference results on the whole counterfactual distribution as they obtain weights that change with each quantile. Finally, these methods implicitly need to assume that each quantile value evolves linearly over time an assumption that is restrictive and rarely met by common stochastic processes. For instance, the transition densities of stochastic differential equations with a nonlinear drift term \citep*{ait1999transition} violate this assumption.
Our method solves all of these issues. First, it provides synthetic control units at the level of interest. Second, it is always geometrically faithful by construction. Third, its implementation is straightforward because it takes the form of a novel linear regression on quantile functions. Fourth, it replicates the whole counterfactual quantile function of the treatment group with a unique set of weights in any setting: for continuous, discrete, or mixed distributions. Fifth, it allows us to provide tests for the counterfactual distribution.
Furthermore, the method is the natural complementary approach to the changes-in-changes estimator \citep*{athey2006identification} in the same sense as the synthetic controls estimator is a complementary approach to the linear difference-in-differences estimator. In particular, we show that the changes-in-changes estimator relies on the exact same mathematical tools as our synthetic controls estimator. The distributional synthetic controls method has three main advantages over the changes-in-changes estimator: (i) it always provides a unique set of weights and counterfactual quantile function, even if the distributions are not continuous, (ii) it does not require the rank-invariance assumption\footnote{Rank-invariance means that the quantiles between the unobservable and the outcome are preserved. For instance, in the classroom setting where the unobservable is ability of the student, rank invariance implies that the student with the highest unobserved ability will receive the highest outcome, the student with the second highest ability will receive the second highest outcome, and so on. The assumption that one can rank individuals based on a one-dimensional criterion and that this ranking stays fixed between distributions can be strong in some settings, see for instance \citet*[p. 131]{dette2016testing}.} on the causal model \citep*{matzkin2003nonparametric} to produce the correct counterfactual distribution, and (iii) it allows for tests of the entire counterfactual quantile function.
The fact that our method is geometrically faithful allows us to perform causal inference on the whole counterfactual quantile function or any function based on it. Rejecting the null-hypothesis of equality of quantile functions is a necessary condition for a causal effect. We also provide large sample results for tests of stochastic dominance and a placebo permutation test in analogy to the classical method \citep*{abadie2010synthetic}: apply the method not just to the actual treatment unit, but also each control unit. If there is a true causal effect, then the change between pre-intervention- and post-intervention period should be largest for the treated unit compared to the control units. The novelty is that one can obtain causal inference on any function of the counterfactual distribution such as averages, quantiles, but also more general concepts such as Lorenz curves \citep*{gastwirth1971general}, Gini-coefficients \citep*{gastwirth1972estimation}, interquartile ranges, or other ROC-curves (\citeauthor*{hsieh1996nonparametric} hsieh1996nonparametric, \citeauthor*{schechtman2019relationship} schechtman2019relationship).
Another contribution of this article is the analysis of a general set of nonlinear causal models for which our method identifies the correct counterfactual distribution. We analyze the nonlinear causal model introduced for the changes-in-changes estimator \citep*[equation (3)]{athey2006identification} through the lens of our synthetic controls approach. The takeaway for applied researchers is that the synthetic controls method, no matter if classical or our distributional extension, works with a causal model that preserves the distance between the outcomes. This follows from the fact that the synthetic controls method depends on an extrapolation step over time: in the pre-period, one finds an optimal combination of the control groups which matches the treatment group as closely as possible. The identifying assumption for causal inference then is that this combination stays optimal in the post-treatment period.
This article complements recently introduced dynamic and nonlinear implementations of the synthetic controls idea, such as \citet*{abadie2017penalized}, \citet*{amjad2018robust}, \citet*{arkhangelsky2019synthetic}, \citet*{athey2018matrix}, \citet*{athey2018design}, \citet*{ben2018new}, \citet*{chernozhukov2018inference}, \citet*{doudchenko2016balancing}, \citet*{imai2019should}, and \citet*{viviano2019synthetic}. The main difference to these approaches is that our method applies the idea of constructing a synthetic control group directly to probability measures and does not consider an aggregate quantity or functionals of the distribution. Also, our method is concerned with estimating the counterfactual distribution for one well-defined treatment intervention and does not deal with the problem of staggered adoption, which some of the other approaches cover. A natural next step for our method is to extend it to staggered adoption settings.
The setup for the proposed method is analogous to the classical synthetic controls approach (\citeauthor*{abadie2003economic} abadie2003economic, \citeauthor*{abadie2010synthetic} abadie2010synthetic, \citeauthor*{abadie2019using} abadie2019using), and the notation is chosen to reflect this fact. Suppose we have data on a set of $J+1$ units, where the first unit $j=1$ is the treated unit and $j=2,\ldots,J+1$ are the potential control units. These units are observed over $T$ time periods, where $T_0<T$ denotes the last time period prior to the treatment intervention in unit $j=1$. In practice, interventions can have an impact prior to their formal implementation, in which case one should define $T_0$ as the last time point before the intervention can have an impact \citep*[p. 494]{abadie2010synthetic}. In accordance with the literature we call $t\leq T_0$ the pre-intervention- or pre-treatment periods and $t> T_0$ the post-intervention- or post-treatment periods.
The classical method of synthetic controls focuses on an aggregated outcome $Y_{jt}$ that is observed for each unit $j=1,\ldots, J$ over the time periods $t=1,\ldots, T$. We denote by $Y_{jt,N}$ the outcome of group $j$ that would be observed at time $t$ in the absence of the intervention; analogously, we denote by $Y_{jt,I}$ the outcome of group $j$ that would be observed at time $t$ if the unit was exposed to the treatment at time $t> T_0$. The standard assumption in this setting is that the intervention has no effect on the outcome before the implementation period, so that we have $Y_{jt,N}=Y_{jt,I}$ for all units $j$ and all pre-intervention periods $t<T_0$.
The key quantity to estimate is $Y_{1t,N}$, the outcome of the treatment unit had it not received the treatment in the post-intervention periods. Based on this, one defines the effect $\alpha_{jt} = Y_{jt,I}-Y_{jt,N}$ of the intervention for unit $j$ at time $t$, so that one can write the observable outcome in terms of the counterfactual notation as \[Y_{jt} = Y_{jt,N} + \alpha_{jt} D_{jt},\] where \[D_{jt} =
\]
The goal in the classical setting is to estimate the treatment effect on the treated group in the post-treatment period, i.e. \[\alpha_{1t} = Y_{1t,I} - Y_{1t,N} = Y_{1t} - Y_{1t,N}\qquad\text{for $t>T_0$,}\] which requires a model for the unobservable $Y_{1t,N}$. \citet*{abadie2010synthetic} introduce a linear factor model \[Y_{jt,N} = \delta_t + \theta_tX_{j}+ w_t\mu_j+\varepsilon_{jt},\] where $\delta_t$ is a univariate factor, $w_t$ is a vector containing factors whose loadings are captured in $\mu_j$, and where $X_j$ are observed covariates of the respective units. The error terms $\varepsilon_{jt}$ are zero mean transitory shocks. The idea for the synthetic controls method is that $Y_{jt,N} = Y_{jt}$ for $t>T_0$ and $j=2, \ldots, J+1$, so that the treatment effect on the treatment group can be obtained by a weighted average \[\hat{\alpha}_{1t} = Y_{1t} - \sum_{j=2}^{J+1}\lambda^*_jY_{jt},\] where $\{\lambda^*_{j}\}_{j=2,\ldots, J}$ is an optimal set of weights.
The classical synthetic controls estimator in this setting then proceeds in two stages.
The idea for the distributional setting considered in this article is analogous to the classical setting, but with $F_{Y_{jt}}$, the probability law of $Y_{jt}$, as the quantity of interest. The idea is identical: we want to approximate the counterfactual distribution $F_{Y_{1t},N}$ by an optimally weighted average of the control distributions $\{F_{Y_{jt}}\}_{j=2,\ldots,J+1}$ that replicates the geometric properties of the target distribution as closely as possible. We choose the Wasserstein barycenter for this task. In particular, we have already argued in the introduction that the classical linear mixture distribution is not a good candidate, as it is not “geometrically faithful”.
Our idea is to instead focus on the associated quantile functions instead of the distributions and to take linear averages of those: \[\sum_{j=2}^{J+1}\lambda_j F^{-1}_{Y_{jt}}(q)\qquad\text{for all $q\in[0,1]$},\] where the quantile function is defined in the usual way as \[F^{-1}(q)\coloneqq\inf_{x\in\mathbb{R}}\left\{F(x)\geq q\right\}.\] This corresponds to the $2$-Wasserstein barycenter in the space of quantile functions and provides the natural extension of the classical synthetic controls estimator. We show this formally in section (ref). Note that the weights $\lambda_j$ do not depend on the quantile, which is the main difference to existing approaches like \citet*{chen2020distributional} or more na\"ive point-wise approaches: the proposed method provides a synthetic control unit at the level of interest, the state-level, not the individual level. This comes with all the advantages mentioned in the introduction.
We now introduce the causal model for the outcome of the groups without treatment, analogous to the linear factor model in the classical setting. This model will be the same causal model as used in the changes-in-changes estimator \citep*{athey2006identification}. In the setting of a few repeated cross-sections, we assume that unobserved heterogeneity $U_j$ for each unit $j$ does not change over time, which leads to the nonseparable model \citep*[equation (3)]{athey2006identification}
Model (ref) is nonseparable in the unobservable $U_j$ and the time period. Here, $U_{j}$ denotes the unobservable characteristics of unit $j$. $h$ is some measure-preserving function which is often called a “production function”. Equation (ref) captures the idea that the outcome of a unit $j$ with $U_j=u$ will be the same in a given time period. The number of time periods can be as low as $2$, with one pre-treatment period and one post-treatment period. In contrast, the classical factor model of the synthetic controls estimator requires an infinite number of pre-treatment periods in order to identify the factors in their model. Model (ref) encompasses the classical linear difference-in-differences model by setting $U_{ij} = \alpha +\gamma\cdot\mathds{1}\{j=1\} +\varepsilon_i$, where $\mathds{1}\{\cdot\}$ is the indicator function taking the value $1$ if the statement in braces is true and $0$ otherwise.
Our setting can also allow for a more general model in which the unobservables are allowed to change over time:
Model (ref) is probably appropriate in settings with short time horizons, i.e. in repeated cross-sections: in this case it is unlikely that the unobservables will change drastically. Model (ref) is more appropriate in a panel setting and for longer time horizons. Our method, just like the changes-in-changes estimator, works in both settings.
The whole question for identification of the counterfactual distribution in this model is what type of functions $h$ provide the correct counterfactual distribution when applying our estimator. In section (ref) we show that one obtains the correct counterfactual for functions $h$ that preserve the Wasserstein distances between the distribution of the unobservables and the outcomes. This is a novel identification result for synthetic control estimators, as we derive it from a nonseparable and nonlinear model.
Before presenting our method, we briefly need to introduce the necessary mathematical machinery, in particular the $2$-Wasserstein space. The $2$-Wasserstein space is a metric space on the set of all probability distributions defined on a general space such as the Euclidean space. In this article we only consider distributions on the real line. The $2$-Wasserstein distance $W_2^2(P_1,P_2)$ between two probability measures $P_1$ and $P_2$ with corresponding distributions $F_1, F_2$ with finite second moments is defined as \citep*[Theorem 2.18]{villani2003topics}: \[W_2(P_1,P_2) = \left(\int_0^1\left\lvert F^{-1}_1(q) - F^{-1}_2(q)\right\rvert^2 dq\right)^{1/2},\] where $F_1^{-1}$ is the quantile function corresponding to $P_1$ and $F_2^{-1}$ is the quantile function corresponding to $P_2$. It is a well-known fact that the quantile function of the $2$-Wasserstein barycenter is simply the weighted average of the quantile functions of the donor distributions (\citeauthor*{agueh2011barycenters} agueh2011barycenters, \citeauthor*{bigot2017geodesic} bigot2017geodesic), i.e. for the barycenter $\bar{P}$ of distributions $\{P_{Y_{jt}}\}_{j=2,\ldots,J+1}$ with weights $\lambda = (\lambda_2,\ldots,\lambda_{J+1})$
The $2$-Wasserstein space in the univariate setting is hence a linear space on quantile functions.
Based on the above concepts, we now introduce the method of synthetic controls. It proceeds exactly as in the classical setting, except that we now assume that we observe the distribution functions $F_{Y_{jt}}$ for $j=1,\ldots, J+1$ instead of aggregated outcomes.
The method of distributional synthetic controls then proceeds in perfect analogy to the classical method. At every time point $t\leq T_0$ we obtain the optimal weights $\vec{\lambda}^*_{t}\coloneqq\left\{\lambda^*_{2t},\ldots,\lambda_{(J+1)t}\right\}\in\Delta^{J-1}$ by finding the $2$-Wasserstein barycenter of the quantile functions $F_{Y_{jt}}^{-1}$ of the control units which is as close as possible (in the squared $2$-Wasserstein distance) to the quantile function $F_{Y_{1t}}^{-1}$ of the target distribution. Just as in the classical setting, $\Delta^{J-1}$ denotes the $J-1$-unit simplex, which is the set of all $J$-dimensional vectors whose elements are non-negative and sum to unity. But we can also allow for extrapolation beyond the unit simplex as in the classical method by allowing the weights to be negative as long as they sum to unity. This immediately yields a regression approach on the Wasserstein space analogous to the Euclidean setting \citep*{abadie2015comparative}.
Formally, we compute \[\vec{\lambda}^*_t=\operatorname*{\arg\!\min}_{\vec{\lambda}\in \Delta^{J-1}}\int_0^1 \left\lvert \sum_{j=2}^{J+1}\lambda_j F_{Y_{jt}}^{-1}(q) - F_{Y_{1t}}^{-1}(q)\right\rvert^2dq,\] which we can write as
where $\mathbb{F}^{-1}_t(q) = \left(F^{-1}_{Y_{2t}}(q),\ldots,F^{-1}_{Y_{(J+1)t}}(q)\right)^T$ is a vector of quantile functions for every $q\in[0,1]$ and $A^T$ denotes the transpose of the matrix $A$. We could also replace the unit simplex $\Delta^{J-1}$ by the set \[\Sigma^{J-1} = \left\{\vec{\lambda}\in\mathbb{R}^J : \vec{\lambda}^T\vec{1} = 1\right\}\] of all weights that sum to unity.\footnote{$\vec{1} = \{1,1,\ldots, 1\}\in\mathbb{R}^J$ denotes the unit vector and $a^Tb$ denotes the inner product in $\mathbb{R}^J$.} This would allow for extrapolation beyond the unit simplex as in the classical method \citep*{abadie2015comparative}, but does not change the method. In this article we focus on the unit simplex as the classical method.
Note that (ref) is a simple convex optimization problem for the weights $\lambda^*_t$ which has a unique solution. In particular, in practice one can compute the integral by randomly drawing a large number $M$ of samples $U_m$ from the uniform distribution on the unit interval $U_m\sim U[0,1]$ and solving
But notice that by the classical inverse transform identity $Y_{jtm} = F^{-1}_{Y_{jt}}(U_m)$, we can write the last expression as a simple linear regression constrained to the unit simplex, i.e.
where $\mathbb{Y}_t$ is the $m\times J$-matrix with entry $Y_{jtm}$ at position $(m,j)$, $\vec{Y}_{1t}$ is the vector of elements $Y_{1tm}$ for $m=1,\ldots, M$, and $\|\cdot\|_2$ is the simple Euclidean norm on $\mathbb{R}^m$. Since we approximate the integral by simulations from a uniform distribution we can make the approximation error as small as desired by choosing a sufficiently large number of simulation samples. Since the objective function is convex---which implies that the optimum is well-separated---and continuous in $\vec{\lambda}$, and since $\vec{\lambda}$ is defined on a convex set, the optimal weights obtained from the approximation converge to the optimal weights of the integral expression as $M\to\infty$ \citep*{newey1994large}.
Note that this approach works for any type of distributions, absolutely continuous, discrete, or mixed, as long as they have finite second moments---otherwise the Wasserstein distance can be infinite and the problem becomes trivial. We can then compute the optimal weights $\vec{\lambda}^*$ as a weighted average of the weights $\vec{\lambda}^*_{t}$ over all pre-intervention periods, i.e. \[\vec{\lambda}^* = \sum_{t\leq T_0} w_{t}\vec{\lambda}^*_t,\qquad\text{for $w_t\geq 0$ and $\sum_{t\leq T_0} w_t = 1$}.\] \citet*{arkhangelsky2019synthetic} provide potential choices of weights $w_t$ that can also be used in our case. In applications, we also found that equal weighs $w_t= \frac{1}{T_0}$ perform well.
At every time point $t>T_0$ in the post-intervention period, we compute the counterfactual quantile distribution for the treatment unit had it not received the treatment as the optimally weighted barycenter in $2$-Wasserstein space via \[F_{Y_{1t,N}}^{-1} = \sum_{j=2}^{J+1}\lambda_{j}^* F_{Y_{jt}}^{-1}.\]
In summary, the whole procedure is as follows.
Just as in the classical setting we will in general obtain “essentially sparse” weights, i.e. weights which only put significantly more weight on a few control units. The reason is that we work in a linear space of quantile functions, so that standard linear concepts like convex hulls of quantile functions and linear projections onto the convex hulls make sense. Whenever the target distribution lies outside the convex hull of the donor distribution we obtain sparse weights just as in the classical setting \citep*{abadie2019using}. We showcase this in Figure (ref) in section (ref).
The method is easy to implement: run a linear regression to obtain the weights and then take a linear weighted average of the quantile functions. Our method provides the counterfactual quantile function $F_{Y_{1t,N}}^{-1}(q)$ which contains almost all of the information of the distribution function. In particular, if the distribution functions are all continuous and strictly increasing, then one can recover the distribution functions from the quantiles by simply taking the inverse. In more general settings, i.e. when the distributions are discrete or mixed, the quantiles still provide the information for most important outcomes of interest. For instance, one can compute the expectation $E[Y_{1t,N}]$ as \[E[Y_{1t,N}] = \int_0^1 F_{Y_{1t,N}}^{-1}(q) dq.\] Based on this one can compute the average treatment effect on the treated unit in time period $t>T_0$ as \[E[Y_{1t}] - E[Y_{1t,N}] = \int_0^1 F_{Y_{1t}}^{-1}(q) dq - \int_0^1 F_{Y_{1t,N}}^{-1}(q) dq.\] Similarly, one can obtain quantile effects for any quantile by \[F_{Y_{1t}}^{-1}(q) - F_{Y_{1t,N}}^{-1}(q).\] Our method allows for more general statistics too like the counterfactual Lorenz curve \citep*{gastwirth1971general} \[L_{Y_{1t,N}}(q) = \frac{\int_0^q F^{-1}_{Y_{1t,N}}(t) dt}{\int_0^1 F^{-1}_{Y_{1t,N}}(t) dt},\] Gini-coefficients \citep*{gastwirth1972estimation}, or interquartile ranges, which could be interesting for questions about inequality. Our approach can provide point-wise, i.e. for fixed $q$, and uniform causal inference on the Lorenz curve using the placebo permutation test Algorithm (ref) in section 5. This algorithm is the exact analogue of the classical synthetic controls estimator.
The proposed method completes the schematic in Figure (ref).
It establishes two existing connections to the existing methods: (i) it is the natural extension of the classical synthetic controls method, and (ii) it is the natural complementary approach to the changes-in-changes estimator in univariate settings.
The proposed method is the correct extension of the classical method in a formal sense: if we apply our distributional method to probability measures supported on just one point, i.e. Dirac measures of the form \[\delta_y(A) =
,\] then we obtain the same results as the classical method.\footnote{We provide a proof in the appendix.} This follows from the fact that the setting of one aggregate value can be phrased as a Dirac measure in probability spaces by the standard embedding $x\mapsto \delta_x$ of $\mathbb{R}$ into $\mathcal{P}_2(\mathbb{R})$. This means that the classical method of synthetic controls is a special case of the proposed distributional method: our method reduces to the classical method when we are only given aggregate values and not distributions.
Note that this does not imply a relation between our method and the classical method applied to moments of a given distribution, such as averages. The classical synthetic controls estimator works with linear averages, which can be computed via \[E[Y_{1t,N}] = \sum_{j=2}^{J+1}\lambda^*_j E[Y_{jt}].\] The $2$-Wasserstein barycenter is the natural analogues of this for quantile functions: \[F^{-1}_{Y_{1t,N}}(q)\coloneqq \sum_{j=2}^{J+1}\lambda^*_j F^{-1}_{jt}(q).\] This means that obtaining the average for given weights $\lambda^*$ in the post-intervention periods our method also implicitly provides an average of the control distributions. In fact, by definition of the expectation as the integral of the quantile function it holds \[E[Y_{1t,N}] = \int_0^1\bar{F}^{-1}_{Y_{1t,N}}(q)dq\coloneqq \sum_{j=2}^{J+1}\lambda^*_j \int_0^1F^{-1}_{Y_{jt}}(q)dq = \sum_{j=2}^{J+1}\lambda^*_j E[Y_{jt}].\] This follows from the fact that the weighted average of quantile functions tries to replicate all moments of the target distribution as closely as possible.
However, in the pre-intervention periods, the DSC estimator will provide different optimal weights than the classical estimator applied to averages. Intuitively, the reason for this is that the distributional synthetic controls method finds optimal weights that replicate all moments of the target distribution as closely as possible. In contrast, applying the classical method to averages of the distributions will find the optimal weights just for the first moment. This can be seen by writing
Taking squares on both sides of the expression it follows by the fact that $dq$ is a probability measure on $[0,1]$ in combination with Jensen's inequality that \[\left( \sum_{j=2}^{J+1}\lambda_j E(Y_{jt}) - E(Y_{1t})\right)^2 \leq \int_0^1\left(\sum_{j=2}^{J+1}\lambda_j F_{Y_{jt}}^{-1}(q) - F^{-1}_{Y_{1t}}(q)\right)^2 dq,\] where the first expression is the objective function for the classical method when applied to the average of a distribution and the second expression is the objective function for the DSC estimator. This shows, unsurprisingly, that the optimal weights can replicate a single number (the average) at least as well as whole distributions. It also gives an indication that in general the weights will be different for the DSC and the classical method applied to expectations, even in binary settings.
Let us now turn to the connection of our method to the changes-in-changes estimator from \citet*{athey2006identification}. The clearest connection is that the changes-in-changes estimator is also based on the mathematical tools used. In particular, the changes-in-changes estimator relies on the theory of optimal transportation and the Wasserstein space, without explicitly mentioning this connection. Recall the approach of the changes-in-changes estimator: in the case of only one pre- and one post-period, as well as only one control group, they construct the counterfactual distribution $F_{Y_{11,N}}$ as $F_{Y_{10}}(F_{Y_{00}}^{-1}(F_{Y_{01}}))$, where the first index is with respect to the treatment group ($0$ for control and $1$ for treatment) and the second index is with respect to the time period ($0$ is pre and $1$ is post).
Note in this respect that the monotone rearrangement $F_{Y_{00}}^{-1}(F_{Y_{01}})$ is the optimal transport map between $F_{00}$ and $F_{01}$ with respect to the $2$-Wasserstein distance. In particular, when $F_{Y_{01}}$ is absolutely continuous, one can compute the $2$-Wasserstein distance between $F_{Y_{01}}$ and $F_{Y_{00}}$ as \citep*[Remark 2.19 (iv)]{villani2003topics} \[W_2(P_{Y_{01}},P_{Y_{00}}) = \left(\int \left\lvert x - F_{Y_{00}}^{-1}(F_{Y_{01}}(x))\right\rvert^2 dF_{Y_{01}}(x)\right)^{1/2}.\] Hence the changes-in-changes estimator relies on the exact same mathematical theory of optimal transportation as our proposed method.
Intuitively, \citet*{athey2006identification} estimate the change of the control unit over time and apply this change to the treatment unit. In contrast, the method of distributional synthetic controls finds the combination of control units that optimally approximates the treatment unit in the pre-intervention periods and then uses this optimal combination of control units to extrapolate the counterfactual distribution in the post period. In a nutshell: the DSC method estimates between units and extrapolates over time, while the changes-in-changes estimator estimates over time and extrapolates between units.
The proposed method has three main advantages over the changes-in-changes estimator in many settings. First, it provides causal inference for the entire counterfactual distribution or any function of it, such as Lorenz curves, Gini-coefficients, or interquartile ranges. Second, it always provides a unique set of optimal weights $\lambda$ in every case: if the distributions are continuous, discrete, or mixed. In contrast, the changes-in-changes estimator can only provide bounds on average-or quantile effects in the discrete case. Third, it does not require the classical rank-invariance assumption \citep*{matzkin2003nonparametric} that the function $h$ is strictly increasing and continuous.
However, both approaches are complementary: the proposed method does require other relatively strong assumptions on the function $h$ which we point out in the next section. These assumption are at the state-level, not at the individual level, however, which might be more reasonable in practice. Moreover, in contrast to the changes-in-changes estimator, the method works with more than a single control group.
We now introduce assumptions on the causal model $h(t,\cdot)$ for which the obtained counterfactual quantile function $F^{-1}_{Y_{1t, N}}$ via our method coincides with the correct counterfactual quantile function. The restrictions on the causal model become apparent when focusing on the simplest case where we have one pre-intervention period ($t=0$) and one post-intervention period ($t=1$). If we face more than two time-periods, we can adjust the setting by simply assuming that there exists one set $\vec{\lambda}^*$ of optimal weights that hold over all time periods, just as in the classical setting \citep*{abadie2010synthetic}. In this setting, we estimate $\vec{\lambda}^*=\sum_{t=1}^{T_0}w_t\vec{\lambda}^*_t$ for some weights $w_t$, as proposed in \citet*{arkhangelsky2019synthetic}.
The structure of the synthetic controls method is an extrapolation over time: we obtain the optimal weights $\vec{\lambda}^*_0\equiv \vec{\lambda}^*$ in the pre-intervention period which provide the closest match of the barycenter of the control units to the target unit; the fundamental implicit assumption of synthetic controls approach, be it classical or our proposed generalization, is that these obtained weights $\vec{\lambda}^*$ from the pre-intervention period are still the optimal weights in the post-intervention periods. This assumption has consequences for the functional form of $h(t,U_j)$.
Figure (ref) depicts a natural sufficient condition on $h(t,\cdot)$ in model (ref): the relative distances between the unobservables $U_{1}, U_{2}, U_{3}, U_{4}$ and the observables $Y_{10}, Y_{20}, Y_{30}, Y_{40}$ remain constant over time. A function $h(t,\cdot)$ that preserves distances is an isometry.
Isometries are maps which conserve the distance between points, so they are continuous and injective by construction. Under the additional requirement that an isometry is surjective, it becomes invertible. In one-dimensional Euclidean space, the Mazur-Ulam theorem \citep*{mazur1932transformations} shows that isometries are affine maps, i.e. compositions of linear functions and shifts. It is here where we can make a connection to the classical method of synthetic controls. In particular, the standard linear factor model for the classical synthetic controls estimator is an affine model which also contains underlying long-run factors. The restriction to isometric/affine models is not surprising as the classical method of synthetic controls requires these restrictions for the extrapolation step going from $t\leq T_0$---where both treatment and control group are observed---to the setting $t>T_0$. In our distributional setting, we require that $h$ is a scaled isometry in the $2$-Wasserstein space, which we look at in more detail below.
Note in this regard that requiring that $h$ be a scaled isometry is only a sufficient assumption. A scaled isometry keeps all distances fixed up to scaling which implies that it preserves barycenters and in particular the weights $\vec{\lambda}_t^*$. However, for our approach to work we only need maps that preserve the weights. These maps can change the overall shape of the convex hull of the control units in general, as long as the target is adjusted so that the optimal weights are preserved. In general, there is not a simple assumption to describe class of functions that preserves the weights $\vec{\lambda}_t^*$, except for isometries. With this assumption we obtain the following result.
The additional requirement in the panel-data case is that the map that changes the unobservables $U_{jt}$ over time is such that the optimal $\vec{\lambda}^*$ stays the same in the post-treatment periods. This map could also be an isometry, but it need not be.
Since isometries are a natural assumption for $h$, we now review some of their properties. In particular, we introduce the properties of isometries in the $2$-Wasserstein space, i.e. the space of all distribution functions with finite second moments equipped with the $2$-Wasserstein distance. As mentioned, in Euclidean space, scaled isometries are affine functions. It would therefore be a parametric assumption on $h$ if we worked in finite-dimensional Euclidean space. We require $h(t,\cdot)$ to be scaled isometries in the $2$-Wasserstein space, however.
What do (scaled) isometries look like in the $2$-Wasserstein space? This has been answered in \citet*{kloeckner2010geometric}. Intuitively, isometries in the $2$-Wasserstein space are surprisingly general in the Wasserstein space on the real line, especially compared to simple affine functions \citep*[Theorem 1.1]{kloeckner2010geometric}. In particular, while isometries on the Euclidean space preserve shapes as depicted in Figure (ref), this is not the case on the Wasserstein space. Intuitively, these isometries on the Wasserstein space only preserve the center of mass of the respective probability distributions and can even deform their support \citep*[section 5]{kloeckner2010geometric}.
So in general, the admissible “production functions” $h(t,\cdot)$ can be more general than simple affine functions for the method of distributional synthetic controls to still work. Assuming linearity in the production function is hence overly restrictive in the setting on the real line (interestingly, in higher dimensions, the admissible isometries are much more restrictive and are only slightly more general than affine functions, see \citet*{kloeckner2010geometric}). This is another extension from the classical synthetic controls method where one assumes a linear model. In particular, the focus on isometries shows that our distributional method can be applied even in the setting with just one pre-intervention and one post-intervention period, analogous to the changes-in-changes estimator \citep*{athey2006identification}. The reason for this is that the assumption of $h$ being scaled isometries is strong enough to not require several time periods for identification.
In most practical settings the distributions $F_{Y_{jt}}$ are not directly known; instead, only $n\in\mathbb{N}$ random draws $Y_{jti}$, $i=1,\ldots,n$, from each distribution are observed. In this case one needs to first estimate $\left\{F_{Y_{jt}}\right\}_{j=1,\ldots,J+1}$ via standard statistical methods, for instance empirical measures or kernel methods. In this section we provide basic statistical analyses for our method in this section. We provide consistency for our general DSC method and also provide the large sample distribution under regularity assumptions on the distributions. This can be used to obtain (uniform) confidence intervals around the counterfactual quantile function.
The empirical analogue to the optimization (ref) for obtaining the weights $\vec{\lambda}^*_t$ for all $t\leq T_0$ is simply the plug-in version
where $\hat{\mathbb{F}}^{-1}_{Y_{tn}}(q) = \left(\hat{F}^{-1}_{Y_{2tn}}(q),\ldots,\hat{F}^{-1}_{Y_{(J+1)tn}}(q)\right)^T$ is a vector of empirical quantile functions and $\hat{F}^{-1}_{Y_{1tn}}(q)$ is the empirical quantile function of the target.
In the second stage, i.e. for $t>T_0$, one then simply has to compute the linear weighted average of the empirical quantile functions. We denote the data as $i=1,\ldots, m$ to show that the weights $\hat{\vec{\lambda}}_{n}^*$ were computed on different data. This leads to the following estimator of the counterfactual quantile function. \[\hat{F}^{-1}_{Y_{1tm,N}}(q) = \sum_{j=2}^{J+1}\hat{\vec{\lambda}}^*_{jn} \hat{F}^{-1}_{Y_{jtm}}(q),\] where \[\hat{\vec{\lambda}}^*_n=(\hat{\lambda}_{2n}^*,\ldots,\hat{\lambda}_{(J+1)n}^*) = \sum_{t\leq T_0} w_t\hat{\vec{\lambda}}^*_{tn}\] are the optimal weights obtained as a weighted average over time.
Based on this counterfactual distribution one can obtain any function of interest such as averages, quantiles, Lorenz curves, or other ROC curves \citep*{hsieh1996nonparametric}.
Consistency of the optimal weights follows under minimal assumptions. In particular, we only need to assume that the second moments of our estimator exist.
Working with distributions instead of mere aggregate values as in the classical method has the added benefit that we can empirically test if the estimated distributions via our method are significantly different from the target distribution. In particular, we can test if the counterfactual distribution $F^{-1}_{Y_{1t,N}}$ is statistically different in the Wasserstein distance from the observed outcome $F^{-1}_{Y_{1t}}$. Focusing on the Wasserstein distance is natural. It is a quantile-quantile plot test in the $L^2$-norm \citep*{ramdas2017wasserstein} and takes into account all geometric properties of the distributions when comparing them.
In order to obtain uniform results for the whole quantile process, we make classical but strong regularity assumptions on the empirical quantile processes (see for instance \citet*{csorgo1993weighted} or chapter 18 in \citet*{shorack2009empirical}).
Assumption (ref) is not the weakest possible. In fact, the assumption that the density is bounded away on its support can be slightly relaxed, see \citep*[Theorem 2.1]{csorgo1990distributions}; this leads to more technical assumptions, however, which distract from the main results in this section. More generally, \citet*{knight2002limiting} provides large-sample result for quantiles under more general assumptions. In this case the large-sample distribution and the rate of convergence can be different form the Gaussian case. The required assumptions are however rather high-level.
Under this assumption we now provide (i) uniform confidence intervals for the quantile functions, (ii) the asymptotic distribution of a two-sample test of equality between the counterfactual distribution and the observed distribution of the treatment unit, and (iii) the large sample distributions for a test of first- and second order stochastic dominance for the observed distribution and the counterfactual distribution. \\
\noindent1. Large-sample distribution for confidence intervals\\ Under the above assumption we obtain the following result.
Four remarks concerning this result are in order. First, in principle this result allows for testing differences in any continuous functional of the quantile functions, such as averages, quantiles, Lorenz curves, and Gini coefficients. In particular, since Lorenz curves are Hadamard differentiable under weak assumptions \citep*{dobler2019bootstrapping}, the rate of convergence will actually stay the same. Second, note that considering $\lambda^*_j$ as fixed quantities in this setting is fine, as these elements were estimated in earlier periods and their estimators are consistent as shown in Proposition (ref). Third, the assumption on the sample sizes $n_{jt}$ is standard and simply requires that the samples are asymptotically balanced. If all sample sizes are the same then the rate of convergence simplifies to $n$. Fourth, note that the large sample distribution depends on the unknown distribution $F^{-1}_{Y_{1t}}$. This means that one cannot directly apply this result for testing the equality of distributions. Since the large-sample distribution is Gaussian, one can rely on the standard bootstrap \citep*[chapter 8]{wasserman2013all} to obtain the confidence intervals for the quantile functions in this setting. \\
\noindent2. Large-sample distribution for test of equality\\ We provide the large sample distribution for two-sided tests of equality of the observed distribution $F_{Y_{1t}}$ and the counterfactual distribution $F_{Y_{1t,N}}$ at every $t>T_0$
We test the hypothesis in the Wasserstein space, i.e. we use the Wasserstein distance to decide if two distributions are equal or not. The asymptotic distribution under the null hypothesis is the following.
The same remarks as above are valid for this result. In particular, the large-sample distribution also depends on the unknown distribution. It is possible to adjust the result in \citet*[Theorem 1]{ramdas2017wasserstein} which provides an asymptotic distribution that does not depend on the unknown distribution in the limit, but this requires that all $f_{Y_{jt}}\circ F^{-1}_{Y_{1t}}$, $j=1,\ldots, J+1$ contain the support of $f_{Y_{1t}}\circ F^{-1}_{Y_{1t}}$, which can be a strong assumption in our setting with potentially many units. Moreover, in contrast to the asymptotic distribution in Proposition (ref) the large-sample distribution in Proposition (ref) is non-Gaussian.
Therefore, the best way to use Proposition (ref) is to perform a Monte-Carlo approach to estimate the asymptotic distribution via a plug-in estimator: estimate all unobserved elements like $f_{Y_{jt}}(F_{Y_jt}^{-1})$ and $\gamma_{jt}$ by their empirical analogues $\hat{f}_{Y_{jtn_{jt}}}(\hat{F}_{Y_jtn_{jt}}^{-1})$ and $\frac{n_{jt}}{\sum_{j=1}^{J+1} n_{jt}}$, sample $B$ random paths from the Brownian bridge process and approximate the integral via some sampling method. For each sample $b=1,\ldots, B$, this provides one value for the test statistic $S_b$. Then one can obtain an approximate $p$-value for the null-hypothesis by comparing the original test statistic $S^0$ to all generated samples $\frac{1}{B}\sum_{b=1}^B\mathds{1}(S^b\geq S^0)$. This provides a simple approximation to the asymptotic $p$-value of the null hypothesis. \\
\noindent3. Large-sample distribution for tests of stochastic dominance\\ We can use the same idea for testing stochastic dominance instead of equality. A distribution $F$ (weakly) dominates another distribution $G$ in the first order if $F(x)\leq G(x)$ for all $x\in\mathbb{R}$. $F$ dominates $G$ in the second order if \[\int_{-\infty}^x\left[G(t)-F(t)\right]dt\geq0.\] The concept of stochastic dominance can be pushed to higher orders, but we only focus on the first two orders. Following the straightforward (convex) duality relations between the integrated distribution functions and the corresponding integrated quantile functions \citep*[Theorem 3.1]{ogryczak2002dual}, we can phrase both first- and second-order stochastic dominance in terms of quantile functions. A test for first-order dominance of the observed distribution over the counterfactual can be implemented as
This hypothesis has become a standard in testing first-order stochastic dominance \citep*{barrett2003consistent}: if the alternative is true then one rejects first-order stochastic dominance. Similarly, we can obtain a test for second order stochastic dominance based on quantile functions \citep*[Theorem 3.2]{ogryczak2002dual} as follows.
Since the null hypotheses in both tests (ref) and (ref) are uniform in $q$, we use the $L^2$-norm over $q$ as before, i.e. \[\int_0^1 \max\left\{F^{-1}_{Y_{1t,N}}(q) - F^{-1}_{Y_{1t}}(q),0\right\}dq\] for the test (ref) and \[\int_0^1 \max\left\{\int_0^q \left[F^{-1}_{Y_{1t,N}}(s) - F^{-1}_{Y_{1t}}(s)\right]ds,0\right\}dq\] for the test (ref). The use of the maximum function for obtaining the large-sample distribution has been proposed in other approaches for stochastic dominance testing, in particular \citet*{schmid1996testing}. Both test statistics are positive if the respective null-hypothesis is violated on a set of positive measure.
This leads to the following result.
The asymptotic distributions in Proposition (ref) are analogous to the goodness-of-fit test. In particular, they also depend on the unknown distributions and are therefore not directly useful for statistical testing of the hypotheses. However, as above, they imply that subsampling routines can be used to perform these tests under Assumptions (ref) and (ref). The requirements for a $\sqrt{n}$ convergence are strong. It is known that the rate of convergence of the expected Wasserstein distance can be quite slow in cases where the distributions are not smooth enough (\citeauthor*{bigot2018upper} bigot2018upper, \citeauthor*{bobkov2019one} bobkov2019one).
In the case where the respective distributions are not absolutely continuous everywhere (for instance when considering financial returns of cross-sectional data of companies at specific time points, where the overall distribution often has atoms at zero), the above results can only be applied at quantiles $q\in(0,1)$ where Assumption (ref) is satisfied, i.e. where $f_{Y_{jt}}(F^{-1}_{Y_{jt}}(q))$ exists and is positive for all $j$. \\
\noindent4. Large sample distribution for equality in the discrete setting\\ Many settings of interest in applied research involve discrete distributions. We therefore now propose the uniform large-sample result for the two-sample goodness-of-fit test in this case. The result is a slight extension of Theorem 2.6 in \citet*{samworth2004convergence}. Note also that recently \citet*{sommerfeld2018inference} derived the large sample distribution for more general underlying spaces, but they are designed for fixed distributions and do not take into account the fact that the counterfactual distribution is computed itself.
Proposition (ref) shows that the rate of convergence in this setting is slower than in the regular continuous case, as the rate for the squared Wasserstein distance is now $\sqrt{n}$, not $n$. The same idea of the proof can be used for stochastic dominance results under minimal changes. Proposition (ref) shows that the asymptotic distribution in the discrete case is also non-standard. Therefore, one should resort to Monte-Carlo approaches as mentioned before in the continuous setting. This is in line with the findings from \citet*{sommerfeld2018inference} who show in a multivariate setting that the Wasserstein distance is only directionally Hadamard differentiable, which implies that the standard bootstrap does not work.
As they stand, the results from the previous section can only be applied in post-intervention periods, as it requires that the optimal weights $\lambda^*$ are known. Since those are computed on independent data-sets in pre-intervention periods via a consistent estimator as shown in Proposition (ref), it is valid to use the estimator $\hat{\lambda}_{n}^*$ computed over the pre-intervention periods.
Confidence bands under Assumption (ref) for the estimated quantile function in the pre-intervention periods can be constructed based on a sample-splitting procedure: since all draws are independently sampled, we can split all samples $\{Y_{jt}\}_{j=1,\ldots,J+1}$ into two independent parts for each $t\leq T_0$. We use the first part to compute the optimal weights $\lambda^*_t$ in each period $t$ and compute the overall optimal $\lambda^* = \sum_{t\leq T_0}w_t\lambda^*_t$ as usual, but only over the first part of the samples that we used to obtain $\lambda^*_t$. On the second sample we can then apply any of the goodness-of-fit tests from the previous sections, in particular the confidence band approach. The downside of this is that the efficiency in both estimating $\lambda^*$ and constructing the confidence bands is reduced due to the sample split, especially compared to the post-intervention periods, where we can use the whole sample in each period $t$. Goodness-of-fit- and stochastic dominance test can be applied via the same sample-splitting idea.
The goodness-of-fit tests from the previous periods are a necessary condition for a causal effect: if we reject the null, then the distributions in the post-treatment are statistically different. It is not a sufficient condition for a true causal effect. In fact, the synthetic controls estimator might not be able to replicate the target well in the pre-intervention periods, so that a significant difference in the post-intervention periods is not an indication for a true causal effect. One way to gather further evidence for a causal effect is a placebo test, analogous to the classical method \citep*{abadie2010synthetic}: apply our method to each control unit as a placebo permutation test. If our model is correct and there is an actual treatment effect only in the treatment group post-intervention, then our method should estimate the effect for the actual treatment unit to be among the most extreme. The pseudo-algorithm Algorithm (ref) provides inference for the counterfactual distribution in a setting with potentially several time periods.
It follows the same underlying idea as the classical placebo-permutation test in \citet*{abadie2010synthetic}. $p_t$ provides the probability of observing a difference between the observable $F^{-1}_{Y_{1t}}$ and the estimated counterfactual $F^{-1}_{Y_{1t,N}}$ given all permutations of the treatment and control groups and hence provides a “significance level” for the permutation test.
The complication compared to the classical method is that our outcome of interest $F^{-1}_{Y_{1t,N}}$ is a functional quantity. Thus, in contrast to the classical setting we use the distances \[d_{\iota t}\coloneqq \int_0^1\left\lvert F^{-1}_{Y_{\iota t,N}}(q)-F^{-1}_{Y_{\iota t}}(q)\right\rvert^2dq\] to rank the difference between the observed $F_{Y_{\iota t}}$ and the computed $F_{Y_{\iota t,N}}$. These distances are always non-negative. In the classical univariate setting one can obtain a direction of the effect. We also rely on the $2$-Wasserstein distance to measure the distance between the counterfactual and the observed distribution. Overall all of these tests can just give an indication of a causal effect, but are not a sufficient condition. It is therefore important to use the method based on a valid model and reasoning to argue that a causal effect exists.
This section provides simulations for the core of our proposed method: finding the optimal weights $\lambda_{jt}^*$ for replicating the target distribution at a fixed time period and computing the synthetic control group as the barycenter. In particular, we want to show how well the “constrained quantile-on-quantile regression” manages to replicate target distributions. We provide evidence for simulated and real data.
We first provide simulation evidence in a setting where all distributions are absolutely continuous. We do this by defining the control distributions to be mixtures of $3$ Gaussians, with uniformly randomly generated means between $-10$ and $10$ and variances between $0.5$ and $6$. To maximize the chance that a target that cannot be perfectly replicated by the control units, we define the target to be a mixture of $4$ Gaussians with uniformly randomly generated means between between $-10$ and $10$ and variances between $0.5$ and $6$. Figure (ref) shows how close the method of distributional synthetic controls can replicate the target for varying numbers of control units in this example.
It is interesting to note that the method manages to replicate the target rather well with already $30$ control units, and replicates it perfectly with $500$. In fact, much fewer control units are necessary for a perfect replication in this case, but the point is to demonstrate that the method works efficiently with many data points and many control units. In particular, the vector of weights $\vec{\lambda}^*$ is “essentially sparse”: roughly $5\%$ of the $500$ control units have a weight of greater than $10^{-4}$. It is not formally sparse, however, i.e. no weight is set to exactly zero. This should not be surprising, as one needs essentially infinite degrees of freedom to replicate an infinite dimensional object like the target, so that even control units with very small weights provide a better fit. However, this essential sparsity shows that one can compress the information of the approach dramatically, which is useful in applied settings.
Even more striking is Figure (ref). It depicts the same setting, but for binomial distributions instead of mixtures of Gaussian.
Here, the control distributions are very different from the target. However, the ways in which the method is able to replicate the target differ dramatically in each case. In the top-left image the method is quite far from the target while in the top-right image it is able to replicate the target rather well. The reason for this is mathematical. On the right the additional control distribution is placed in such a way that the target lies close to the convex hull of the control distributions. We can define the convex hull as \[co_W = \left\{F^{-1}: \thickspace F^{-1} = \sum_{j=2}^{J+1}\lambda_j F^{-1}_{Y_{jt}},\quad \thickspace\lambda=(\lambda_2,\ldots, \lambda_{J+1})\in\Delta^{J-1}\right\}.\] It is therefore possible to replicate the target well by the given control distributions even if those distributions are not similar to the target. In the left image the target lies far outside the convex hull of the distributions and it is not possible to accurate replicate it. In the setting with $1000$ control units, the replication is even better, but not perfect, especially compared to the continuous setting in Figure (ref). This reflects the fact that the barycenter has a slower rate of convergence in the discrete case (Proposition (ref)) compared to the continuous case (Proposition (ref)).
This replication behavior is exactly analogous to the classical method of synthetic controls and relies on the fact that we require all weights to lie in the unit simplex. This means that we only allow for interpolation in our setting, i.e. we only use the information provided by the control units. One can straightforwardly extend our method to the case where the weights are allowed to be negative as long as they sum to unity. This relaxation then corresponds directly to a linear regression on quantile functions and is perfectly analogous to the regression approach in the classical setting \citep*{abadie2015comparative}. The interpretation of such an extrapolation is an open question and is not as simple to interpret as the interpolation approach.
In this section we provide a real-world example for our method. We use a subset of the data provided in \citet*{dube2019minimum} on minimum wages in the United States. The data consist of all $50$ states and the District of Columbia. The outcome of interest is the distribution of equalized family income from wages and salary, defined as multiples of the federal poverty threshold as in \citet*{dube2019minimum}. We pick the final year, $2013$, in the data for illustration. We designate the state of Michigan as the “treated” unit we want to replicate and all other $49$ states and DC as control units. Figure (ref) captures this.
The control units manage to replicate the target perfectly, except for the extreme value for the largest quantile. Figure (ref) also provides the $99\%$ confidence intervals, but they are so tight that they are hard to distinguish from the estimated line. We estimated this example like in a “pre-intervention” period with a $50-50$ data split. On the first sample we estimate the optimal weights $\hat{\vec{\lambda}}$. On the second sample we use those weights to estimate the barycenter and the confidence intervals via a standard bootstrap procedure following the large-sample distribution result from Proposition (ref). Here we also observe the “essential sparsity phenomenon”: $7$ states received weights larger than $0.0001$. $0.286$ for New Mexico, $0.179$ for Illinois, $0.167$ for Alaska, $0.164$ for Ohio, $0.111$ for Tennessee, $0.084$ for Maryland, and $0.009$ for Massachusetts.
This article introduced the method of distributional synthetic controls. It is a natural generalization of the classical and widely-used method of synthetic controls (\citeauthor*{abadie2003economic} abadie2003economic and \citeauthor*{abadie2010synthetic} abadie2010synthetic) which allows the researcher to replicate whole quantile functions and functional data. The key criterion is to obtain geometrically faithful estimators, i.e. estimators that approximate the global geometric properties of the target distribution, such as the support.
We achieve this by relying on the concept of barycenters in Wasserstein spaces \citep*{agueh2011barycenters}. The idea is to find weights such that the corresponding barycenter of the control distributions for these weights is as close as possible to the target. Just like the classical synthetic controls estimator, it requires the weights to lie in the unit simplex, which corresponds to interpolation. One can straightforwardly extend the method to the extrapolation case, analogously to \citep*{abadie2015comparative}, but the interpretability is not clear for these settings. In the latter case, our approach provides the most natural definition of a linear regression of outcomes on Wasserstein space, which defines a mathematical concept that could be of interest in itself. In this article we focus on the method as a tool for causal inference and show in simulations that it performs extraordinarily well in practice, without the need for any tuning parameters.
Our approach allows us to perform causal inference on many outcomes of interest beyond mere averages or quantiles. For instance, we can perform inference on measures of inequality such as Lorenz curves or Gini-coefficients. Working with distributions instead of aggregate measures allows us to perform goodness-of-fit tests of the hypothesis that the counterfactual is equal to the observed distribution of the treatment. This is a necessary condition for a causal effect, but it is not sufficient. We also provide a placebo permutation approach analogous to the classical synthetic controls method \citep*{abadie2010synthetic}. In the supplementary material we briefly argue that current existing conformal prediction approaches are not applicable in our setting. It could hence be fruitful to construct conformal prediction approaches for whole probability measures. Another interesting extension of our approach is to allow for staggered adoption.