Distributionally Robust Synthetic Control: Ensuring Robustness Against Highly Correlated Controls and Weight Shifts

abstractThe synthetic control method estimates the causal effect by comparing the treated unit’s outcomes to a weighted average of control units that closely match its pre-treatment outcomes, assuming the relationship between treated and control potential outcomes remains stable before and after treatment. However, the estimator may become unreliable when these relationships shift or when control units are highly correlated. To address these challenges, we introduce the {\bf D}istributionally {\bf Ro}bust {\bf S}ynthetic {\bf C}ontrol (DRoSC) method, which accommodates potential shifts in relationships and addresses high correlations among control units. The DRoSC method targets a novel causal estimand defined as the optimizer of a worst-case optimization problem considering all possible weights compatible with the pre-treatment period. When the identification conditions for the classical synthetic control method hold, the DRoSC method targets the same causal effect as the synthetic control; when these conditions are violated, we demonstrate that this new causal estimand is a conservative proxy for the non-identifiable causal effect. We further show that the DRoSC estimator's limiting distribution is non-normal and propose a novel inferential approach. We demonstrate its performance through numerical studies and an analysis of the economic impact of terrorism in the Basque Country.

bibunit\section{Introduction} The synthetic control (SC) method abadie2003economic, abadie2010synthetic, abadie2015comparative is increasingly used in social sciences for its transparency and interpretability. The method estimates the counterfactual outcome by constructing a weighted average of control units that closely matches the treated unit’s pre-treatment trajectory, approximating its potential outcome under control and estimating the causal effect via comparison with the observed post-treatment outcome. The SC framework has inspired a wide range of methodological developments under various structural assumptions, including linear prediction models li2020statistical, chernozhukov2021exact, cattaneo2021prediction, chernozhukov2025debiasing, factor models xu2017generalized, shi2021theory, ben2021augmented, quantile functions gunsilius2023distributional, and matrix completion approaches amjad2018robust, bai2021matrix, athey2021matrix. The SC method is recognized as a key contribution to the policy evaluation literature athey2017state; see abadie2021using for a comprehensive review. Despite its widespread use, the SC method has limitations that can compromise the reliability of constructing a synthetic control unit. In particular, \begin{itemize} • Weight instability: with highly correlated controls, multiple weight configurations may yield similar pre-treatment fits, leading to unstable counterfactual predictions and treatment effect estimates. • Weight shifts: changes in the treated--control relationship between the pre- and post-treatment periods, which may render the synthetic control a biased approximation of the treated unit's counterfactual. \end{itemize} When either challenge arises, the treatment effect is no longer point-identifiable via the SC method. To address this, we define a novel causal estimand through the lens of distributionally robust optimization (DRO). When neither of the aforementioned issues occurs, we show that this estimand coincides with the identifiable treatment effect. In contrast, in the presence of the aforementioned challenges, the proposed estimand provides a conservative lower bound for the treatment effect while preserving its sign. \subsection{Our Results and Contributions} In this paper, we introduce the Distributionally Robust Synthetic Control (DRoSC) estimator, which targets a novel causal estimand–the weight-robust treatment effect–defined as the solution to a DRO problem within the SC framework. This estimand is identifiable even when the true treatment effect is not. Rather than assuming a uniquely identifiable post-treatment weight (as in the standard SC framework), we consider a class of plausible weights arising from weight shifts or highly correlated controls and define the robust causal estimand as the treatment effect that minimizes the worst-case risk over this class, yielding an estimand that remains meaningful even when standard SC assumptions fail. As our main result, Theorem (ref) shows that the weight-robust treatment effect is the optimal value of a degenerate constrained convex optimization problem and, importantly, it is uniquely identified as the most conservative treatment effect over all post-treatment weights. Although meaningful, statistical inference for the weight-robust treatment effect is challenging because the DRoSC estimator may have a non-standard limiting distribution, making conventional asymptotics unreliable. We address this with a perturbation-based method for constructing valid confidence intervals (CIs) that separates regular uncertainty from an irregular component induced by the geometry of the constrained optimization. To quantify the uncertainty arising from the nonregular component, we generate a collection of perturbed optimization problems and show that at least one perturbed problem nearly recovers the population optimization problem. We leverage this observation to quantify the uncertainty from the nonregular component and construct a valid CI even when the DRoSC estimator’s limiting distribution is not standard (e.g., non-normal). We evaluate our proposal across diverse data-generating processes, including highly correlated controls and post-treatment weight shifts. In regimes with non-regular asymptotics, the standard Wald CIs suffer from undercoverage, whereas our perturbation-based CIs attain nominal coverage; see Section (ref). We further illustrate the practical utility of DRoSC with a reanalysis of the Basque Country case study abadie2003economic. To summarize, our main contributions are as follows: \begin{itemize} • We introduce DRoSC as a generalization of the standard SC method. The DRoSC estimator targets a new causal estimand, the weight-robust treatment effect, which remains interpretable even when the treatment effect is not point-identified. • We propose a perturbation-based inference procedure that yields valid CIs even under non-regular limiting distributions {and is of independent interest for non-regular inference in convex optimization with non-unique minimizers.} \end{itemize} \subsection{Other Related Works} In addition to the works discussed above, we review other relevant literature by topic, with additional references on non-regular inference provided in Appendix (ref). {\bf Sensitivity analysis and DRO.} While recent work has incorporated sensitivity analyses into SC to address identification violations, it focuses on different forms of weight misspecification. zeitler2023non analyze bias from distributional shifts in latent causes but do not address highly correlated controls. ferguson2020assessing study sensitivity to model misspecification, allowing post-treatment weights outside the simplex, but do not formally incorporate statistical uncertainty. In contrast, we focus on identification failure driven by weight shifts and highly correlated controls. Moreover, we adopt a DRO-based framework that yields a single, interpretable causal estimand, unlike conventional sensitivity analyses that deliver sets of plausible values manski1990nonparametric. To our knowledge, the current work is the first to explicitly connect DRO and sensitivity analysis in the SC literature, with Theorem (ref) providing geometric intuition: the weight-robust treatment effect is the most conservative effect over all post-treatment weights in the uncertainty class. {\bf Non-unique SC weights.} The most relevant work on highly correlated controls and non-unique synthetic weights is abadie2021penalized, who consider multiple treated units and propose a penalized SC method to promote uniqueness. Their approach penalizes covariate discrepancies but is not applicable when only outcomes are available doudchenko2016balancing,amjad2018robust,chernozhukov2021exact. In contrast, rather than mitigating non-uniqueness through penalization, we address non-uniqueness by introducing an uncertainty class that fully accounts for all possible synthetic weights. {\bf Notation.} For any $v\in\mathbb{R}^d$, $v_j$ denotes its $j$-th element, $\|v\|_{q} = (\sum_{i=1}^d v_i^q)^{1/q}$ for $q\ge 1$, and $\|v\|_{\infty}=\max_{i}|v_i|$; $\mathbf{1}_d$ and $\mathbf{0}_d$ denote $d$-dimensional vectors of ones and zeros. For any $n\times d$ matrix $M$, $M^{\mkern-1.5mu\mathsf{T}}$ is its transpose, $M_{i,j}$ its $(i,j)$ entry, and $\|M\|_{\max}=\max_{i,j}|M_{i,j}|$. Let $\mathbf{I}$ denote the identity matrix; for symmetric $M\in\mathbb{R}^{d\times d}$, $\|M\|_2$ and $\lambda_{\min}(M)$ denotes its largest and smallest eigenvalues. For a sequence $x_n$, $x_n\to x$, $x_n\xrightarrow{p} x$, and $x_n\xrightarrow{d} x$ denote convergence, convergence in probability, and convergence in distribution. For positive sequences $a_n$, $b_n$, $a_n\lesssim b_n$ means $a_n\le C b_n$ for some constant $C>0$, and $a_n\asymp b_n$ if $a_n\lesssim b_n$ and $a_n\gtrsim b_n$. For a random vector $X_n$, $X_n=O_p(1)$ denotes boundedness in probability. i.i.d.\ means independent and identically distributed, and $\mathbb{I}(\cdot)$ is the indicator function. \section{Synthetic Control: Essential Assumptions and Challenges} We review the standard SC setup with $N+1$ units observed over $T$ time periods. Let $T_0\le T-1$. Unit $1$ is untreated for $t=1,\dots,T_0$ (pre-treatment) and treated for $t=T_0+1,\dots,T$ (post-treatment), while units $2,\dots,N+1$ remain untreated throughout. Let $Y_{j,t}^{(0)}$ and $Y_{j,t}^{(1)}$ denote the potential outcomes of unit $j$ at time $t$ under control and treatment, respectively. Accordingly, the observed outcomes satisfy $Y_{1,t}=Y_{1,t}^{(0)}$ for $t\le T_0$, $Y_{1,t}=Y_{1,t}^{(1)}$ for $t>T_0$, and $Y_{j,t}=Y_{j,t}^{(0)}$ for $2\le j\le N+1$ and $1\le t\le T$. We define the average treatment effect on the treated (ATT) at time $t$ shi2021theory as $\tau_t = \mathbb{E}[Y_{1,t}^{(1)} - Y_{1,t}^{(0)}]$, where the expectation is taken over the randomness of $Y_{1,t}^{(1)}$ and $Y_{1,t}^{(0)}$ under a super-population framework imbens2015causal. Thus, we write \begin{align} Y_{1,t}^{(1)}-Y_{1,t}^{(0)} = \tau_t+v_t \quad for\quad t=T_0+1,\ldots,T, \end{align} where $v_t$ is a mean-zero error term. We further define the time-averaged ATT: \begin{align} \bar{\tau} = \frac{1}{T_1} \sum_{t=T_0+1}^{T} \tau_t \quad with\quad T_1=T-T_0. \end{align} Our primary focus is on inference for $\bar{\tau}$ arkhangelsky2021synthetic, liu2024proximal, or its conservative proxy when $\bar{\tau}$ is not identifiable. We allow for non-constant effects $\tau_t$, generalizing the constant-effect assumption in li2020statistical and shi2021theory. Inference for a single $\tau_t$ is more challenging since only a single treated unit is observed at $t$, typically requiring additional assumptions such as a static effect ($v_t=0$) abadie2010synthetic or a parametric model for $\tau_t$ park2025single. Without such assumptions, researchers turn to constructing a prediction interval for $Y_{1,t}^{(1)} - Y_{1,t}^{(0)}$ cattaneo2021prediction. In contrast, focusing on $\bar{\tau}$ enables valid inference without imposing additional assumptions on $\tau_t$ or $v_t$, or resorting to prediction intervals. \subsection{Essential Assumptions for the Synthetic Control Method} We discuss the key assumptions under which the SC method identifies $\bar{\tau}$. Throughout, we write $X_t = (Y_{2,t}, \dots, Y_{N+1,t})^{\mkern-1.5mu\mathsf{T}}$ and consider the following models between unit 1 and controls chernozhukov2021exact, shen2023same: \begin{align} Y_{1,t}^{(0)} = \begin{cases} X_t^{\mkern-1.5mu\mathsf{T}} \beta^{(0)} + u_t^{(0)}\ & \text{for } t=1, \dots, T_0, \\ X_t^{\mkern-1.5mu\mathsf{T}} \beta^{(1)} + u_t^{(1)}\ & \text{for } t=T_0+1, \dots, T, \end{cases} \quad \text{with}\quad \beta^{(0)},\beta^{(1)} \in \Delta^{N}, \end{align} where $\Delta^{N} = \{ \beta : \beta_j \geq 0, \ \mathbf{1}_N^{\mkern-1.5mu\mathsf{T}} \beta = 1 \}$, and $\{u_t^{(0)}\}_{t=1}^{T_0}$ and $\{u_t^{(1)}\}_{t=T_0+1}^{T}$ are sequences of mean-zero error terms satisfying $\mathbb{E}[X_t u_t^{(0)}] =\mathbb{E}[X_t u_t^{(1)}]= \mathbf{0}_N$. We now state the two critical assumptions for SC to identify $\bar{\tau}$ under model (ref). \begin{enumerate} • $\beta^{(0)}$ is the unique minimizer of the following constrained least squares: \begin{align} \beta^{(0)}=\underset{\beta \in \Delta^{N}}{\operatorname*{arg\,min}}\; \frac{1}{T_0}\sum_{t=1}^{T_0}\mathbb{E}\left[\left(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta\right)^2\right]. \end{align} • There is no weight shift before and after the treatment: $\beta^{(0)} = \beta^{(1)}.$ \end{enumerate} Under conditions (E1) and (E2), the SC method identifies ${\beta}^{(0)}$ via (ref) and $\bar{\tau}$ as $T_1^{-1}\sum_{t=T_0+1}^{T}\mathbb{E}[Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta^{(0)}]$, which motivates the following SC estimators of $\beta^{(1)}$ and $\bar{\tau}$: \begin{align} \widehat{\beta}^{\rm SC}=\underset{\beta \in \Delta^{N}}{\operatorname*{arg\,min}}\; \frac{1}{T_0}\sum_{t=1}^{T_0}\left(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta\right)^2\quad\text{and}\quad \widehat{\tau}^{\rm SC} = \frac{1}{T_1}\sum_{t=T_0+1}^{T}\left(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{\rm SC}\right). \end{align} \subsection{Identification Challenges: Non-uniqueness and Weight Shift} We discuss how the identification conditions (E1) and (E2) may fail to hold in practice. First, the minimizer in (ref) may not be unique when control units' pre-treatment outcomes are highly correlated. For example, in the Basque study abadie2003economic, most correlations between selected and all control units are near one (Figure (ref)). Second, the treated-control relationship may not remain stable, as the treatment itself can alter it. \begin{figure}[ht] \caption{Correlation plot from the Basque study. The vertical axis denotes control units selected from SC, and the horizontal axis denotes all control units.} \end{figure} To demonstrate that the SC estimator can be unreliable when (E1) and (E2) fail, we conduct semi-real data analyses using the Basque study abadie2011synth. Detailed implementation steps are provided in Appendix (ref). We first investigate the instability of $\widehat{\beta}^{\rm SC}$ by adding small random noise to the pre-treatment data and analyzing the resulting estimates. Specifically, we add noise with standard error equal to $c$ times that of the corresponding pre-treatment data with $c \in \{0.05, 0.1, 0.15\}$. For each $c$, Figure (ref) reports the proportion of control units selected by SC across 1000 simulated pre-treatment data, displaying only those units that were selected more than 20% of the simulations when $c=0.015$. While Madrid and Baleares remain consistently selected, Rioja’s selection frequency declines as noise increases, whereas similar units—Asturias and Cataluna—are selected more often. This plot illustrates selection uncertainty arising from a violation of (E1): with highly correlated control units, multiple nearly equivalent weight configurations can yield comparable pre-treatment fits. \begin{figure}[ht] \caption{ The proportion of control units being selected by the SC method out of 1000 perturbed data sets. The variable \texttt{c} indicates the noise level applied to the dataset to generate the pre-treatment data.} \end{figure} Next, we examine how weight shifts affect the SC estimator’s performance. To generate semi-real data, we set $\beta^{(0)}=\widehat{\beta}^{\rm SC}$ and construct $\beta^{(1)}$ by shifting weights away from $\beta^{(0)}$ across pairs of similar regions (Baleares-Cataluna and Rioja-Asturias), with shift magnitude controlled by $\kappa\in[0.05,0.4]$, as shown in the left panel of Figure (ref). We generate 1,000 semi-real datasets by adding noise to the control units and computing the treated unit’s outcomes using (ref) and (ref) with $\beta^{(0)}$ and $\beta^{(1)}$. For each dataset, we compute $\widehat{\tau}^{\rm SC}$ in (ref). The right panel of Figure (ref) shows that $\widehat{\tau}^{\rm SC}$ fails to accurately estimate $\bar{\tau}$ as $\kappa$ grows, reflecting bias induced by weight shifts and thus a violation of (E2). \begin{figure} \raisebox{-10px}{ \begin{minipage}[t]{.15\linewidth} \begin{tabular}[b]{c|cc} Region &$\beta^{(0)}$ &$\beta^{(1)}$\\ \hline Madrid & 0.483 & 0.483 \\ Baleares & 0.311 & $(1-\kappa)\cdot0.311$\\ Rioja & 0.206 & $(1-\kappa)\cdot0.206$\\ Cataluna & 0 & $\kappa\cdot0.311$\\ Asturias & 0 & $\kappa\cdot0.206$ \\ \hline \end{tabular} \end{minipage} } \begin{minipage}[t]{.6\linewidth} \end{minipage} \caption{The left table reports the pre- and post-treatment weights used to generate the model (ref). The right panel shows violin plots of $\widehat{\tau}^{\rm SC}$ in (ref) across 1000 simulations for each $\kappa \in \{0.05, 0.1, 0.2,0.3, 0.4\}$; the blue dashed line denotes the true $\bar{\tau}$.} \end{figure} \section{Distributionally Robust Synthetic Control} As discussed in Section (ref), identification of the time-averaged ATT $\bar{\tau}$ fails when key identification conditions are violated. To address this, we introduce a new causal effect via distributionally robust optimization (DRO) in Section (ref). In Sections (ref) and (ref), we establish its identification and compare it with sensitivity analysis, respectively. \subsection{Novel Causal Estimand via Distributionally Robust Optimization} In the following, we introduce a novel causal estimand as a proxy for $\bar{\tau}$ in (ref). To motivate the definition, we express the time-averaged ATT $\bar{\tau}$ as the solution to an optimization problem and then generalize its definition borrowing the strength of DRO. We start with the objective function of the optimization problem and show that $\bar{\tau}$ is the maximizer of the following optimization problem (ref). For any given weight vector $\beta$, we define the reward function $R_{\beta}(\tau)$ associated with the weight $\beta$ and the treatment effect $\tau$ as \begin{equation} R_{\beta}(\tau) \coloneqq \frac{1}{T_1}\sum_{t=T_0+1}^{T}\mathbb{E}\left[ \left(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta\right)^2 - \left(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta-\tau\right)^2\right]. \end{equation} For given $\beta$ and $\tau$, $R_{\beta}(\tau)$ compares the prediction error under a null treatment effect to that under a constant treatment effect $\tau$, representing the improvement in fit when $\tau$ is introduced to quantify the treatment effect in the post-treatment period. Thus, for a given $\beta$, we aim to maximize $R_{\beta}(\tau)$, as larger values indicate greater reduction in prediction error. When the oracle knowledge of $\beta^{(1)}$ is available, we write the time-averaged ATT $\bar{\tau}$ as the solution to the following optimization problem, \begin{equation} \bar{\tau} = \operatorname*{arg\,max}_{\tau \in \mathbb{R}} R_{\beta^{(1)}}(\tau). \end{equation} Note that, for a given $\beta \in \Delta^{N}$, the maximizer of $R_{\beta}(\tau)$ is given by \begin{align} \tau(\beta) = \mu_Y-\mu^{\mkern-1.5mu\mathsf{T}}\beta, \quad \text{with}\quad \mu_Y = \frac{1}{T_1} \sum_{t=T_0+1}^{T} \mathbb{E}[Y_{1,t}]\quad\text{and}\quad \mu = \frac{1}{T_1} \sum_{t=T_0+1}^{T} \mathbb{E}[X_t]. \end{align} It follows from the definition of (ref) and the model (ref) that $\bar{\tau}=\tau(\beta^{(1)}).$ Although we write $\bar{\tau}$ as the solution to the optimization problem in (ref), the identification challenge for $\bar{\tau}$ remains, since identifying $\beta^{(1)}$ still hinges on (E1) and (E2). To address this challenge, instead of attempting to recover the true $\beta^{(1)}$ by imposing conditions such as (E1) and (E2), we introduce a new estimand motivated by distributionally robust optimization (DRO) ben2009robust,duchi2021learning. For $\lambda \geq 0$, we define the following uncertainty class: \begin{align} \Omega(\lambda)\! \coloneqq \left\{\beta \in \Delta^{N}\!: \left\|\frac{1}{T_0}\sum_{t=1}^{T_0}\mathbb{E}\!\left[X_t(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta)\right]\right\|_{\infty}\!\leq \lambda\right\}=\left\{\beta \in \Delta^N\!: \|\gamma - \Sigma\beta\|_{\infty}\!\leq \lambda\right\}, \end{align} where $\Sigma = T_0^{-1}\sum_{t=1}^{T_0}\mathbb{E}[X_tX_t^{\mkern-1.5mu\mathsf{T}}]$ and $\gamma = T_0^{-1}\sum_{t=1}^{T_0}\mathbb{E}[X_tY_{1,t}]$. The uncertainty class $\Omega(\lambda)$ consists of all weight vectors $\beta$ for which the time-averaged covariance between $X_t$ and the residuals $Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta$ over the pre-treatment period is small. This criterion is motivated directly by the logic of the population pre-treatment SC weights $\beta^{(0)}$: in the population pre-treatment SC problem, any optimal $\beta^{(0)}$ satisfies the orthogonality condition: after fitting with $\beta^{(0)}$, the pre-treatment residual $Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\beta^{(0)}$ should not exhibit systematic association with any control series. When control units are highly correlated, $\beta^{(0)}$ may not be uniquely identified; however, this does not invalidate the criterion. Rather, $\Omega(0)$ contains \emph{all} population-optimal pre-treatment weights, i.e., all candidate $\beta^{(0)}$ achieving equally good pre-treatment fit. The parameter $\lambda$ controls the allowable deviation from this pre-treatment criterion. When $\lambda = 0$, the set reduces to weights that exactly satisfy the moment condition (in the population sense), and it contains all candidate $\beta^{(0)}$. As $\lambda$ increases, $\Omega(\lambda)$ admits larger deviations and therefore a wider range of post-treatment weights, with the goal of containing the true post-treatment weight $\beta^{(1)}$ when weight shifts are present. Thus, all weight vectors with deviation level up to $\lambda$ are treated as potential candidates for $\beta^{(1)}$. The choice of $\lambda$ reflects the user’s belief about the extent of weight shift, and in practice one may conduct a sensitivity-style analysis by varying $\lambda$ rosenbaum2002observational. Since $\beta^{(1)}$ can be any element of $\Omega$, we consider all weights in $\Omega$ and define the worst-case reward for a candidate treatment effect $\tau$ over $\Omega$ as $\min_{\beta \in \Omega} R_{\beta}(\tau)$. By taking the minimum over all admissible post-treatment weights, this formulation evaluates $\tau$ under the most adverse (yet plausible) weight configuration in $\Omega$. Similar to (ref), we define the weight-robust treatment effect as the optimizer of this worst-case reward: \begin{align} \tau^*(\Omega) \coloneqq \operatorname*{arg\,max}_{\tau \in \mathbb{R}}\left[\min_{\beta \in \Omega}R_{\beta}(\tau)\right]. \end{align} When (E1) and (E2) hold so $\bar{\tau}$ is identifiable by the SC method, the new estimand $\tau^*(\Omega)$ reduces to $\bar{\tau}$ by using $\Omega$ with $\lambda=0$. Even when $\bar{\tau}$ is not identifiable, $\tau^*(\Omega)$ remains identifiable and serves as a conservative proxy for $\bar{\tau}$ as established in Theorem (ref). Our choice of the reward function $R_{\beta}(\tau)$ in (ref) ensures that $\tau^*(\Omega)$ in (ref) admits a meaningful interpretation and a connection to sensitivity analysis, as established in Section (ref). We can interpret $\tau^*(\Omega)$ through a game-theoretic lens blackwell1979theory: nature adversarially selects the worst-case post-treatment weight from $\Omega$, while the decision-maker chooses $\tau$ to maximize the resulting reward. Thus, $\tau^*(\Omega)$ represents a treatment effect that is robust to adversarially chosen (yet plausible) post-treatment weights. \subsection{Identification and Interpretation} In this subsection, we present the identification theorem of $\tau^*(\Omega)$ defined in (ref). \begin{Theorem} $\tau^*(\Omega)$ defined in (ref) is uniquely identified as \begin{align} \tau^*(\Omega) = \mu_Y-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*(\Omega) \quad \text{where} \quad \beta^*(\Omega) = \operatorname*{arg\,min}_{\beta \in \Omega}\left[\mu_Y-\mu^{\mkern-1.5mu\mathsf{T}}\beta\right]^2. \end{align} \end{Theorem} Theorem (ref) provides a method for explicitly computing $\tau^*(\Omega)$ by first identifying the adversarial weight $\beta^*(\Omega)$ through solving a quadratic program and then computing $\tau^*(\Omega)=\mu_Y-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*(\Omega)$. Intuitively, $\beta^*(\Omega)$ is the post-treatment weight in $\Omega$ that drives the time-averaged ATT closest to zero, making $\tau^*(\Omega)$ the most conservative treatment effect. This identification does not rely on the key assumptions of the SC method. When clear from context, we denote $\tau^*(\Omega)$ and $\beta^*(\Omega)$ as $\tau^*$ and $\beta^*$, respectively. We shall remark that the quadratic program in (ref) is degenerate since $\mu\mu^{\mkern-1.5mu\mathsf{T}}$ is rank one. Consequently, while the optimal value $\tau^*$ is unique, $\beta^*$ may not be unique, complicating estimation, inference, and theory. In particular, this degeneracy leads to a slower convergence rate than the parametric rate; see Theorem (ref). Building on Theorem (ref), we present a theorem that interprets our proposed causal estimand. \begin{Theorem} For $\tau^*$ defined in (ref), we attain the following equivalent expression: \begin{align*} \tau^* = \begin{cases} \underset{\beta\in\Omega}{\min} \ \tau(\beta) & \text{if} \ \tau(\beta) > 0 \ \text{for all} \ \beta\in\Omega, \\ \underset{\beta\in\Omega}{\max} \ \tau(\beta) & \text{if} \ \tau(\beta) < 0 \ \text{for all} \ \beta\in\Omega, \\ 0 & \text{if there exists} \ \beta\in\Omega \ \text{such that} \ \tau(\beta) = 0, \end{cases} \end{align*} where $\tau(\beta)$ is defined in (ref). \end{Theorem} Theorem (ref) establishes that $\tau^*$ is the point in the range of $\{\tau(\beta)\}_{\beta \in \Omega}$ closest to the origin. Intuitively, this means that we consider all possible time-averaged ATTs with weights in $\Omega$ and then the new causal estimand $\tau^*$ represents the most conservative time-averaged ATT. Theorem (ref) implies the following result, relating $\tau^*(\Omega)$ to $\bar{\tau}$ under the condition that $\beta^{(1)}\in \Omega$. \begin{Theorem} If $\beta^{(1)} \in \Omega$, then $\tau^*(\Omega)$ does not have an opposite sign to $\bar{\tau}$ and $|\tau^*(\Omega)|\leq |\bar{\tau}|$. \end{Theorem} Theorem (ref) shows that $\tau^*(\Omega)$ is a conservative proxy for $\bar{\tau}$ whenever $\Omega$ is large enough to contain $\beta^{(1)}$. When (E1) and (E2) hold so $\bar{\tau}$ is identifiable, setting $\lambda=0$ yields $\Omega=\Omega(0)$ and hence $\tau^*(\Omega)=\bar{\tau}$. When identification of $\bar{\tau}$ fails due to a failure of (E1) or (E2), $\tau^*(\Omega)$ remains conservative: its magnitude is bounded above by that of $\bar{\tau}$, and its sign cannot be opposite to that of $\bar{\tau}$. In particular, whenever $\tau^*(\Omega)\neq 0$, $\tau^*(\Omega)$ and $\bar{\tau}$ agree in sign. Figure (ref) illustrates how $\tau^*(\Omega)$ operates in both identifiable and non-identifiable settings. \begin{figure}[ht] \caption{Relationship between $\tau^*$ and $\bar{\tau}$ when $\Omega$ contains $\beta^{(1)}$. The circles labeled E1 and E2 denote the settings where conditions (E1) and (E2) are satisfied, respectively, while the rectangle represents the general case, including violations of these conditions. } \end{figure} \subsection{Comparison to Sensitivity Analysis} Theorems (ref) and (ref) also connect our framework to sensitivity analysis. Using the sensitivity parameter $\lambda$, we specify a family of plausible post-treatment weights $\Omega(\lambda)$ in (ref). If the true post-treatment weight satisfies $\beta^{(1)}\in\Omega(\lambda)$ as required in Theorem (ref), then the population sensitivity interval $[\min_{\beta\in\Omega(\lambda)}\tau(\beta),\ \max_{\beta\in\Omega(\lambda)}\tau(\beta)]$ contains $\bar{\tau}=\tau(\beta^{(1)})$, yielding partial identification of $\bar{\tau}$ manski1990nonparametric. Importantly, $\tau^*(\Omega(\lambda))\neq 0$ implies that the sensitivity interval excludes zero and therefore identifies $\mathrm{sgn}(\bar{\tau})$, aligning with Theorem (ref). While sensitivity analysis emphasizes this set-valued conclusion, our DRO formulation complements it in two ways. First, it defines a principled point estimand in the partial identification regime, $\tau^*(\Omega(\lambda))$, which admits a game-theoretic interpretation and provides an informative summary of the sensitivity interval. Second, it enables direct statistical inference for this robust target: we {shall} develop in Section (ref) a perturbation-based procedure that yields valid CIs for $\tau^*(\Omega(\lambda))$ even when the conventional asymptotics are unreliable. More broadly, the same perturbation idea provides guidance for inference in sensitivity analysis, where inference for the end points of the sensitivity interval may inherit similar non-regularity. \section{Estimation Procedure} We devise a data-dependent estimator of $\beta^*$ and $\tau^*$ based on the identification established in Theorem (ref), beginning with a data-dependent estimator of the uncertainty class $\Omega$ in (ref): \begin{align} \widehat{\Omega}(\lambda) = \left\{\beta \in \Delta^{N}: \|\widehat{\gamma}-\widehat{\Sigma}\beta\|_{\infty} \leq \lambda+ \rho \right\}, \end{align} where $\widehat{\Sigma} = T_0^{-1}\sum_{t=1}^{T_0}X_tX_t^{\mkern-1.5mu\mathsf{T}}$ and $\widehat{\gamma} = T_0^{-1}\sum_{t=1}^{T_0}X_tY_{1,t}$ and $\rho$ is a tuning parameter used to account for the estimation errors in $\widehat{\Sigma}$ and $\widehat{\gamma}$. We discuss selection of the tuning parameter $\rho$ in practice at the end of the section. When unambiguous, we denote $\widehat{\Omega}(\lambda)$ as $\widehat{\Omega}$. Building on the identification in Theorem (ref), we estimate $\beta^{*}$ and $\tau^*$ by \begin{align} \widehat{\beta}(\widehat{\Omega})\coloneqq\operatorname*{arg\,min}_{\beta \in \widehat{\Omega}}\left[\widehat{\mu}_Y-\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\beta\right]^2,\quad \widehat{\tau}(\widehat{\Omega}) = \widehat{\mu}_Y - \widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}(\widehat{\Omega}), \end{align} where we estimate $\Omega$ by $\widehat{\Omega}$ in (ref) and estimate $\mu_Y$ and $\mu$ by the corresponding sample averages $\widehat{\mu}_Y = T_1^{-1}\sum_{t=T_0+1}^{T} Y_{1,t}$ and $\widehat{\mu} = T_1^{-1}\sum_{t=T_0+1}^{T}X_t$. We denote $\widehat{\beta}(\widehat{\Omega})$ and $\widehat{\tau}(\widehat{\Omega})$ as $\widehat{\beta}$ and $\widehat{\tau}$ when there is no confusion. We refer to our procedure as the Distributionally Robust Synthetic Control (DRoSC) method, summarized in Algorithm (ref) in Appendix (ref). {\bf Tuning Parameter Selection.} To estimate $\Omega$, we replace $\gamma$ and $\Sigma$ in (ref) with their empirical counterparts $\widehat{\gamma}$ and $\widehat{\Sigma}$. This plug-in step introduces additional estimation uncertainty, so we introduce a tuning parameter $\rho$ to enlarge the uncertainty set and account for this error. For i.i.d.\ data, our theory suggests a data-dependent choice of $\rho$: \begin{align} \rho &= C\left[\widehat{\sigma}\cdot\max_{2 \leq j \leq N+1} \left(\frac{1}{T_0}\sum_{t=1}^{T_0}Y_{j,t}^2\right)^{\frac{1}{2}}+\lambda\right]\frac{{\log (\max\{T_0,N\})}^{1/2}}{\sqrt{T_0}}, \end{align} where $\widehat{\sigma}^2 = T_0^{-1}\sum_{t=1}^{T_0}(Y_{1,t}-X_t^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{\rm SC})^2$ and $C>0$ is a constant; see Appendix (ref) for the justification. To determine $C$ in (ref), we initialize $C$ as a small constant (e.g., $C=0.01$). In some instances, this choice can make (ref) infeasible when using $\widehat{\Omega}$ in (ref). To address this, we iteratively increase $C$ by a factor of $1.25$ until (ref) admits a feasible solution. This procedure yields the smallest $\rho$ (equivalently, the smallest $C$) for which feasibility is attained. We discuss the choice of $\rho$ for non-i.i.d.\ data in Appendix (ref). \section{DRoSC Inference: Perturbation-based Methods} We turn to the statistical inference for the estimand $\tau^*$. We demonstrate the related inference challenge in Section (ref) and devise a novel perturbation-based inference in Section (ref). \subsection{Inference Challenge: Non-regularity and Instability} The inference challenge arises since the estimator $\widehat{\tau}$ in (ref) may not admit a standard limiting distribution. The estimation error decomposes as $\widehat{\tau}-\tau^* =\widehat{\mu}_Y-\mu_Y-(\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*).$ While $\widehat{\mu}_Y-\mu_Y$ is asymptotically normal by the central limit theorem, we demonstrate in Figure (ref) that $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$ may exhibit a non-regular behavior due to the boundary constraint on $\beta^*$ and the highly correlated controls. In Figure (ref), we plot histograms of $\widehat{\tau}$ and $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}$ based on 500 simulations; see Section (ref) for the data generating details. The left panel exhibits near-normal behavior, whereas the other panels show deviations from normality, leading to undercoverage of normality-based CIs. \begin{figure} \begin{subfigure}{.72\textwidth} \end{subfigure} \begin{subfigure}{.72\textwidth} \end{subfigure} \caption{Histograms of $\widehat{\tau}$ and $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}$ based on 500 simulations from setting (S2) in Section (ref). From left to right, the panels correspond to $\tau^*\approx-0.6$, $\tau^*\approx 1.35$, and $\tau^*\approx 0.05$ with $T_0 = T_1 = 25$. In the top and bottom panels, the red solid line denotes $\tau^*$ and $\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$. The blue dashed lines in both panels denote the sample averages across the 500 simulations.} \end{figure} We explain that the nonregular behavior of $\widehat{\tau}-\tau^*$ in Figure (ref) is driven primarily by that of $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$. First, the simplex (boundary) constraint on $\beta^*$ can make the active set change across samples, which in turn yields a non-standard (typically mixture) limiting distribution for $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$. Such boundary constraints lead to non-regular asymptotics self1987asymptotic,andrews1999estimation,drton2009likelihood. As a consequence, standard procedures that rely on asymptotic normality or bootstrap/sub-sampling may fail andrews2000inconsistency,wasserman2020universal, guo2023causal, xie2024repro, kuchibhotla2024hulc. Second, with highly correlated controls, $\Omega$ can become nearly flat along certain directions, so small estimation errors in $\gamma$ or $\Sigma$ may lead to a substantial discrepancy between $\widehat{\Omega}$ and $\Omega$. When $\beta^*$ lies near the boundary of $\Omega$, such discrepancies can cause $\widehat{\beta}$ to cross between regions where different constraints are active, thereby inducing both non-regularity and instability. \subsection{Perturbation-based Inference} We now propose a novel perturbation method to address the inferential challenge highlighted in Section (ref). We begin with the intuition, recalling the identification result from Theorem (ref): \begin{align} \beta^*=\operatorname*{arg\,min}_{\beta \in \Omega}\left(\mu_Y-\mu^{\mkern-1.5mu\mathsf{T}}\beta\right)^2\quad\text{with}\quad\Omega = \left\{\beta \in \Delta^N:\|\gamma-\Sigma\beta\|_{\infty}\leq \lambda\right\}. \end{align} The data-driven estimator $\widehat{\beta}$ presented in (ref) is to replace $\{\Sigma,\gamma,\mu_Y,\mu\}$ in (ref) with their sample analogs $\{\widehat{\Sigma},\widehat{\gamma},\widehat{\mu}_Y,\widehat{\mu}\}$. Our main idea is to add perturbation to the sample-based optimization problem in (ref) and create a collection of perturbed optimization problems, with the hope that one of these perturbed optimization problems almost recovers (ref). We generate $M$ perturbed quantities $\{\widehat{\Sigma}^{[m]},\widehat{\gamma}^{[m]},\widehat{\mu}^{[m]}_Y,\widehat{\mu}^{[m]}\}_{m=1}^M$ by adding perturbations to the sample estimator $\{\widehat{\Sigma},\widehat{\gamma},\widehat{\mu}_Y,\widehat{\mu}\}$. We use each set of perturbed quantities to define a perturbed optimization problem, and solve for the corresponding weight vector $\widehat{\beta}^{[m]}$ as \begin{align} \widehat{\beta}^{[m]} =\operatorname*{arg\,min}_{\beta \in \widehat{\Omega}^{[m]}(\lambda)}\left[\tau^{[m]}(\beta)\right]^2, \quad\text{with}\quad\tau^{[m]}(\beta) = \widehat{\mu}_Y^{[m]}-(\widehat{\mu}^{[m]})^{\mkern-1.5mu\mathsf{T}}\beta, \end{align} where the perturbed uncertainty class $\widehat{\Omega}^{[m]}(\lambda)$, defined in the following (ref), is constructed using $\widehat{\gamma}^{[m]}$ and $\widehat{\Sigma}^{[m]}$. We show that there exists $m^*$ such that (ref) with $m=m^*$ nearly recovers (ref) and hence $(\widehat{\mu}^{[m^*]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m^*]}$ is nearly the same as $\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$; see Theorem (ref) in the Appendix. With such a nice $(\widehat{\mu}^{[m^*]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m^*]}$, the remaining uncertainty lies primarily in estimating $\mu_Y$, which can be addressed using standard inference methods, such as asymptotic normality. Figure (ref) presents a workflow of the perturbation-based inference procedure. \begin{figure}[ht] \caption{Workflow of the perturbation-based inference procedure.} \end{figure} We now provide full details of our proposal. Standard asymptotic results ensure that the estimators $\widehat{\Sigma}$ and $\widehat{\gamma}$ in (ref), as well as $\widehat{\mu}_Y$ and $\widehat{\mu}$ in (ref) admit the following asymptotic approximations, which can be established under mild conditions (e.g., i.i.d.\;data): \begin{align*}{\rm vecl}(\widehat{\Sigma}-\Sigma) \overset{d}{\approx} \mathcal{N}(0,\widehat{\mathbf{V}}_{\Sigma}),\quad \widehat{\gamma}-\gamma \overset{d}{\approx} \mathcal{N}(0,\widehat{\mathbf{V}}_{\gamma}),\quad \widehat{\mu}_Y-\mu_Y\overset{d}{\approx}\mathcal{N}(0,\widehat{{\rm V}}_Y),\quad \widehat{\mu}-\mu \overset{d}{\approx} \mathcal{N}(0,\widehat{\mathbf{V}}_{\mu}), \end{align*} where ${\rm vecl}(\widehat{\Sigma}-\Sigma)$ is the vector formed by stacking the columns of the lower triangular part of $\widehat{\Sigma} - \Sigma$, and $\overset{d}{\approx}$ denotes approximate equality in distribution, and $\widehat{\mathbf{V}}_{\Sigma}$, $\widehat{\mathbf{V}}_{\gamma}$, $\widehat{{\rm V}}_Y$, and $\widehat{\mathbf{V}}_{\mu}$ denote scaled estimated covariance matrices. We construct these estimators under an i.i.d.\;assumption as in (ref) in Appendix (ref), which also discusses covariance estimation in more general settings. In the following, we detail our two-step proposal: perturbation and aggregation. {\bf Step 1: Perturbation.} We begin by generating perturbed quantities related to the uncertainty class $\Omega$ in (ref) and objective function $\tau(\beta)$ in (ref). Specifically, conditioning on the observed data, we generate i.i.d.\ samples $\{\widehat{\Sigma}^{[m]},\widehat{\gamma}^{[m]}\}_{m=1}^M$, following \begin{align} &{\rm vecl}(\widehat{\Sigma}^{[m]}) \sim \mathcal{N}\left({\rm vecl}(\widehat{\Sigma}),\widehat{\mathbf{V}}_{\Sigma}+\|\widehat{\mathbf{V}}_{\Sigma}\|_{\max}\mathbf{I}\right), \quad \widehat{\gamma}^{[m]}\sim \mathcal{N}\left(\widehat{\gamma},\widehat{\mathbf{V}}_{\gamma}+\|\widehat{\mathbf{V}}_{\gamma}\|_{\max}\mathbf{I}\right). \end{align} To ensure symmetry of $\widehat{\Sigma}^{[m]}$, we impute the upper triangle part of each perturbed matrix $\widehat{\Sigma}^{[m]}$ by setting $\widehat{\Sigma}^{[m]}_{k,l} = \widehat{\Sigma}^{[m]}_{l,k}$ for $1 \leq l < k \leq N$. Also, conditioning on the observed data, we generate i.i.d.\ samples $\{\widehat{\mu}^{[m]},\widehat{\mu}_Y^{[m]}\}_{m=1}^M$, related to the objective function $\tau(\beta)$, as follows: \begin{align} &\widehat{\mu}_Y^{[m]} \sim \mathcal{N}\left(\widehat{\mu}_Y,\widehat{{\rm V}}_Y\right),\quad \widehat{\mu}^{[m]} \sim \mathcal{N}\left(\widehat{\mu},\widehat{\mathbf{V}}_{\mu}+\|\widehat{\mathbf{V}}_{\mu}\|_{\max}\mathbf{I}\right). \end{align} We add a diagonal matrix to the corresponding covariance matrix in the above generating process ensuring that the covariance matrix is positive definite. Specifically, we slightly enlarge the covariance matrices $\widehat{\mathbf{V}}_{\Sigma}$, $\widehat{\mathbf{V}}_{\gamma}$, and $\widehat{\mathbf{V}}_{\mu}$ to $\widehat{\mathbf{V}}_{\Sigma} + \|\widehat{\mathbf{V}}_{\Sigma}\|_{\max} \mathbf{I}$, $\widehat{\mathbf{V}}_{\gamma} + \|\widehat{\mathbf{V}}_{\gamma}\|_{\max}\mathbf{I}$, and $\widehat{\mathbf{V}}_{\mu} + \|\widehat{\mathbf{V}}_{\mu}\|_{\max} \mathbf{I}$, respectively. This adjustment mitigates numerical instability arising from near-singular covariance matrices, especially when $N$ is relatively large relative to $T_0$ or $T_1$.\footnote{While the main paper focuses on the regime where $T_0$ and $T_1$ are large relative to a fixed $N$, in practice, it is possible for $N$ to exceed either $T_0$ or $T_1$.} Throughout the paper, we use $p=1+N(N+5)/2$ to represent the total dimensionality of the quantities ${\rm vecl}(\Sigma)$, $\gamma$, $\mu_Y$, and $\mu$ that are used in the population optimization problem in (ref). For $1\leq m\leq M$, we substitute $\widehat{\Sigma}$ and $\widehat{\gamma}$ in (ref) with the perturbed ones $\widehat{\Sigma}^{[m]}$ and $\widehat{\gamma}^{[m]}$ and construct the perturbed uncertainty class $\widehat{\Omega}^{[m]}(\lambda)$ as \begin{align} \widehat{\Omega}^{[m]}(\lambda) = \left\{\beta \in \Delta^N: \|\widehat{\gamma}^{[m]} - \widehat{\Sigma}^{[m]}\beta\|_{\infty}\leq \lambda+\rho_M\right\}, \end{align} where $\rho_M \asymp\left[\log(\min\{T_0,T_1\})/M\right]^{1/p}/\sqrt{T_0}$ is a tuning parameter. We provide the details for the data-dependent selection of the tuning parameter $\rho_M$ after providing the full procedure of our proposed method. When there is no confusion, we denote $\widehat{\Omega}^{[m]}(\lambda)$ as $\widehat{\Omega}^{[m]}$. Using $\widehat{\mu}^{[m]}_Y$ and $\widehat{\mu}^{[m]}$ for the objective function $\tau^{[m]}(\beta)$, we construct the perturbed weights $\{\widehat{\beta}^{[m]}\}_{m=1}^M$ as the optimizer of the minimization problem (ref). We show that for sufficiently many perturbations, there exists an index $m^*$ such that the perturbed quantities, $\widehat{\Sigma}^{[m^*]}$, $\widehat{\gamma}^{[m^*]}$, $\widehat{\mu}^{[m^*]}_Y$, and $\widehat{\mu}^{[m^*]}$, closely retrieve the true population quantities $\Sigma$, $\gamma$, $\mu_Y$, and $\mu$, respectively; see Proposition (ref) in the Appendix. Thus, the perturbed optimization problem in (ref) becomes nearly equivalent to the population-level problem in (ref) when $m = m^*$. Finally, we construct the $m$-th perturbed estimator of $\tau^*$ as \begin{equation} \widehat{\tau}^{[m]} = \widehat{\mu}_Y - (\widehat{\mu}^{[m]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m]}, \end{equation} with $\widehat{\beta}^{[m]}$ defined in (ref). In (ref), we replace $\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}$ with its perturbed version $(\widehat{\mu}^{[m]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m]}$ while retaining $\widehat{\mu}_Y$ from $\widehat{\tau}$ in (ref). We decompose the estimation error $\widehat{\tau}^{[m]} - \tau^*$ as follows: \begin{align} \widehat{\tau}^{[m]} - \tau^* &= (\widehat{\mu}_Y-\mu_Y)-\left[(\widehat{\mu}^{[m]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m]}-\mu^{\mkern-1.5mu\mathsf{T}}\beta^*\right]. \end{align} This reveals that the estimation error consists of the asymptotically normal term $\widehat{\mu}_Y-\mu_Y$ and a perturbation error that becomes negligible for the index $m^*$ such that $(\widehat{\mu}^{[m^*]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m^*]}$ is nearly the same as $\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$. Simulation studies further confirm that such an index $m^*$ exists for which the non-regular term becomes negligible; see Appendix (ref). {\bf Step 2: Filtering and Aggregation.} We now discuss filtering out some inaccurate perturbations. Even though it is impossible to identify the best index $m^*$, our goal is to retain the perturbation $m^*$ and exclude those perturbations that are unlikely to be $m^*$. We screen out a small proportion of perturbations if they appear on the tails of the distributions in (ref) and (ref). We define normalized perturbed statistics as $\widehat{T}^{[m]}= \widehat{S}\widehat{U}^{[m]}$ for $m=1,\ldots,M$ where $\widehat{U}^{[m]}=(\widehat{\mu}^{[m]}_Y-\widehat{\mu}_Y, (\widehat{\mu}^{[m]}-\widehat{\mu})^{\mkern-1.5mu\mathsf{T}},({\rm vecl}(\widehat{\Sigma}^{[m]}-\widehat{\Sigma}))^{\mkern-1.5mu\mathsf{T}},(\widehat{\gamma}^{[m]}-\widehat{\gamma})^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$, $\widehat{S}= \text{diag}(\widehat{{\rm V}}_Y^{-1/2},D(\widehat{\mathbf{V}}_{\mu}) ,D(\widehat{\mathbf{V}}_{\Sigma}),D(\widehat{\mathbf{V}}_{\gamma}))$. Here, $D(A)$ denotes a vector of length equal to the diagonal dimension of $A$, with $[D(A)]_j=(A_{j,j}+\|A\|_{\max})^{-1/2}$. The vector $\widehat{U}^{[m]}$ collects the centered $m$-th perturbed quantities generated from the distributions in (ref) and (ref), and {$\widehat{S}$ is the corresponding diagonal scaling matrix.} Thus, $\widehat{T}^{[m]}$ represents the vector of normalized deviations between the perturbed quantities and their associated estimators. With $\widehat{T}^{[m]}$, we introduce the following index set $\mathbb{M}$ as \begin{align} &\mathbb{M}=\left\{1\leq m\leq M: \lambda_{\min}(\widehat{\Sigma}^{[m]})\geq 0, \quad \left\|\widehat{T}^{[m]}\right\|_{\infty} \leq 1.1z_{\alpha_0/(2p)}\right\}, \end{align} where $z_{q}$ is the upper $q$ quantile of the standard normal distribution, $\alpha_0 \in (0,0.01]$ is a prespecified constant to exclude extreme tail perturbations from (ref) and (ref), and the factor 1.1 adjusts for estimation error (any value greater than 1 may be used). In (ref), the index set $\mathbb{M}$ further excludes the $m$-th perturbation if {the minimum eigenvalue of $\widehat{\Sigma}^{[m]}$ is negative or} the maximum of the test statistics exceeds a specified threshold, which is chosen to adjust for multiple comparisons using the Bonferroni correction. Since $\lambda_{\min}(\Sigma)\geq 0$, it is reasonable to filter out $\widehat{\Sigma}^{[m]}$ when it has a negative eigenvalue. {Given the significance level $\alpha > \alpha_0$ and for each $m \in \mathbb{M}$,} we construct the $m$-th interval as \begin{align} {\rm Int}^{[m]} = \left[\widehat{\tau}^{[m]}-z_{\alpha'/2}\widehat{{\rm V}}_Y^{1/2},\;\widehat{\tau}^{[m]}+z_{\alpha'/2}\widehat{{\rm V}}_Y^{1/2}\right], \end{align} where $\alpha' = \alpha - \alpha_0$, $\widehat{\tau}^{[m]}$ is defined in (ref), and $\widehat{{\rm V}}_Y$ denotes the estimator of the variance of $\widehat{\mu}_Y$. If no feasible solution exists for (ref), we set ${\rm Int}^{[m]} = \varnothing.$ For each $m\in\mathbb{M}$, ${\rm Int}^{[m]}$ quantifies the uncertainty of $\widehat{\mu}_Y$ at the confidence level $\alpha$ using a standard inference method, while treating $(\widehat{\mu}^{[m]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m]}$ as being fixed. By the decomposition (ref) and discussion after that, there exists an index $m^* \in \mathbb{M}$ such that $(\widehat{\mu}^{[m^*]})^{\mkern-1.5mu\mathsf{T}}\widehat{\beta}^{[m^*]}$ as being nearly the same as $\mu^{\mkern-1.5mu\mathsf{T}}\beta^*$ and ${\rm Int}^{[m^*]}$ nearly serves as a level-$\alpha'$ CI for $\tau^*$. Since the specific identity of $m^*$ is unknown, we take the union and construct the CI for $\tau^*$: \begin{align} {\rm CI}_{\alpha} = \bigcup_{m\in \mathbb{M}}{\rm Int}^{[m]}, \end{align} with ${\rm Int}^{[m]}$ defined in (ref). We refer to ${\rm CI}_{\alpha}$ as a confidence interval even though $\cup_{m\in \mathbb{M}}{\rm Int}^{[m]}$ may not be an interval. We summarize our proposal in Algorithm (ref) in Appendix (ref). {\bf Tuning Parameter Selection.} Since we replace $\gamma$ and $\Sigma$ in (ref) with $\widehat{\gamma}^{[m]}$ and $\widehat{\Sigma}^{[m]}$ to construct $\widehat{\Omega}^{[m]}$ in (ref), we introduce the tuning parameter $\rho_M$ in (ref) to adjust for the error by this substitution. Proposition (ref) in the Appendix suggests selecting $\rho_M$ in a data-dependent way as $\rho_M = C_1[\log(\min\{T_0,T_1\})/M]^{1/p}/\sqrt{T_0}$ for some constant $C_1>0$. In practice, however, the exact value of $C_1$ is unknown. Analogous to selecting $C$ in $\rho$ in (ref), we initialize $C_1$ with a small value (e.g., $C_1=0.01$). However, a small value of $\rho_M$ may cause feasibility issues, as similar to those in choosing $\rho$ in (ref). We iteratively increase $C_1$ by a factor of 1.25 until a prespecified proportion of perturbed optimization problems (e.g., 10% by default) are feasible. Numerical studies in Appendix (ref) show robustness to this prespecified proportion: using 20% or 30% instead of our default 10% yields similar CI coverage and length. \section{Theoretical Justification} In this section, we provide theoretical justification for our methods. To facilitate theoretical analysis, we let $T_0$ and $T_1$ grow, and while we can consider growing $N$, we focus on a fixed-$N$ regime throughout the paper. In the main paper, we present results for $\lambda>0$, deferring the $\lambda=0$ case to Appendix (ref). We introduce assumptions on pre- and post-treatment data and error terms to characterize the convergence rates of $\widehat{\tau}$, beginning with the pre-treatment period. \begin{Assumption} For the pre-treatment control units' outcomes and errors $\{X_{t},u_t^{(0)}\}_{t=1}^{T_0}$ {in (ref)}, there exist constants $C_0,b>0$, independent of $T_0$ and $N$, such that as $T_0\to\infty$, \begin{align*} \mathbb{P}\left(\sup_{\beta\in\Delta^N}\left\|\frac{1}{T_0}\sum_{t=1}^{T_0}\left[X_tu_t^{(0)}+(X_tX_t^{\mkern-1.5mu\mathsf{T}}-\mathbb{E} X_tX_t^{\mkern-1.5mu\mathsf{T}})(\beta^{(0)}-\beta)\right]\right\|_{\infty}\leq \frac{C_0[\log(\max\{T_0,N\})]^{\frac{1+b}{2b}}}{\sqrt{T_0}}\right)\to 1. \end{align*} \end{Assumption} Assumption (ref) requires that the empirical averages of $X_tu_t^{(0)}$ and the deviation of the sample covariance matrix from its expectation remain uniformly controlled over $\beta^{(0)}-\beta$ with $\beta\in\Delta^N$. Since $\mathbb{E} X_tu_t^{(0)}=0$, this ensures these fluctuations vanish as $T_0$ grows, so that $\widehat{\Omega}$ behaves like its population counterpart $\Omega$. If $\{X_t,u_t^{(0)}\}_{t=1}^{T_0}$ are i.i.d.\ with fixed $N$, the assumption holds for all $b>0$. More generally, it holds for $\beta$-mixing pre-treatment data with exponential decay, where $b$ corresponds to the order of the $\beta$-mixing coefficients chernozhukov2021exact. Similar conditions appear in the SC literature for controlling prediction error without weight shifts chernozhukov2021exact,ben2021augmented. Next, we introduce an assumption on post-treatment error terms by defining the control-unit outcome errors $\nu_t = X_t - \mathbb{E}X_t$ for $t=T_0+1,\ldots,T$. \begin{Assumption} As $T_1\to\infty$, $\mu=T_1^{-1}\sum_{t=T_0+1}^{T}\mathbb{E} X_t$ is bounded and the error terms $\{ \epsilon_t\}_{t=T_0+1}^{T}$ satisfy $T_1^{-1/2}\sum_{t=T_0+1}^{T}\epsilon_{t} = O_p(1)$ with $ \epsilon_t= (\nu_{t}^{\mkern-1.5mu\mathsf{T}},v_t,u_t^{(1)})^{\mkern-1.5mu\mathsf{T}}$ where $v_t$ and $u_t^{(1)}$ are defined in (ref) and (ref), respectively. \end{Assumption} Assumption (ref) holds for i.i.d.\ post-treatment data. Assumption (ref) may hold under more general settings with dependent structures, such as strong mixing, provided that suitable conditions are satisfied billingsley2017probability. A similar assumption regarding the post-treatment data is used in the SC literature li2020statistical. The following theorem establishes the convergence rate of the estimator $\widehat{\tau}$ defined in (ref). \begin{Theorem} Suppose Assumptions (ref) and (ref) hold, and the tuning parameter $\rho$ used in (ref) satisfies $\rho =C[\log(\max\{T_0,N\})]^{\frac{1+b}{2b}}/{\sqrt{T_0}}$ for some positive constant $C\geq C_0$ with $C_0$ and $b$ in Assumption (ref). Then, for $\lambda>0$, $\widehat{\tau}$ in (ref) satisfies the following: \begin{align} \lim_{T_0,T_1\to\infty}\mathbb{P}\left(|\widehat{\tau}-\tau^*|\lesssim \left[\frac{[\log (\max\{T_0,N\})]^{\frac{1+b}{2b}}}{\sqrt{T_0}\cdot \lambda}\right]^{1/2} + \left[\frac{1}{\sqrt{T_1}}\right]^{1/2}\right)= 1. \end{align} \end{Theorem} The convergence rate (ref) of $|\widehat{\tau}-\tau^*|$ in Theorem (ref) has two components: the first term captures the error from estimating $\Omega$ by $\widehat{\Omega}$ and the second term reflects the error in estimating $\tau(\beta)$ by $\widehat{\mu}_Y-\widehat{\mu}^{\mkern-1.5mu\mathsf{T}}\beta$. The rate also applies as $\lambda\to 0$, but consistency of $\widehat{\tau}$ requires $\lambda$ large enough so that $[\log(\max\{T_0,N\})]^{\frac{1+b}{2b}}/(\sqrt{T_0}\cdot\lambda)\to 0$ as $T_0,T_1\to\infty$. The rate is slower than $1/\sqrt{\min\{T_0,T_1\}}$ since the objective in (ref) is convex but not strictly convex: $\mu\mu^{\mkern-1.5mu\mathsf{T}}$ is rank one, yielding a degenerate quadratic objective. This prevents the use of standard $M$-estimation theory and precludes the usual parametric rate, which is analogous to the slow rates in high-dimensional regression when restricted eigenvalue or strong convexity conditions fail buhlmann2011statistics,wainwright2019high. Finally, we theoretically justify our perturbation-based inference method. We impose the following assumption on the limiting distributions of ${\rm vecl}(\widehat{\Sigma})$, $\widehat{\gamma}$, $\widehat{\mu}_Y$, and $\widehat{\mu}$, as well as on the consistency of the corresponding covariance estimators. \begin{Assumption} ${\rm vecl}(\widehat{\Sigma})$, $\widehat{\gamma}$, $\widehat{\mu}_Y$, and $\widehat{\mu}$ admit the following asymptotic distributions: \begin{equation} \begin{aligned} &\mathbf{V}_{\Sigma}^{-1/2}\left({\rm vecl}(\widehat{\Sigma})-{\rm vecl}(\Sigma)\right)\xrightarrow{d}\mathcal{N}(0,\mathbf{I}),\quad \mathbf{V}_{\gamma}^{-1/2}\left(\widehat{\gamma}-\gamma\right)\xrightarrow{d}\mathcal{N}(0,\mathbf{I}),\quad \text{as}\quad T_0\rightarrow \infty, \\ &{\rm V}_Y^{-1/2}\left(\widehat{\mu}_Y-\mu_Y\right)\xrightarrow{d}\mathcal{N}(0,1),\quad \mathbf{V}_{\mu}^{-1/2}\left(\widehat{\mu}-\mu\right)\xrightarrow{d}\mathcal{N}(0,\mathbf{I}) \quad \text{as}\quad T_1\rightarrow \infty, \end{aligned} \end{equation} for some positive constant ${\rm V}_Y$ and positive definite matrices $\mathbf{V}_{\Sigma}$, $\mathbf{V}_{\gamma}$, and $\mathbf{V}_{\mu}$, where $T_1{\rm V}_Y$ and the elements of $T_0\mathbf{V}_{\Sigma}$, $T_0\mathbf{V}_{\gamma}$, and $T_1\mathbf{V}_{\mu}$ are bounded above and bounded away from zero. Furthermore, the rescaled covariance estimators $\widehat{\mathbf{V}}_{\Sigma}$, $\widehat{\mathbf{V}}_{\gamma}$, $\widehat{{\rm V}}_Y$, and $\widehat{\mathbf{V}}_{\mu}$ used in (ref) and (ref) satisfy consistency: $\|T_0(\widehat{\mathbf{V}}_{\Sigma}-\mathbf{V}_{\Sigma})\|_2\xrightarrow{p}0$ and $ \|T_0(\widehat{\mathbf{V}}_{\gamma}-\mathbf{V}_{\gamma})\|_2\xrightarrow{p}0$ as $T_0\to\infty$, as well as $|T_1(\widehat{{\rm V}}_Y-{\rm V}_Y)|\xrightarrow{p}0$ and $\|T_1(\widehat{\mathbf{V}}_{\mu}-\mathbf{V}_{\mu})\|_2\xrightarrow{p}0$ as $T_1\to\infty$. \end{Assumption} Assumption (ref) ensures the asymptotic normality of estimators for ${\rm vecl}(\widehat{\Sigma})$, $\widehat{\gamma}$, $\widehat{\mu}_Y$, and $\widehat{\mu}$, as well as the consistency of the corresponding covariance estimators. The central limit theorem (ref) holds when the pre- and post-treatment data are i.i.d.\;in the fixed $N$ regime, and more generally, under $\alpha$-mixing and stationarity billingsley2017probability. For consistency, the covariance estimators must match the data regime: the sample covariance estimators are appropriate under i.i.d.\;sampling, while HAC estimators accommodate weak dependence and stationarity newey1987simple,andrews1991heteroskedasticity; see (ref) and (ref) in Appendix (ref) for concrete formulas for covariance estimators under both data regimes. The following theorem establishes the coverage property of the proposed CI. \begin{Theorem} Suppose Assumption (ref) holds, and the tuning parameter $\rho_M$ used in (ref), satisfies $\rho_M = C_1[\log(\min\{T_0,T_1\})/M]^{1/p}/\sqrt{T_0}$ for some positive constant $C_1\geq (2/c^*(\alpha_0))^{1/p}$ with constants $c^*(\alpha_0)$ in (ref) in Appendix (ref) and $\alpha_0\in(0,0.01]$. Then, for $\alpha\in(\alpha_0,1)$ and $\lambda>0$, ${\rm CI}_{\alpha}$ in (ref) satisfies $\liminf_{T_0,T_1\rightarrow\infty}\liminf_{M\to\infty}\mathbb{P}(\tau^* \in {\rm CI}_{\alpha})\geq 1-\alpha.$ \end{Theorem} We note that Theorem (ref) establishes only a one-sided coverage guarantee, as our proposed perturbation-based method involves taking a union over $\mathbb{M}$. We now consider the length of our proposed CI. For controlling the interval length, we define a refined index set $\tilde{\mathbb{M}}$ based on $\mathbb{M}$ in (ref) as follows, \begin{align} \tilde{\mathbb{M}} = \mathbb{M}\cap\left\{1\leq m \leq M: \text{$\widehat{\Omega}^{[m]}(0)$ is non-empty}\right\}, \end{align} where $\widehat{\Omega}^{[m]}(\lambda)$ is defined in (ref) for $\lambda \geq 0$. The rationale for this refinement is that the best index $m^*$ should yield $\widehat{\Omega}^{[m^*]}(0)$ that nearly recovers $\Omega(0)$. Since $\Omega(0)$ contains the true pre-treatment weight $\beta^{(0)}$, $\widehat{\Omega}^{[m^*]}(0)$ should also contain $\beta^{(0)}$ and thus be non-empty. Thus, we filter out indices that violate (ref), as they are unlikely to correspond to $m^*$. Because $m^* \in \tilde{\mathbb{M}}$, the CI constructed using $\tilde{\mathbb{M}}$ retains the same coverage property in Theorem (ref). The following theorem presents the result on the length of the proposed CI. \begin{Theorem} Suppose Assumptions (ref), (ref), and (ref) hold, {and the tuning parameters $\rho$ used in (ref) and $\rho_M$ used in (ref) satisfy the same conditions in Theorems (ref) and (ref), respectively. For $\alpha\in(\alpha_0,1)$ with $\alpha_0\in[0,0.01)$, suppose ${\rm CI}_{\alpha}$ in (ref) is constructed using the refined index set $\tilde{\mathbb{M}}$ in (ref). Then, for $\lambda>0$, $\mathbf{L}({\rm CI}_{\alpha})$, the length of ${\rm CI}_{\alpha}$, satisfies the following:} $$\liminf_{T_0,T_1\to\infty}\liminf_{M\to\infty}\mathbb{P}\left(\mathbf{L}\left({\rm CI}_{\alpha}\right)\lesssim \left[\frac{[\log (\max\{T_0,N\})]^{\frac{1+b}{2b}}}{\sqrt{T_0}\cdot \lambda}\right]^{1/2} + \left[\frac{1}{\sqrt{T_1}}\right]^{1/2}\right)\geq 1-\alpha_0,$$ where $b>0$ is defined in Assumption (ref). \end{Theorem} Consistent with Theorem (ref), the CI length $\mathbf{L}({\rm CI}_{\alpha})$ in Theorem (ref) similarly does not attain the $1/\sqrt{\min\{T_0,T_1\}}$ rate due to the lack of strict convexity in the optimization problem (ref). It is unclear whether the constructed interval achieves the optimal precision property. Importantly, however, the validity of ${\rm CI}_{\alpha}$ does not rely on asymptotic normality of $\widehat{\tau}$. In Figure (ref), we evaluate its finite-sample performance and compare it to normality-based methods. The proposed CI is slightly longer but overall comparable to the oracle benchmark. \section{Simulation Studies} In this section, we evaluate DRoSC inference via numerical studies. Section (ref) describes the simulation settings, and Section (ref) assesses the proposed CIs from Section (ref). DRoSC estimation results are reported in Appendix (ref). \subsection{Simulation Setup} We describe the simulation design used to assess the accuracy of our inference procedure. Control units' outcomes $X_t$ are generated via a stationary AR(1) model with coefficient $\phi \in [0,1)$. Pre-treatment outcomes $\{X_t\}_{t=1}^{T_0}$ have mean $\mu_0=(0.8,1.2,...,0.8,1.2)^{\mkern-1.5mu\mathsf{T}}$ and equi-correlation covariance $\Sigma_0=(1 - \rho_0)\mathbf{I}_N + \rho_0 \mathbf{1}_N \mathbf{1}_N^{\mkern-1.5mu\mathsf{T}}$, while post-treatment outcomes $\{X_t\}_{t = T_0 + 1}^{T}$ have mean $\mu$ and covariance $\mathbf{I}_N$; see Appendix (ref) for details. The treated unit's potential outcomes follow (ref) with $\beta^{(0)} = (1/3\cdot \mathbf{1}_3^{\mkern-1.5mu\mathsf{T}}, \mathbf{0}_{N-3}^{\mkern-1.5mu\mathsf{T}})^{\mkern-1.5mu\mathsf{T}}$ and i.i.d.\ errors $u_t^{(0)}, u_t^{(1)} \sim \mathcal{N}(0, 1)$. We set $Y_{1,t}^{(1)} - Y_{1,t}^{(0)} = \tau + v_t$ with $v_t \overset{i.i.d.}{\sim} \mathcal{N}(0, 0.25^2)$, thus $\bar{\tau} = \tau$. We consider three settings by varying $\mu$, $\rho_0$, and $\beta^{(1)}$. Here, we present setting (S2), which violates (E1) due to high correlations and (E2) via small weight shifts: \begin{enumerate} • $\rho_0 = 0.95$, $\mu = \mu_0 + (0.6, 0.4, 0.2, \mathbf{0}_{N-3})^{\mkern-1.5mu\mathsf{T}}$, and $\beta^{(1)} = \beta^{(0)} + 0.05 \cdot (-1, \mathbf{0}_{N-2}^{\mkern-1.5mu\mathsf{T}}, 1)^{\mkern-1.5mu\mathsf{T}}$. \end{enumerate} The additional two settings (S1) and (S3) are detailed in Appendix (ref). We vary $\bar{\tau} \in \{-1.5,-1.4,\ldots,1.5\}$, $(T_0, T_1) \in \{25, 50\}\times\{25, 50\}$, and $\phi \in \{0, 0.5\}$, with fixed $N = 10$. In the main paper, we present results for setting (S2) with $\phi = 0$ (i.i.d.\ data) and $T_0 = T_1 = 25$. Additional results are provided in Appendix (ref) and are qualitatively similar to those in the main paper. \subsection{Inference} We evaluate the perturbation-based method ($1-\alpha = 0.95$) using 500 simulations for setting (S2). For cases sharing the same $\tau^*$ (e.g., $\tau^*=0$ when $\bar{\tau} \in\{-0.4,\ldots,0.1\}$), we report the minimum coverage and maximum CI length of each method. Our procedure is implemented with $M = 500$, and we denote our proposed CI in (ref) as Perturbed. We compare the Perturbed CI with CIs based on normality assumptions on DRoSC estimator $\widehat{\tau}$ in (ref): $(\widehat{\tau} - {\tau}^*) / {\rm SE}(\widehat{\tau}) \xrightarrow{d} \mathcal{N}(b^*, 1)$, where ${\rm SE}(\widehat{\tau})$ and $b^*$ denotes the standard error and bias of $\widehat{\tau}$, respectively. With $b^* = 0$, a valid oracle CI from the asymptotic normality is \begin{align} \left[\widehat{\tau} - z_{\alpha/2}\widehat{{\rm SE}}(\widehat{\tau}),\ \widehat{\tau} + z_{\alpha/2}\widehat{{\rm SE}}(\widehat{\tau}) \right], \end{align} where $\widehat{{\rm SE}}(\widehat{\tau})$ is the empirical standard error of $\widehat{\tau}$ across the 500 simulation replicates. We refer to this as the Normality CI. When $b^*\neq 0$, we construct an oracle bias-aware (OBA) CI following (6) of armstrong2023bias as the new benchmark: \begin{align} \left[\widehat{\tau} - \chi^*_{\alpha},\ \widehat{\tau} + \chi^*_{\alpha} \right],\quad \text{with}\quad \chi^*_{\alpha} = \widehat{{\rm SE}}(\widehat{\tau}) \left[{\text{cv}_{\alpha}\left(|\widehat{\mathbb{E}}\widehat{\tau} -\tau^*|^2 / \widehat{{\rm SE}}(\widehat{\tau})^2\right)}\right]^{1/2}, \end{align} where $\text{cv}_{\alpha}(B^2)$ is the $1-\alpha$ quantile of a non-central $\chi^2(1)$ distribution with non-centrality parameter $B^2$, $\widehat{\mathbb{E}}\widehat{\tau}$ is the empirical mean of $\widehat{\tau}$ across the 500 simulations, and $\widehat{{\rm SE}}(\widehat{\tau})$ is as in (ref). The OBA CI uses oracle knowledge of the bias $|\mathbb{E}\widehat{\tau} - \tau^*|$ as a rescaled term of the bias $b^*$, and is therefore infeasible, but we adopt it as the benchmark since it represents the best achievable CI given oracle information on the bias and standard error of $\widehat{\tau}.$ We present inference results with empirical coverages and mean lengths of CIs in Figure (ref). The Normality CI's coverage is near 0.95 when $\tau^*<0$ but drops to $0.9$ for some positive $\tau^*$ and below $0.9$ near zero, reflecting the non-regularity discussed in Section (ref). However, the OBA and Perturbed CIs remain uniformly valid, with the Perturbed CI somewhat conservative but comparable in length to the benchmark OBA. \begin{figure}[ht] \caption{Empirical coverages and interval lengths for setting (S2). $x$-axis plots the values of $\tau^*$. \texttt{Normality} refers to the Normality CI in (ref), \texttt{OBA} denotes the OBA CI in (ref), and \texttt{Perturbed} refers to the perturbation-based CI in (ref).} \end{figure} \section{Real Data Applications} In this section, we reanalyze the Basque Country case study of abadie2003economic using SC and our proposed method. The study examined the economic impact of terrorism on per capita GDP for the Basque Country and $N=16$ Spanish regions from 1955–1997 ($T=43$). Since terrorism affected only the Basque Country, it is considered treated, with the remaining regions as controls. The pre-treatment period is 1955–1969 ($T_0=15$); the post-treatment period is 1970–1997 ($T_1=28$). The original SC analysis matched the Basque Country to a weighted average of control regions using both pre-treatment outcomes and additional covariates; see abadie2003economic. As discussed in Section (ref), highly correlated controls and weight shifts raise concerns about SC’s stability and causal conclusions. To address these concerns, we apply DRoSC by estimating $\tau^*$ and constructing its 95% CIs using $M=500$, varying $\lambda\in\{0,0.001,\ldots,0.06\}$; for comparison, we also report the outcome-only SC point estimate (ref). Figure (ref) presents the results. For $\lambda=0$, our estimate is $\widehat{\tau}\approx -0.76$, compared to the SC estimate $\widehat{\tau}^{\rm SC}\approx -0.89$, yielding a more conservative effect even without weight shifts. As $\lambda$ increases, $\widehat{\tau}$ moves toward zero, reaching and remaining at zero for $\lambda\geq 0.054$. This highlights that (i) even without weight shifts, highly correlated controls permit alternative weights yielding more conservative estimates, and (ii) even small weight shifts can make the estimated effect zero. Consistently, for all $\lambda\ge0$, our CIs include zero, so the null hypothesis of no effect cannot be rejected; the slightly wider CIs are attributable to high correlations (Figure (ref)), which inflate the covariance matrix used in (ref) and (ref), yielding more dispersed perturbations and wider aggregated CIs in (ref). \begin{figure}[ht] \caption{Reanalysis of the Basque study. The black solid line shows $\widehat{\tau}$ in (ref) and the orange intervals are 95% CIs in (ref) for $\lambda\in\{0,0.015,\dots,0.06\}$; the blue horizontal dashed line shows $\widehat{\tau}^{\rm SC}$, and the red vertical dashed line marks 0.054, where $\widehat{\tau}$ first reaches 0. } \end{figure} \section{Conclusion and Discussion} The SC method is widely used for estimating the causal effects, but its validity relies on key assumptions: control units are not highly correlated and the treated-control relationship remains stable after treatment. To relax these conditions, we introduce a DRO-based causal estimand $\tau^*$ that equals {the time-averaged ATT $\bar{\tau}$} when $\bar{\tau}$ is identifiable otherwise serves as a conservative proxy for $\bar{\tau}$. We propose the DRoSC estimator, establish its convergence rate, and develop a perturbation-based method for valid inference despite non-standard limiting distributions. We robustify SC to accommodate weight shifts and ill-posed weight learning, and we propose the weight-robust treatment effect as an interpretable target under partial identification. We also bridge sensitivity analysis and DRO in the SC setting. Several avenues for future research remain. One direction is to extend the framework to incorporate covariates. While some applications use only outcomes, covariates can improve counterfactual construction abadie2022synthetic. In our framework, this can be done by incorporating covariates into the uncertainty class so that weights balance both outcomes and covariates. Additionally, extending the framework to staggered adoption settings, where units receive treatment at different times, is both relevant and challenging athey2022design, ben2022synthetic, cattaneo2025uncertainty. The same identification challenges—stemming from highly correlated controls and weight shifts—persist. While our analysis focuses on the time-averaged ATT for a single treated unit within a DRO framework, the approach naturally extends to staggered adoption settings by defining an alternative reward function that averages across varying adoption timing. \section*{Appendix} The appendices contain all proofs, additional methods, theories, and numerical results. \setstretch{1.73} \putbib