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.
58,102 characters · 12 sections · 42 citation commands
Asymptotic Properties of the Distributional Synthetic Controls
Key words: Distributional synthetic control, Quantile functions, Asymptotic optimality
Causal inference is a pivotal undertaking in social science research, with the synthetic control (SC) method, proposed by abadie2003economic and abadie2010synthetic, serving as a fundamental tool for assessing the causal effects of policies and interventions in settings with observational data. However, the original method of synthetic controls predominantly focuses on point estimates of causal effects on aggregate units, neglecting the heterogeneity present in the distributional characteristics. In many cases, researchers and policy makers want to identify causal effects of policy changes on a treated unit at an aggregate level while having access to data at a finer granularity. A classic example is the evaluation of the effects of minimum wage policies. In this context, the intervention occurs at the state level, yet researchers have access to individual-level data within a state (e.g., card1994minimum, card1994minimum; neumark2000, neumark2000; dube2019, dube2019). These additional data enable the estimation of heterogeneous treatment effects, shedding light on the varied causal impacts of the policy change across different segments of the population within a state. Considering the distributional characteristics, in a seminal paper, gunsilius2023distributional proposes the distributional synthetic control (DSC) estimator. The idea of DSC method is to reconstruct the quantile function associated with the treated unit through a weighted average of quantile functions of the control units, and this weighted average is employed to construct the counterfactual quantile function of the treated unit had it not received treatment.
Compared to the classical SC method, DSC method possesses a notable advantage. It is capable of providing estimates for the effects at different quantiles, enabling researchers to comprehensively understand the impact of interventions. With this distributional information, researchers can estimate the quantile treatment effect (QTE), which offers certain advantages over the average treatment effect (ATE). While ATE often provides a limited perspective on the impact of a treatment, QTE can reveal more comprehensive insights. It is common for a treatment to leave the mean of the outcome distribution unchanged while affecting its dispersion or altering its shape. This granularity is particularly valuable when the treatment effects are not uniform across the population. For instance, a policy intervention might have a larger impact on lower-income individuals compared to higher-income ones. ATE would only show the overall average effect, potentially overlooking the impacts on the lower-income group. In contrast, QTE would reveal these differences by examining the effects at specific points in the distribution, such as the median or the lower and upper quartiles. Therefore, both academic research and practical applications place great importance on understanding the treatment's impact on the entire distribution of outcomes. As stated by tang2020some, rather than focusing on the ATE, applied economists and policymakers are increasingly interested in the distributional treatment effect or QTE. In summary, the QTE function serves as a powerful tool for summarizing the causal effect of a treatment or policy on the marginal distribution of the outcome variable of interest.
Although gunsilius2023distributional has demonstrated the notable performance of the DSC method and explored its identification, there remain other important properties that warrant further investigation. Therefore, in this paper, we provide the asymptotic properties of the DSC estimator. First, we establish the DSC estimator's asymptotic optimality, in the sense that it achieves the lowest possible squared prediction error among all potential treatment effect estimators from averaging quantiles of control units, as the number of draws $M$ goes to infinity. Second, we show that the DSC weight converges to a limiting weight that minimizes the averaged 2-Wasserstein distance of post-treatment periods. Additionally, we quantify the rate of this convergence. We find that an enhanced fit before and after the treatment both facilitate the convergence of the DSC weight. Moreover, a larger number of control units is linked to a slower convergence rate. We also show that increasing $M$ tightens the bound through the term $M^{-1/4}J$. Finally, we provide a data-driven diagnostic by estimating $\xi_t$ from pre-treatment periods. Additionally, the asymptotic property of the DSC estimator, established in this paper, does not rely on the model structure. In other words, it does not need to assume the DGP of the potential outcomes, our asymptotic property holds in a model-free setup. Thus, our work includes the factor model used in many studies as a special outcome model. In the synthetic control literature, zhang2022asymptotic and chen2023synthetic study the large sample properties of SC estimators.
The rest of the paper is organized as follows. Section 2 introduces the DSC estimator and describes the implementation of the DSC method. Section 3 establishes asymptotic properties of the DSC estimator. Section 4 discusses the assumptions required for asymptotic properties. Section 5 reports the results of Monte Carlo experiments. Section 6 draws some conclusions and briefly points to a natural extension of the DSC method, where mixtures of quantile functions are replaced by mixtures of distribution functions; the details of this extension are presented in Appendix F. Technical proofs are given in the Appendix.
In order to facilitate a better understanding of the reader, we will provide a detailed introduction of DSC method below. The methodology and writing style we used in this paper are consistent with gunsilius2023distributional, and the setup and notation closely resemble the classical synthetic controls approach (abadie2003economic, abadie2003economic; abadie2010synthetic, abadie2010synthetic; abadie2021using, abadie2021using).
We possess data pertaining to a set of $J + 1$ units, with the first unit ($j = 1$) designated as the treated unit and the subsequent units ($j = 2, \ldots, J + 1$) designated as the potential control units. The observations span $T$ time periods, and $T_{0}~(T_{0} < T)$ represents the last time period before the treatment in unit $j = 1$. Define $\mathcal{T}_0 = \{ 1, \ldots, T_0 \}$ as the pre-intervention or pre-treatment periods and $\mathcal{T}_1 = \{ T_0+1, \ldots, T_0+T_1 \}$ as the post-intervention or post-treatment periods. All vectors are marked in bold in this paper.
Before delving into the DSC method, we will briefly introduce the classical setting in the literature of SC method. The classical setting focuses on an aggregated outcome, denoted as $Y_{j t}$, observed for each unit $j = 1, \ldots, J+1$ across the time periods $t = 1, \ldots, T$. Potential outcomes are denoted by $Y_{j t, I}$ when unit $j$ is treated at time $t$ and by $Y_{j t, N}$ when no treatment is applied. The standard assumption in this setting is that the intervention has no effect on the outcome before the treatment period, ensuring $Y_{j t, N}=Y_{j t, I}$ for all units $j$ and all pre-treatment periods $t \in \mathcal{T}_0$.
The treatment effect, $\alpha_{j t} = Y_{j t, I}-Y_{j t, N}$, for unit $j$ at time $t$ is defined, allowing the observable outcome to be expressed in terms of counterfactual notation as $Y_{j t} = Y_{j t, N}+\alpha_{j t} D_{j t}$, where $D_{j t}=1$ if $j=1$ and $t\in \mathcal{T}_1$, 0 otherwise. In the classical setting, the interest is to estimate the treatment effect $\alpha_{1t}$ of the treated unit for $t\in \mathcal{T}_1$, that is, $\alpha_{1t} = Y_{1 t, I}-Y_{1 t, N}=Y_{1 t}-Y_{1 t, N}$ for $t\in \mathcal{T}_1$. Thus, the crucial quantity to estimate is $Y_{1 t, N}$ $(t \in \mathcal{T}_1)$, representing the outcome of the treated unit had it not received the treatment in the post-treatment periods.
The distributional setting in gunsilius2023distributional is similar to the classical setting, but with the quantile function $F_{Y_{j t}}^{-1}$ of $Y_{j t}(q)$ as the quantity of interest. The quantile function is formally defined as $$F^{-1}(q):= \inf _{y \in \mathbb{R}}\{F(y) \geq q\}, \quad q \in(0,1),$$ where $F(y)$ is the corresponding cumulative distribution function.
Analogous to the classical setting, the quantiles of potential outcomes are denoted by $F_{Y_{j t, I}}^{-1}(q)$ when unit $j$ is treated at time $t$ and by $F_{Y_{j t, N}}^{-1}(q)$ when no treatment is applied. We define $\alpha_{1 t, q} = F_{Y_{1 t, I}}^{-1}(q) - F_{Y_{1 t, N}}^{-1}(q)$ as the treatment effect of the treated unit for each quantile $q~(q\in (0,1))$ in the DSC setting, corresponding to $\alpha_{1 t}$ at the classical SCM setting, and $\widehat{\alpha}_{1 t, q}(\mathbf{w}) = F_{Y_{1 t, I}}^{-1}(q) - \widehat{F}_{Y_{1 t, N}}^{-1}(q)$ is the estimator of $\alpha_{1 t, q}$, where $\mathbf{w} = (w_2, \ldots, w_{J+1})^{\top}$ is the weight vector belong to the set $$\mathcal{H} = \left\{\mathbf{w}=\left(w_2, \ldots, w_{J+1}\right)^{\top} \in[0,1]^J \mid \sum_{j=2}^{J+1} w_j = 1\right\}.$$ The purpose of DSC method is to estimate the counterfactual quantile function $F_{Y_{1 t, N}}^{-1}(q)$ of the treated unit had it not received treatment. The DSC method constructs the estimation of $F_{Y_{1 t, N}}^{-1}(q)$ through an optimally weighted average of the control quantile functions $F_{Y_{j t}}^{-1}(q), ~j=2, \ldots, J+1$, that is, $$ \widehat{F}_{Y_{1 t, N}}^{-1}(q) = \sum_{j=2}^{J+1} w_j F_{Y_{j t}}^{-1}(q) \quad \text { for all } q \in(0,1). $$
The question moves to how DSC method determines $\mathbf{w}$. This question will be addressed in the subsequent exposition. The estimator $\widehat{\mathbf{w}}$ of $\mathbf{w}$ is a weighted average of the weights $\widehat{\mathbf{w}}_{t}$ over all pre-treatment periods, where $\widehat{\mathbf{w}} = \left(\widehat{w}_{2}, \ldots, \widehat{w}_{(J+1)} \right)^{T}$. Thus, to derive the weight $\widehat{\mathbf{w}}$, we need to first obtain weights $\widehat{\mathbf{w}}_t$ at each time period $t$ ($t \in \mathcal{T}_{0}$), where $\widehat{\mathbf{w}}_t = \left(\widehat{w}_{2t}, \ldots, \widehat{w}_{(J+1)t} \right)^{T}$. For each time period $t \in \mathcal{T}_{0}$, the weights $\widehat{\mathbf{w}}_t \in \mathcal{H}$ are determined to ensure that the weighted average of quantile functions for the control units closely approximates that of the treated unit. To quantify the accuracy of approximation mathematically, gunsilius2023distributional choose the 2-Wasserstein distance as it can simplify the task of determining the weights $\widehat{\mathbf{w}}_t$ into a straightforward regression problem. Following gunsilius2023distributional, the 2-Wasserstein distance, denoted $W_{2}(P_{1}, P_{2})$, between two probability measures $P_{1}$ and $P_{2}$ with finite second moments is defined as $$ W_2\left(P_1, P_2\right)=\left(\int_0^1\left|F_1^{-1}(q)-F_2^{-1}(q)\right|^2 d q\right)^{1 / 2}, $$ where $F_1^{-1}(q)$ and $F_2^{-1}(q)$ are the quantile functions corresponding to $P_1$ and $P_2$, respectively.
Consequently, the method involves determining the weights to minimize the distance between $\sum_{j=2}^{J+1} \widehat{w}_{jt} F_{Y_{j t}}^{-1}(q)$ and the target $F_{Y_{1 t}}^{-1}(q)$ in the 2-Wasserstein space. For each $t \in \mathcal{T}_0$, gunsilius2023distributional determine the weights by solving
The optimization ((ref)) is a convex problem for the weights $\mathbf{w}_{t}$, guaranteeing a unique solution. In practice, the integral can be approximated by randomly sampling a considerable number $M$ of draws $\left\{V_m\right\}_{m=1}^M$ from the uniform distribution on the unit interval, i.e., $V_m \sim$ $U[0,1]$ and solving
It is worth noting that we allow sequence $\{V_m,~m = 1, \ldots, M\}$ to be dependent, aligning more closely with actual data characteristics, while gunsilius2023distributional necessitates their independence.
When the quantile functions $F_{Y_{j t}}^{-1}$ are known, it becomes possible to construct an artificial sample $\widetilde{Y}_{j t m} = F_{Y_{j t}}^{-1}(V_m)$ indexed by $m$ with the number of draws $M$, where the choice of $M$ in the approximation is determined by the researcher. Analogously, the quantiles of potential outcomes are denoted by $\widetilde{Y}_{j t m, I}$ when unit $j$ is treated at time $t$ and by $\widetilde{Y}_{j t m, N}$ when no treatment is applied. In practice, however, the quantile functions $F_{Y_{j t}}^{-1}(q)$ are unknown and necessitate estimated from available data. The empirical quantile functions $\widehat{F}_{Y_{j t n_{j}}}^{-1}(q)$, based on the samples $\{Y_{l,jt}\}_{l=1}^{n_j}$ for $j=1,\ldots,J+1$ and $t\in\{1,\ldots,T\}$, are used as the estimator of $F_{Y_{j t}}^{-1}(q)$, where the subscript $n_{j}$ in $\widehat{F}_{Y_{j t n_{j}}}^{-1}(q)$ denotes that it is based on $n_{j}$ samples. We view $\{Y_{l,jt}\}_{l=1}^{n_j}$ as a sample from the unit--time outcome distribution $F_{Y_{jt}}$. An approach commonly employed for this estimation is through order statistics: $\widehat{F}_{Y_{j t n_{j}}}^{-1}(q) = Y_{t, n_{j} (k)}$, where $k$ is selected such that $(k-1)/n_{j} < q < k/n_{j}$, $Y_{t, n_{j} (k)}$ represents the order statistics of the sample $\{Y_{l, j t}\}, l = 1, \ldots, n_{j}$, $j = 1, \ldots, J+1$, and subscript $n_{j}(k)$ represents the $k$-th sample of the samples $\{Y_{l,j t}\}_{l=1}^{n_j}$ after sorting $\{Y_{l, j t}\}_{l=1}^{n_j}$ in ascending. Correspondingly, in practice, we let $\widehat{\widetilde{Y}}_{j t m} = \widehat{F}_{Y_{j t n_{j}}}^{-1}\left(V_m\right)$ and choose $M$ such that $M=Cn$, where $C$ is a constant and $n=\min\{n_{1}, n_{2}, \ldots, n_{J+1}\}$.
We define the loss function for each $t \in \mathcal{T}_0$ as $L_{t}(\mathbf{w}_t) = M^{-1} \sum_{m=1}^M |\sum_{j=2}^{J+1} w_{j t} \widehat{\widetilde{Y}}_{j t m} - \widehat{\widetilde{Y}}_{1 t m} |^2$.
One can then write the expression ((ref)) as a linear regression, that is,
where $\widehat{\widetilde{\mathbb{Y}}}_t$ is the $M \times J$-matrix with entry $\widehat{\widetilde{Y}}_{(j+1) t m}$ at position $(m, j)$, $\widehat{\widetilde{\boldsymbol{Y}}}_{1 t}$ is the vector of elements $\widehat{\widetilde{Y}}_{1 t m}$ for $m=1, \ldots, M$, and $\|\cdot\|_2$ denotes the Euclidean norm on $\mathbb{R}^M$.
Subsequently, the DSC weight $\widehat{\mathbf{w}}$ can be calculated as a weighted average of the weights $\widehat{\mathbf{w}}_{t}$ over all pre-treatment periods, that is,
$$\widehat{\mathbf{w}} = \sum_{t \in \mathcal{T}_0} \lambda_{t} \widehat{\mathbf{w}}_{t} \quad \text{for } \lambda_{t}\geq 0 \text{ and } \sum_{t \in \mathcal{T}_0} \lambda_{t} = 1.$$ Regarding the choice of weights $\lambda_{t}$, viable options are provided by arkhangelsky2021synthetic, which are also applicable in this case. At every time point $t \in \mathcal{T}_{1}$ within the post-treatment period, the counterfactual quantile function for the treated unit had it not received the treatment is calculated by $\widehat{F}_{Y_{1t n_1, N}}^{-1} = \sum_{j=2}^{J+1} \widehat{w}_{j} \widehat{F}_{Y_{j t n_{j}}}^{-1}$.
In summary, the algorithm for DSC method is shown in Algorithm 1.
In this section, we will list some assumptions and present our theoretical results. Our first result is the asymptotic optimality of the DSC estimator in the sense that it achieves the lowest possible averaged 2-Wasserstein distance of post-treatment periods among all possible averaging estimators over control units, when the $M$ goes to infinity. Additionally, we establish the convergence of DSC weights to the infeasible optimal weights that minimize the averaged 2-Wasserstein distance of post-treatment periods. Unless specified otherwise, all limiting properties hold as $M \rightarrow \infty$.
To facilitate the theoretical analysis, we define the corresponding risk function for each $t \in \mathcal{T}_0$ as $R_{t}(\mathbf{w}_t) = \mathbb{E}_{V_m} [L_{t}(\mathbf{w}_t) \mid \mathcal D]$, where $\mathcal D \equiv \{Y_{l,jt}: l=1,\ldots,n_j;\ j=1,\ldots,J+1;\ t=1,\ldots,T\}$ denotes the collection of observed samples. Note that randomness in our setup arises from two sources: the draws $\{V_m\}_{m=1}^M$ used to approximate the integral, and the observed sample $\mathcal D$ used to construct the empirical quantile functions. To avoid ambiguity, we formulate risks conditionally on $\mathcal D$ and take expectations only with respect to $V_m$. For notational convenience, throughout the paper, we use the shorthand $\mathbb E[\cdot]\equiv \mathbb E_{V_m}[\cdot\mid \mathcal D]$, i.e., $\mathbb E[\cdot]$ denotes expectation with respect to the $V_m$ conditional on $\mathcal D$. Then $R_t(\mathbf w_t)$ can be written equivalently as \[ R_{t}(\mathbf{w}_t) = \frac{1}{M} \sum_{m=1}^M \mathbb{E} \left|\sum_{j=2}^{J+1} w_{j t} \widehat{\widetilde{Y}}_{jtm} - \widehat{\widetilde{Y}}_{1tm}\right|^2. \]
To evaluate the performance of the DSC estimator, we consider the average of the 2-Wasserstein distance at post-treatment period for some weight $\mathbf{w} \in \mathcal{H}$, defined as $$\bar{R}_{T_1}(\mathbf{w}) = \frac{1}{T_1} \sum_{t \in \mathcal{T}_1} \int_0^1 \left( \sum_{j=2}^{J+1} w_{j} F_{Y_{j t}}^{-1}(q) - F_{Y_{1 t, N}}^{-1}(q) \right)^2 dq.$$
Define $\xi_{t} = \inf _{\mathbf{w_t} \in \mathcal{H}} R_{t}(\mathbf{w}_{t})$, and $\bar{\xi}_{T_1}=\inf _{\mathbf{w} \in \mathcal{H}} \bar{R}_{T_1}(\mathbf{w})$. To show the asymptotic optimality and convergence, the following assumptions are used. All explanations of these assumptions are given in the next section.
Denote $e_{t, m, \widehat{\widetilde{Y}}_N}^{(i)} = \widehat{\widetilde{Y}}_{itm,N} - \widehat{\widetilde{Y}}_{1tm,N}$ for $i \in\{2, \ldots, J+1\}$, $m \in \{1, 2, \ldots, M\}$ and $t \in \mathcal{T}_0 \cup \mathcal{T}_1$.
Theorem (ref) establishes the asymptotic optimality of the DSC estimator. Specifically, ((ref)) shows that the DSC weight is asymptotically optimal among all possible weighting combinations in the sense that the averaged 2-Wasserstein distance of the DSC estimator is asymptotically identical to those of the infeasible but best estimator. Moreover, the result in Theorem (ref) can be understood as consisting of two conceptually distinct steps. The first step is a purely statistical result establishing the asymptotic optimality of the proposed estimator with respect to the pre-treatment risk. The second step relies on the identification condition in Assumption 2, which links the pre-treatment and post-treatment risks and allows the post-treatment optimality result in Theorem (ref) to be derived from the pre-treatment optimality result. A more detailed discussion of this decomposition is provided in Appendix A.2.
This optimality statement is closely related to the classical SCM theory of zhang2022asymptotic, who establish an analogous risk-optimality property under the mean-squared prediction error (MSPE) criterion for aggregate-level outcomes. Our result can be viewed as a distributional analogue: we evaluate prediction performance through the squared $2$-Wasserstein distance between outcome distributions, and show that the DSC weights achieve asymptotic optimality. Moreover, in the distributional setting, the loss depends on estimated quantile functions and hence incorporates within-unit sampling variability. Accordingly, our limiting argument is driven by $M\to\infty$ (allowing $T_0$ to be fixed), rather than relying primarily on a long pre-treatment time series. This regime is natural for applications where each unit-time cell contains many micro-level observations, and it clarifies that DSC remains asymptotically optimal even when the number of pre-treatment periods is limited.
We define $\mathbf{\Sigma}_t = M^{-1} \mathbb{E}(\widehat{\widetilde{\mathbb{Y}}}_t^{\top} \widehat{\widetilde{\mathbb{Y}}}_t)$ for $t \in \mathcal{T}_0$, where $\mathbb E[\cdot]$ here denotes expectation over $V_m$ conditional on $\mathcal D$, and use $\lambda_{min}(\cdot)$ and $\lambda_{max}(\cdot)$ to represent the minimum and maximum eigenvalue of a matrix.
The optimal weight vector for a given $T_1$ is defined as the minimizer of the $\bar{R}_{T_1}(\mathbf{w})$, i.e.,
Theorem 2 provides a bound for the Euclidean norm of the difference between the DSC weight and the infeasible optimal weight. If we aim to make the bound on the right-hand side of equation ((ref)) converge to 0, it is not only necessary for $M \rightarrow \infty$, but also for both $\xi_{t}$ for $t \in \mathcal{T}_0$ and $\bar{\xi}_{T_1}$ to tend towards 0; therefore, if these requirements are satisfied, Theorem 2 will establish both the convergence of DSC weight $\widehat{\mathbf{w}}$ to the infeasible optimal weight $\mathbf{w}_{T_1}^{\text {opt}}$ and quantifies the rate of convergence, which depends on $\bar{\xi}$, $\bar{\xi}_{T_1}$ and $J$ as $M\rightarrow \infty$. We discuss the roles of $\bar{\xi}$, $\bar{\xi}_{T_1}$, $J$ and $M$ in turn. First, a faster rate of $\bar{\xi}$ and $\bar{\xi}_{T_1}$ going to zero implies quicker convergence of $\widehat{\mathbf{w}}$. Since $\bar{\xi}$ is the weighted average of $\xi_{t}~(t \in \mathcal{T}_0)$, where $\xi_{t}$ serves as a measure of the fit of the quantiles of the control units to the quantile of the treated unit for each $t \in \mathcal{T}_0$, the theorem establishes a link between good pre-treatment fit and accurate weight estimation. Then, from the term $M^{-1/4} J$, we know that a larger $J$ is linked to a slower convergence rate. Anticipated reductions in estimation accuracy are expected with an increase in the dimension of parameters, given that $J$ corresponds to the number of weight parameters to be estimated. Finally, the role of $M$ is transparent from the term $M^{-1/4}J$: holding $J$ fixed, increasing $M$ decreases this component at rate $M^{-1/4}$.
Theorem 2 is also connected to the weight-convergence theory for classical SCM developed in zhang2022asymptotic. Under an MSPE criterion for aggregate-level outcomes, they provide an explicit upper bound for $\|\widehat{\mathbf w}-\mathbf w^{\mathrm{opt}}_{T_1}\|$ in terms of pre- and post-treatment approximation errors and a term that increases with the donor dimension $J$, with the asymptotic argument primarily driven by the length of the pre-treatment period $T_0$. Our Theorem 2 provides a distributional analogue under the squared $2$-Wasserstein loss: the approximation errors $\bar\xi$ and $\bar\xi_{T_1}$ are defined via quantile-function fitting, and the relevant sample-size dimension is governed by $M$, the within-unit sample size underlying empirical quantiles. Accordingly, Theorem 2 shows that consistent estimation of the DSC weights requires not only a sufficiently large within-unit sample size, but also vanishing approximation errors, so that the treated-unit distribution can be well approximated by convex combinations of donor distributions both in the pre-treatment period and for the post-treatment counterfactual.
This section discusses Assumptions 1-5 and their roles in our theoretical analysis. Since Assumption 2 is stated at a relatively high level, we provide model-based illustrations under concrete and stylized settings to aid interpretation. These conditions are presented as interpretable sufficient conditions and are not additional assumptions required for our main results.
Assumption 1 imposes restrictions on the relative rate of several quantities approaching infinity, i.e., $\xi_{t},~M$ and $J$. It is crucial to highlight that this assumption implies $\xi_{t} \neq 0$, a crucial assumption for establishing the asymptotic optimality of the DSC weight. Intuitively, $\xi_{t} \neq 0$ means that for each $t \in \mathcal{T}_0$, it is impossible to achieve a perfect fit for the pre-treatment quantiles of the treated unit using a linear combination of the quantiles of the control units, and this situation is referred to as an imperfect pre-treatment fit. Assumption 1 can be connected to the classical SC setting (e.g. abadie2010synthetic, abadie2010synthetic), as it mirrors the requirement in the classical SC setting to obtain a unique set of weights. Specifically, if $\xi_t = 0$, then the pre-treatment quantiles of the treated unit lie within the (geodesic) convex hull of the quantiles of the control units. By Carathéodory’s theorem, $\xi_t = 0$ implies that there may exist multiple sets of weights that achieve an exact fit, resulting in the DSC weights failing to be identified. By ruling out this perfect fit scenario, \(\xi_t>0\) ensures the uniqueness of the optimal DSC weight.
Furthermore, $\xi_t$ plays a pivotal role in determining the practical applicability of the theoretical properties, as Assumption 1 -- required for guarantees -- assumes that $\xi_t>0$. The magnitude of $\xi_t$ directly influences the convergence rate of the weight $\widehat{\mathbf{w}}$, making it relevant for both theoretical and empirical considerations. Importantly, in practice, \(\xi_t\) can be estimated from pre-treatment data. This provides a simple, data-driven diagnostic for assessing whether the condition implicit in Assumption 1 is plausible in a given application. We recommend approximating it using the sample analogue
This diagnostic enables researchers to empirically learn about the magnitude and behavior of $\xi_t$. Therefore, reporting $\widehat\xi_t$ can be informative in empirical work as a complementary diagnostic alongside standard pre-treatment fit measures.
We next relate Assumption 1 to the implementation choice of $M$ used in practice, which clarifies that the assumption effectively imposes a requirement on within-unit sample sizes. Assumption 1 is stated in terms of $(\xi_t,M,J)$, whereas our empirical implementation fixes $M$ to be proportional to the minimal within-unit sample size. Specifically, we choose $M=Cn$. Under this choice, Assumption 1 is equivalently rewritten as the sample-size condition $$\xi_t^{-1} n^{-1/2} J^{2} = o(1)~ for~ t \in \mathcal{T}_0.$$ Hence, in the pre-treatment imperfect-fit regime $\xi_t>0$, Assumption 1 imposes a restriction on the sample size. This is natural because, while $M$ controls the discretization used to approximate the integral representation of the Wasserstein loss, the quantile functions $F^{-1}_{Y_{jt}}(q)$ are unknown and must be estimated from $n_j$ individual-level observations. Accordingly, $\widehat F^{-1}_{Y_{jtn_j}}(q)$ is an empirical quantile function whose effective information content is governed by $n_j$. Consequently, choosing $M$ far larger than $n$ mainly re-samples essentially the same set of order statistics more densely and increases computational cost, without delivering commensurate additional information. These considerations motivate letting $M$ grow proportionally with $n$, i.e., $M=Cn$, which aligns the discretization level with the effective sample size while keeping computation manageable.
Assumption 2 constrains the difference between the fits in each pre-treatment period $t$ and the fits in the post-treatment periods. This implies that the primary distinction between the quantiles of the outcomes at each pre-treatment period $t$ and post-treatment period is entirely attributable to the treatment effect. A similar assumption has been discussed in hansen2012jackknife.
In the distributional setting, the quantity we are interested in is the quantile function $F_{Y_{j t}}^{-1}$ of $Y_{j t}$ instead of $Y_{j t}$. To provide sufficient conditions under which Assumption 2 holds, we consider a stylized quantile factor structure similar to that studied in chen2021quantile to generate the potential outcome of the quantile version. This quantile factor representation serves as a quantile-level analogue of the traditional factor models commonly employed in the SCM literature, and is introduced purely as an analytical device rather than as a structural model of distributional dynamics. Specifically, suppose that the potential outcomes $\widetilde{Y}_{itm, N}$ are generated from the following quantile factor model:
where the subscript $m$ has the same meaning as in Section 2, and $v_{m}$ is the observation of $V_{m}$, $\boldsymbol{f}_{t, m} $ signifies an $F_{m} \times 1$ vector of unobserved random common factors, $\boldsymbol{\lambda}_{i,m}$ is an $F_{m} \times 1$ vector of non-random factor loadings.
Since the quantile function $F_{Y_{j t}}^{-1}$ is estimated using the order statistic, we have $\widehat{F}_{Y_{j t n_{j}}}^{-1}(V_{m}) = Y_{t, n_{j} (k)}$. Let $V_{m}^{*} \in (0,1)$ denote the population quantile level such that $F_{Y_{j t}}^{-1}\left(V_m^*\right) = Y_{t, n_{j} (k)}$. That is, $V_{m}^{*}$ is the true quantile level at which the sample quantile $Y_{t, n_{j} (k)}$ lies. Note that $V_{m}^{*}$ depends on $m$, and converges to $V_m$ as the sample size increases. To simplify notation, we write $F_{Y_{j t}}^{-1}\left(V_m^*\right) = \widetilde{Y}_{j t m^{*}}.$ Accordingly, $\widehat{\widetilde{Y}}_{itm, N}$ can be expressed in the form of a factor model: $$\widehat{\widetilde{Y}}_{itm, N} = \boldsymbol{\lambda}_{i, m^{*}}^{\top} \boldsymbol{f}_{t, m^{*}}.$$
Then, Assumption 2 can be derived from more general assumptions as follows. The detailed proof is provided in Appendix C.
\setcounter{assumption}{2}
Assumption 2.1 serves to simplify the proof, and a similar assumption is also employed in ferman2021synthetic and ferman2021properties. Assumption 2.2 implies that for each $t \in \mathcal{T}_0$, the common factors may differ from the average of the common factors across all post-treatment periods. However, this discrepancy gradually diminishes as more samples are used, indicating that the variation of the common factors should not vary significantly after treatment. This assumption means that the treatment effect alone accounts for the majority of the difference between the quantiles of the outcomes at each pre-treatment period $t$ and post-treatment period. Assumption 2.3 requires the uniform boundedness of the factor loadings, and the same assumption is also employed in ferman2021properties. Assumption 2.4 requires that there be no substantial discrepancy between the common factors $\boldsymbol{f}_{t,m}$ and $\boldsymbol{f}_{t, m^{*}}$.
To further illustrate Assumption 2, consider the model introduced in the appendix of gunsilius2023distributional, which is given by
where $\alpha_t$ and $\beta_{t}$ are unknown parameters, and $U_{j,m}$ are independent and identically distributed draws from the unobservable distribution $F_{U_j}$, i.e., $U_{j,m} = F^{-1}_{U_j}(V_m)$. Similarly, we can express $\widehat{\widetilde{Y}}_{itm,N}$ in the same functional form as ((ref)): $\widehat{\widetilde{Y}}_{itm, N} = \alpha_t + \beta_t U_{j,m^{*}}$.
The appendix of gunsilius2023distributional provides an identification motivation suggesting that stable identification of the counterfactual distribution naturally points to an affine (scaled-isometry) relationship between pre- and post-treatment distributional objects (see the discussion and figure therein). Model ((ref)) provides a concrete affine structure at the quantile level. Moreover, ((ref)) exhibits a one-factor structure at the quantile level and is closely aligned with the factor-type illustration in ((ref)).
We now show that Assumption 2 can be derived from more general assumptions as follows. The proof of how these subassumptions jointly imply Assumption 2 is provided in Appendix D.
\noindentAssumption 2.1$^\prime$. {\it There exists a constant $C_3$ such that $\mathbb{E}(U_{j, m^{*}}) < C_3$ for $j \in \{1,\ldots, J+1\}$ and $m\in\{1, 2, \ldots, M\}$.}
\noindentAssumption 2.2$^\prime$. {\it $T_1^{-1} \sum_{t_1\in \mathcal{T}_1} (\beta^2_t - \beta^2_{t_1}) = 0$ for $t \in \mathcal{T}_0$.}
\noindentAssumption 2.3$^\prime$. {\it $M^{-1} \sum_{m=1}^M \mathbb{E} (U_{i,m^*} U_{j,m^*}) - M^{-1} \sum_{m=1}^M \mathbb{E} (U_{i,m} U_{j,m}) = O (n^{-1/2})$ for any $i,j \in \{1,\ldots, J+1\}$.}
Assumption 2.1$^\prime$ imposes a uniform upper bound on the expectation of the latent variables $U_{j,m^*}$. Assumption 2.2$^\prime$ requires that the squared coefficient $\beta^2_t$ in any pre-treatment period $t$ equals the average of $\beta^2_{t_1}$ over all post-treatment periods $t_1$. Notably, this condition is trivially satisfied when the coefficients are time-invariant, i.e., $\beta_t = \beta_{t_1}$ for $t \in \mathcal{T}_0 \cup \mathcal{T}_1$. We emphasize that Assumption 2.2$^\prime$ is introduced to facilitate the derivation of Assumption 2 under the model ((ref)), and thus serves a technical purpose. If the objective is solely to establish the asymptotic optimality of the proposed estimator, this condition can be relaxed. A detailed discussion of this relaxation is provided in Appendix D. Assumption 2.3$^\prime$ requires that there be no substantially discrepancy between the second-moment structure of the latent variables $U_{i,m}$ and $U_{i,m^*}$.
For completeness, an alternative discussion of Assumption 2 under a dynamic panel quantile autoregression framework is provided in Appendix E.
Assumption 3 imposes constraints on the dependency of the potential outcomes $\widehat{\widetilde{Y}}_{itm, N}$ across quantile draws $m$. We illustrate the mildness of Assumption 3 in two scenarios. First, conditional on the observed data $\{Y_{l,jt}\}$, the empirical quantile function $\widehat{F}^{-1}_{Y_{jtn_j}}$ is entirely determined, so $\widehat{\widetilde{Y}}_{itm,N} = \widehat{F}^{-1}_{Y_{itn_i}}(V_m)$ is a deterministic transformation of $V_m$ alone. Thus, the mixing requirement of $\{\widehat{\widetilde{Y}}_{itm,N}\}$ is reduced entirely to a condition of $\{V_m\}$. Assumption 3 holds trivially, when $V_m \overset{iid}{\sim} U[0,1]$ as in gunsilius2023distributional; when $\{V_m\}$ is generated by a standard pseudo-random number generator, Assumption 3 also holds because the dependence decays rapidly by construction. Furthermore, when we consider the randomness of the observed data $\{Y_{l,jt}\}$, note that the randomness of $\widehat{F}_{Y_{itn_i}}$ shrinks when $n_i\to\infty$ and thus the dependence of $\widehat{\widetilde{Y}}_{itm,N} = \widehat{F}^{-1}_{Y_{itn_i}}(V_m)$ across $m$ becomes asymptotically weak (typically of order $1/n_i$). Hence, in large samples the variables $\widehat{\widetilde{Y}}_{itm,N}$ behave nearly as independent draws from ${F}^{-1}_{Y_{it}}(\cdot)$, for any fixed $i\in\{1,\cdots,J+1\}$ and $t\in\mathcal{T}_0\cup\mathcal{T}_1$; therefore, Assumption 3 can be satisfied.
We next illustrate how Assumptions 3 translate into more general assumptions under the two concrete models ((ref)) and ((ref)) introduced above.
\noindentAssumption 3 under the quantile factor model ((ref)). In this case, Assumption 3 can be ensured by imposing a weak-dependence condition on the common factors across $m$: \noindentAssumption 3.1. {\it For any $i \in\{1, \ldots, J+1\}$ and $t \in \mathcal{T}_0 \cup \mathcal{T}_1$, $\{\boldsymbol{f}_{t, m}\}_{m=1}^M$ is either $\alpha$-mixing with the mixing coefficient $\alpha = -r /(r-2)$ or $\phi$-mixing with the mixing coefficient $\phi = -r /(2 r-1)$ for $r \geq 2$. }
Assumption 3.1 ensures that the weak dependence among the common factors $\boldsymbol{f}_{t,m}$ decays sufficiently fast, enabling the application of uniform laws of large numbers and central limit theorems.
\noindentAssumption 3 under the simple linear model ((ref)) of gunsilius2023distributional. In model ((ref)), the sequence $\{U_{j,m}: m=1,\ldots,M\}$ consists of independent and identically distributed draws from the unobservable distribution $F_{U_j}$. Hence, the required weak-dependence condition in Assumption 3 is automatically satisfied in this setting.
Assumption 4 consists of two parts. Assumption 4 (i) implies that the fourth moments of all $\widehat{\widetilde{Y}}_{itm, N}$ can be uniformly bounded. Assumption 4 (ii) concerns the difference between the quantiles of potential outcomes $Y_{jt, N}$ of the treated and control units, ensuring that these variances do not degenerate as $M$ increases. We now illustrate how Assumptions 4 can be ensured under the same two concrete models ((ref)) and ((ref)).
\noindentAssumption 4 under the quantile factor model ((ref)). Define $\mathbf{\Sigma}_{f,t} = \mathbb E(f_{t,m}f_{t,m}^\top)$ for $t\in \mathcal{T}_0 \cup \mathcal{T}_1$ and $m\in\{1, 2, \ldots, M\}$. Assumption 4 can then be further specified as follows:
\noindentAssumption 4.1.
{\it (i) There exists a constant $C_f$ such that $\mathbb{E} \{\|\boldsymbol{f}_{t, m}\|^4\} \leq C_f < \infty$ for $m \in \{1, 2, \ldots, M\}$ and $t \in \mathcal{T}_0 \cup \mathcal{T}_1$.}
{\it (ii) There exist a constant $\kappa_f$ such that $\lambda_{\min}(\mathbf{\Sigma}_{f,t})\geq \kappa_f>0$ for $t \in \mathcal{T}_0 \cup \mathcal{T}_1$ and $m\in\{1, 2, \ldots, M\}$.}
Assumption 4.1 (i) ensures that the fourth moments of the common factors are uniformly bounded, while Assumption 4.1 (ii) requires that the factor covariance matrix $\mathbf{\Sigma}_{f,t}$ is uniformly positive definite over time, preventing the latent factors from becoming degenerate.
\noindentAssumption 4 under the simple linear model ((ref)) of gunsilius2023distributional. Define $e_{U,m}^{(i)} = U_{i,m} - U_{1,m}$ for $i \in \{2,\ldots,J+1\}$ and $m \in \{1,\ldots,M\}$. Assumption 4 can be decomposed as the follows:
\noindentAssumption 4.1$^\prime$.
{\it (i) There exists a constant $C_U$ such that $\mathbb{E} \{U_{j,m}^4\} \leq C_U < \infty$ for $j \in\{1, \ldots, J+1\}$ and $m \in \{1, 2, \ldots, M\}$.}
{\it (ii) There exists a constant $C_4$ such that $\operatorname{var} (M^{-1/2} \sum_{m=1}^M e_{U,m}^{(i)} e_{U,m}^{(j)}) \geq C_4 > 0$ for $M$ sufficiently large and for any $i, j \in\{2, \ldots, J+1\}$.}
Assumption 4.1$^\prime$ (i) ensures that the fourth moments of all $U_{j,m}$ can be uniformly bounded. Assumption 4.1$^\prime$ (ii) concerns the difference between the $U_{j,m}$ of the treated and control units, ensuring that these variances do not degenerate as $M$ increases.
Assumption 5 imposes both lower and upper bounds on the variability of the quantiles of outcomes for each pre-treatment period $t$ of control units. This assumption ensures that the variation among the outcome quantiles of the control units is neither too small nor too large, and it plays a crucial role in establishing the convergence of the DSC weight.
In this section, Monte Carlo simulations are conducted in both model-free and quantile factor model setups to verify Theorems 1 and 2. First, we examine the asymptotic optimality of the DSC estimator and subsequently verify the convergence of DSC weight.
We consider the following simulation design to validate the theoretical results in Section 3. For the time periods $t = 1, \ldots, T$, and $m = 1, 2, \ldots, M$, $\widetilde{Y}_{1 t m}$ are drawn from $ \chi^2(\mu_{1})$, where $\mu_{1}=2$, and $\widetilde{Y}_{j t m}$ are drawn independently from $\mathcal{N}(\mu_{j}, \sigma_{j}^{2})$ for $j = 2, \ldots, J + 1$, where $\mu_{j} \sim U(3, 10)$ and $\sigma_{j} = 3$ ($j$ is odd) or 2.5 ($j$ is even). To allow for dependence across the sampled ranks, we generate $\{V_m\}_{m=1}^M$ in dependent pairs. Specifically, we first draw $M/2$ ranks $V_k^{(1)}\sim U[0,1]$ and construct paired ranks \[ V_k^{(2)}=
\] where $\delta=0.01$, so that the resulting $M$ ranks $\{V_m\}_{m=1}^M$ exhibit dependence across $m$ and remain in $(0,1)$. This construction is one convenient way to introduce within-sample dependence across $m$; other dependence-generating schemes could be used as well. As mentioned above, we set $j=1$ as the treated unit, and $j = 2, \ldots, J + 1$ are the control units. We set $J\in \{20, 50\}$, $M\in \{50, 100, 200, 400\}$, the number of pre-treatment periods $T_0 = 10$ and the number of post-treatment periods $T_1 = 5$. The number of replications is $R = 1000$.
In order to investigate the asymptotic optimality of the DSC estimator, we need to know $\bar{R}_{T_1}(\mathbf{w})$. One can compute $\bar{R}_{T_1}(\mathbf{w})$ as follows: $$\bar{R}_{T_1}(\mathbf{w}) = \frac{1}{T_1} \sum_{t \in \mathcal{T}_1} \int_0^1 \left( \sum_{j=2}^{J+1} w_{jt} F_{Y_{j t}}^{-1}(q) - F_{Y_{1 t}}^{-1}(q) \right)^2 dq.$$
Since any weight $\lambda_{t}~(t\in \mathcal{T}_0)$ satisfying $\lambda_{t}\geq 0$ and $\sum_{t \in \mathcal{T}_0} \lambda_{t} = 1$ can be used, for the sake of simplicity, we use equal weights $\lambda_{t} = 1/ {T_0}$. The weights $\mathbf{w}_t$ in each pre-treatment period $t \in \mathcal{T}_0$ are estimated by equation ((ref)) and the optimal weight vector $\mathbf{w}_{T_1}^{\text {opt}}$ for a given $T_1$ is obtained by $\mathbf{w}_{T_1}^{\text {opt}} = {\arg \min }_{\mathbf{w} \in \mathcal{H}} \bar{R}_{T_1}(\mathbf{w})$.
Figure 1 plots the ratio $\bar{R}_{T_1}(\widehat{\mathbf{w}})/\inf_{\mathrm{w} \in \mathcal{H}} \bar{R}_{T_1}(\mathbf{w})$, under $J=20$ (solid line) and $J=50$ (dashed line), averaged over 1000 replications, as $M$ increases. The curves of the ratio under $J=20$ and $J=50$ both monotonically decrease toward 1 as $M$ increases. This observation indicates that the averaged 2-Wasserstein distance of post-treatment periods of the DSC estimators converges to the lowest possible averaged 2-Wasserstein distance of post-treatment periods as $M$ increases. This result aligns with the asymptotic optimality stated in Theorems 1.
To investigate the convergence of the DSC weight, Figure 2 plots vector norm of the difference between the $\widehat{\mathbf{w}}$ and $\mathbf{w}_{T_1}^{\text {opt}}$ under $J=20$ (solid line) and $J=50$ (dashed line), averaged over 1000 replications, as $M$ increases. We can find that no matter $J=20$ or 50, $\| \widehat{\mathbf{w}} - \mathbf{w}_{T_1}^{\text {opt}} \|$ is monotonically decreasing as $M$ increases, which agrees with the convergence result in Theorem 2. At the same time, comparing the values obtained under the different $J$, we find that $\widehat{\mathbf{w}}$ converges faster when $J = 20$ than $J = 50$, which again agrees with Theorem 2 that the convergence rate slows down when $J$ increases.
We generate the data from the following quantile factor structure: $$\widetilde{Y}_{itm, N} = \lambda_{1, i, m} f_{1, t, m} + \lambda_{2, i, m} f_{2, t, m},$$ where the common factors $f_{s, t, m}$, $s\in\{1,2\}$ are drawn independently from $\mathcal{N}(\mu_{t}, 3^{2})$ for $t = 1, \ldots, T$, and $\mu_{t}\sim \mathcal{N}(0, 1)$, and the factor loadings $\lambda_{s, i, m}$, $s\in\{1,2\}$ are drawn independently from $\mathcal{N}(\mu_{j}, \sigma_{j}^{2})$, where $\mu_{1}=2$, $\mu_{j} \sim U(2, 10)$ for $j = 2, \ldots, J + 1$, and $\sigma_{j} = 2.7$ ($j$ is odd) or 3 ($j$ is even). Similarly, we set $j=1$ as the treated unit, and $j = 2, \ldots, J + 1$ are the control units. We set $J\in \{10, 20\}$, $M\in \{100, 200, 300, 400\}$, the number of pre-treatment periods $T_0 = 10$ and the number of post-treatment periods $T_1 = 5$. The number of replications is $R = 1000$.
As in the previous setup, to investigate the asymptotic optimality and the convergence of the DSC weight, we present the results in Figure 3 and Figure 4. Figure 3 plots the ratio $\bar{R}_{T_1}(\widehat{\mathbf{w}})/\inf _{\mathrm{w} \in \mathcal{H}} \bar{R}_{T_1}(\mathbf{w})$, under $J=10$ (solid line) and $J=20$ (dashed line), averaged over the 1000 replications, as $M$ increases. The curves of the ratio under $J=10$ and $J=20$ both monotonically decrease toward 1 as $M$ increases. This observation indicates that the averaged 2-Wasserstein distance of post-treatment periods of the DSC estimators converges to the lowest possible averaged 2-Wasserstein distance of post-treatment periods as $M$ increases. This result aligns with the asymptotic optimality stated in Theorems 1.
Figure 4 plots vector norm of the difference between the $\widehat{\mathbf{w}}$ and $\mathbf{w}_{T_1}^{\text {opt}}$ under $J=10$ (solid line) and $J=20$ (dashed line), averaged over the 1000 replications, as $M$ increases. We can find that no matter $J=10$ or 20, $\| \widehat{\mathbf{w}} - \mathbf{w}_{T_1}^{\text {opt}} \|$ is monotonically decreasing as $M$ increases, which agrees with the convergence result in Theorem 2. At the same time, comparing the values obtained under the different $J$, we find that $\widehat{\mathbf{w}}$ converges faster when $J = 10$ than $J = 20$, which again agrees with Theorem 2 that the convergence rate is slower when $J$ increases.
In this paper, we investigate the asymptotic properties of the DSC estimator as $M\rightarrow \infty$. We establish the asymptotic optimality of the DSC estimator, in the sense that it achieves the lowest possible averaged 2-Wasserstein distance of post-treatment periods among all possible averaging estimators that are based on an average of quantiles of control units. Furthermore, we show that the DSC weight converges to a limiting weight that minimizes the averaged 2-Wasserstein distance of post-treatment periods. At the same time, we quantify the rate of convergence, providing a better understanding of how the pre- and post-treatment fit, the number of control units and the number of draws $M$ influence the convergence rate. Moreover, we present a natural extension of the DSC method in Appendix F, in which mixtures of quantile functions are replaced by mixtures of distribution functions. This alternative formulation may offer computational or interpretative advantages in certain applications, especially when working directly with estimated distribution functions.