Synthetic Control Methods
Synthetic control methods estimate what would have happened to a treated unit without treatment by constructing a weighted comparison from untreated units. This guide develops the method around a single running example, California’s Proposition 99, and uses simulated panels with known answers to show what good and bad designs actually look like. It covers the design logic, diagnostics, a reproducible Synth workflow, generalized synthetic control in fect, augmented synthetic controls, and bounded or proportion outcomes.
- synthetic control
- causal inference
- program evaluation
- counterfactuals
- panel data
- policy evaluation
Synthetic Control Methods
Synthetic control methods estimate what would have happened to a treated unit without treatment by constructing a weighted comparison from untreated units. This guide develops the method around a single running example, California’s Proposition 99, and uses simulated panels with known answers to show what good and bad designs actually look like. It covers the design logic, diagnostics, a reproducible Synth workflow, generalized synthetic control in fect, augmented synthetic controls, and bounded or proportion outcomes.
Keywords
synthetic control, causal inference, program evaluation, counterfactuals, panel data, policy evaluation
Setup: R packages and Stata/Python engines
# Reproducibility: pacman installs any missing CRAN packages, then loads them all.
if (!requireNamespace("pacman", quietly = TRUE)) install.packages("pacman")
pacman::p_load(
Synth, fect, tidysynth, # synthetic control estimators (CRAN)
dplyr, tidyr, ggplot2, # data wrangling and plotting
Statamarkdown, # Stata engine for knitr
reticulate # Python engine for knitr
)
# augsynth is GitHub-only. Install it ONCE outside the render (see the
# Reproducibility appendix); we only attach it here, with no network call.
if (!requireNamespace("augsynth", quietly = TRUE)) {
stop("Package 'augsynth' is not installed. It is not on CRAN; install it once with\n",
' pacman::p_load_gh("ebenmichael/augsynth")\n',
' # or remotes::install_github("ebenmichael/augsynth")\n',
"See the Reproducibility appendix at the end of this document.")
}
library(augsynth)
# Python engine: this guide ships parallel R / Python / Stata code in every
# tabset (see the "Cross-language translation" appendix). A global
# RETICULATE_PYTHON in ~/.Rprofile would otherwise override use_virtualenv(), so
# clear it and bind this guide's own ./.venv (which has pysyncon, etc.).
Sys.unsetenv("RETICULATE_PYTHON")
reticulate::use_virtualenv(file.path(getwd(), ".venv"), required = TRUE)
# The Stata tabs write intermediate .dta handoffs; keep them out of the project
# root by pointing them at a scratch folder inside synthetic-control_files/.
dir.create("synthetic-control_files/stata", showWarnings = FALSE, recursive = TRUE)
# Stata graphs are exported to sc_figs/ and included as images (Statamarkdown
# does not auto-embed them). `graph export` errors with r(603) if the directory
# does not already exist, and sc_figs/ is gitignored, so create it here.
dir.create("sc_figs", showWarnings = FALSE, recursive = TRUE)
theme_set(theme_minimal(base_size = 12))
PAL <- c(treated = "#B22222", synth = "#1F4E79", donor = "grey70")Setup: Python imports (shared by every Python tab)
# Imported once; reticulate keeps a single Python session for the whole document,
# so every Python tab below can use np / pd / plt / pysyncon and any object a
# prior Python tab defined (fit_good, pl_keep, ...).
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pysyncon import Dataprep, Synth, AugSynth
plt.rcParams["figure.figsize"] = (8, 4.5)
PAL = {"treated": "#B22222", "synth": "#1F4E79", "donor": "0.7"}
# Simulation constants, mirroring the R "Simulation design" blocks.
J, T0, TT = 25, 20, 35Introduction
In 1988, California passed Proposition 99. The state raised its cigarette excise tax by 25 cents a pack and funded an anti-smoking campaign with the revenue. Cigarette sales in California then fell.
That last sentence is not yet a finding. Sales were already falling in California before 1988, and they were falling across the country. The question that matters is counterfactual: how many packs would Californians have bought if Proposition 99 had never passed? That quantity is never observed, and everything hard about this design lives in estimating it.
Why the Obvious Comparisons Fail
The natural instinct is to find a comparison state. The problem is that the treatment is idiosyncratic: it happened to exactly one unit, and that unit is unlike every other one.
- Pick one similar state. New York? It is more urban, has different tobacco politics, and its own tax history. Nevada? Different demographics entirely, and a tourism economy that moves cigarette sales for reasons unrelated to health policy. There is no second California sitting in the data, and any single choice is one referee question away from collapsing.
- Average all the other states. The average of 38 states is a perfectly well-defined object, and it is not California either. California smoked less than the national average and was on a steeper downward path well before 1988. Compare California to the national average and you will “find” an effect that started years before the policy did.
This is the recurring situation for aggregate interventions. The treatment is clean and unambiguous. The counterfactual is a mess, because every unit is unique in its own way.
The Synthetic Control Idea
Synthetic control methods resolve this by refusing to pick a comparison and refusing to weight everyone equally. Instead, they build one.
The method searches for a set of weights across the untreated states and forms a weighted average, a “synthetic California,” subject to one demand: the weighted average must reproduce California’s own pre-1988 path. Not California’s level in one year, and not its average, but its trajectory over the years leading up to the policy.
The logic is a matching argument moved from units to time. If a particular blend of Utah, Nevada, Montana, Colorado, and Connecticut tracked California’s cigarette sales year after year through the 1970s and 1980s, when neither California nor those states had Proposition 99, then that blend is a credible stand-in for what California was going to do next. After 1988, the blend keeps running as a no-treatment California. The two paths then differ by one thing that was not true before: California has the policy and the blend does not.
Synthetic Control in One Sentence
Build a comparison unit that reproduces the treated unit’s path before treatment, then ask whether the two paths separate after treatment.
Two features of this are worth naming immediately, because they drive everything that follows.
First, the pre-treatment fit is where the identifying argument becomes visible, not a preliminary step. In a regression you assume the comparison is valid. Here you display it. A reader can look at the pre-1988 fit and judge for themselves whether the counterfactual earned the right to speak about 1989. A good fit is necessary but not sufficient, a point we sharpen once we meet an estimator flexible enough to fit any past.
Second, the method can visibly fail. If no weighted blend of donors reproduces the treated unit’s pre-treatment path, you can see it in the figure. That is not a software problem; it is the design telling you that your data do not contain a plausible untreated version of the treated unit. Later in this guide we simulate exactly that case and watch the estimate go wrong in a specific, diagnosable way.
Is SCM the Right Design?
| Use SCM when… | Pause or use something else when… |
|---|---|
| One unit, or a small number of units, receives a well-defined intervention. | Treatment is vague, anticipated, reversible, or badly timed. |
| The unit is aggregate: a country, state, city, school district, firm, project area, or region. | Your data are primarily individual-level and rich in observed confounders. |
| You have a donor pool of untreated units observed before and after treatment. | Likely donors are affected by spillovers, related policies, or the treatment itself. |
| No single comparison unit is credible. | The treated unit is unlike every plausible donor before treatment. |
| The pre-treatment record is long enough to assess fit. | You have three pre-treatment periods and a hope. |
SCM is not magic. It is a way to learn whether your data contain a plausible untreated version of the treated unit. Sometimes the answer is no, and learning that early is a useful result.
Where SCM Gets Used
The same structure recurs across fields: one aggregate unit gets something, the counterfactual is contested, and the donor pool is the whole argument.
| Domain | Treated unit | Outcome | Donor pool | What makes the counterfactual hard |
|---|---|---|---|---|
| Tobacco policy | California, 1988 | Cigarette packs sold per capita | Other US states without comparable programs | California was already on a distinctive downward trend (Abadie, Diamond, and Hainmueller 2010) |
| Political conflict | Basque Country, from the 1970s | Real GDP per capita | Other Spanish regions | The region was economically atypical before the conflict (Abadie and Gardeazabal 2003) |
| Macro shocks | West Germany, 1990 | GDP per capita | Other industrialized economies | No country is “Germany without reunification” (Abadie, Diamond, and Hainmueller 2015) |
| Carbon offsets | Forest project area | Forest loss or deforestation share | Comparable non-project forest areas | The credited baseline is the estimand, and nearby areas may absorb displaced logging (West et al. 2023) |
| Economic shocks | The UK after the 2016 Brexit referendum | Real output | Other advanced economies | Markets and firms respond to the expected policy before it takes effect (Born et al. 2019) |
Theory
Potential Outcomes Setup
Let units be indexed by \(i = 1, \ldots, J + 1\). Unit \(1\) is treated; units \(2\) through \(J + 1\) are untreated donors. Time runs from \(t = 1, \ldots, T\), and treatment begins after \(T_0\). The periods \(1,\ldots,T_0\) are pre-treatment.
Each unit has two potential outcomes:
- \(Y_{it}(1)\) is the outcome under treatment.
- \(Y_{it}(0)\) is the outcome without treatment.
After treatment begins, we observe \(Y_{1t}(1)\) for the treated unit but never observe \(Y_{1t}(0)\). The period-specific treatment effect is therefore:
\[\alpha_{1t} = Y_{1t}(1) - Y_{1t}(0), \quad t > T_0.\]
SCM estimates the missing untreated outcome with a weighted average of donor outcomes:
\[\widehat{Y}_{1t}(0) = \sum_{j=2}^{J+1} w_jY_{jt},\]
so the estimated gap is:
\[\widehat{\alpha}_{1t} = Y_{1t} - \sum_{j=2}^{J+1} w_jY_{jt}, \quad t > T_0.\]
The gap plot is the figure everyone wants to see. A post-treatment gap is persuasive only when the pre-treatment fit makes it believable.
How Classical SCM Chooses a Counterfactual
Let \(X_1\) be a vector of pre-treatment predictors for the treated unit, and let \(X_0\) contain the same predictors for the \(J\) donor units. Classical SCM chooses donor weights \(W = (w_2, \ldots, w_{J+1})\) to minimize pre-treatment imbalance:
\[\widehat{W}(V) = \arg\min_W (X_1 - X_0W)^\top V(X_1 - X_0W),\]
subject to:
\[w_j \geq 0 \quad \text{and} \quad \sum_{j=2}^{J+1}w_j = 1.\]
\(V\) assigns importance to different predictors. The donor weights \(W\) tell us which units build the synthetic control. \(V\) decides what counts as a good match; \(W\) is the match.
The nonnegative, sum-to-one restrictions matter more than they look. They make the synthetic control a convex combination of real donors, which means the counterfactual can never leave the range spanned by units you actually observed. No extrapolation into territory no donor occupies. We return to this when outcomes are bounded, where the restriction stops being a philosophical preference and starts preventing impossible numbers.
That buys transparency, not automatic credibility. If synthetic California is mostly Utah and Nevada, a reader can decide whether that is substantively sensible.
A Predictor Set Can Defeat Itself
Feeding the optimizer the full pre-treatment outcome path alongside other covariates is a common instinct: more pre-treatment fit should be better. It can backfire. When the outcome path dominates, the \(V\)-optimization can drive the importance weights on the other covariates toward zero, so the “balance on predictors” diagnostic ends up checking a predictor set the optimizer never actually used (Kaul et al. 2022). A closely related failure mode is cherry-picking the predictor set or donor pool after seeing which choices produce an appealing gap (Ferman, Pinto, and Possebom 2020). Fix the predictor set and donor pool before looking at post-treatment results.
Identification: What You Must Defend
The identifying claim is simple to state and hard to earn: a donor combination that mimics the treated unit before treatment is a good approximation to what the treated unit would have done after treatment without treatment.
Under the latent-factor view, untreated outcomes reflect common time shocks, observed predictors, and unobserved factors with unit-specific loadings. Matching a long pre-treatment path can also match the treated unit’s exposure to those latent factors (Abadie 2021; Xu 2017). It is a helpful intuition, not a permission slip: the argument requires a pre-treatment window that is long relative to the number of latent factors, and a fit that is close to exact, not merely close-looking on a plot. As fit quality degrades, so does the bound on the resulting bias (Abadie 2021).
Before interpreting an estimate, defend these conditions:
| Design condition | Practical question |
|---|---|
| Clear treatment timing | When did behavior, policy, or anticipation actually begin? |
| Clean donor pool | Could any donor be affected by spillovers, related policies, displacement, or anticipation? |
| Plausible support | Can donors reproduce the treated unit before treatment without heroic extrapolation? |
| Stable measurement | Is the outcome measured consistently across units and over time? |
| Informative pre-period | Is there enough pre-treatment history to see persistent fit, not just accidental overlap? |
Diagnostics Are the Method
SCM lives or dies by diagnostics. The estimate is not one number; it is a design with receipts.
| Diagnostic | What it checks | What a bad result means |
|---|---|---|
| Pre-treatment RMSPE | Average mismatch before treatment | The synthetic control does not track the treated unit well enough. |
| Path plot | Treated and synthetic outcomes over time | Pre-period divergence undermines the post-period story. |
| Gap plot | Treated minus synthetic outcome | A post-treatment gap is only meaningful relative to pre-period gaps. |
| Predictor balance | Whether matching improves relevant predictors | Balance helps, but it cannot rescue a bad outcome path. |
| Donor weights | What the synthetic unit is made of | Large weights identify donors that need substantive scrutiny. |
| Placebo tests | Whether the treated gap is unusual among fake-treated donors | The observed pattern may be ordinary design noise. |
| Leave-one-out checks | Dependence on high-weight donors | The conclusion may rest on one donor. |
| Donor-pool sensitivity | Dependence on defensible exclusion rules | The result may be an artifact of one convenient pool. |
RMSPE Is a Diagnostic, Not a Universal Threshold
The pre-treatment root mean squared prediction error is:
\[RMSPE_{pre} = \sqrt{\frac{1}{T_0}\sum_{t=1}^{T_0} \left(Y_{1t} - \widehat{Y}_{1t}(0)\right)^2}.\]
Smaller is better. There is no sacred cutoff that converts a weak design into a strong one. Judge the error on the outcome’s scale, compare it with placebo fits, and show the path plot.
Placebos and Permutation Logic
For a placebo test, temporarily label each eligible donor as treated at the original treatment date, build its synthetic control, and compare its post-treatment gap with the treated unit’s gap. The point is not to manufacture a familiar regression p-value. It is to ask whether the treated unit looks unusual relative to units that received fake treatment.
For a chosen test statistic \(S\), a simple randomization-style summary is:
\[p = \frac{1 + \sum_{j \in \mathcal{D}} I\left\{S_j \geq S_{treated}\right\}} {1 + |\mathcal{D}|},\]
where \(\mathcal{D}\) includes placebo units with acceptable pre-treatment fit. State the statistic, the fit-screening rule, and whether the comparison is one- or two-sided.
Placebo permutation is not the only principled option, and it depends on having a donor pool large enough to support comparison. Chernozhukov, Wüthrich, and Zhu develop an exact conformal-inference procedure for synthetic control settings that does not rely on a large cross-section of placebo units (Chernozhukov, Wüthrich, and Zhu 2021). Cattaneo, Feng, and Titiunik develop prediction intervals for the synthetic control estimate itself, implemented in the scpi package (Cattaneo, Feng, and Titiunik 2021).
Seeing It Work and Fail
Everything above is a claim about what diagnostics buy you. This section tests those claims on simulated data, where we know the true answer and can check whether the diagnostics actually catch a bad design.
We generate panels from a factor model, the data-generating process the identification argument assumes:
\[Y_{it}(0) = \mu + \alpha_i + \lambda_i' f_t + \varepsilon_{it},\]
with two latent factors \(f_t\) and unit-specific loadings \(\lambda_i\). Two scenarios differ in exactly one respect:
- Good support: the treated unit’s factor loadings sit inside the range spanned by the donors, so a convex combination can reproduce them.
- Poor support: the treated unit’s loadings sit far outside the donor range, so no convex combination can.
In both cases we inject the same known treatment effect: a ramp that grows by \(-1.2\) units per period after \(T_0\), averaging \(-9.6\) over the post-treatment window. The design question is whether the diagnostics tell us when to trust the estimate.
The panels are generated once in R and written to synthetic-control_files/data/, so the R, Python, and Stata analyses below all read the same numbers. The generator is R-only (it depends on R’s random-number stream); expand the box if you want to see exactly how the data were built.
Simulation design (click to expand)
set.seed(20260720)
J <- 25 # number of donors
TT <- 35 # total periods
T0 <- 20 # last pre-treatment period
EFF <- -1.2
f1 <- scale(cumsum(rnorm(TT, 0, 1)))[, 1] # a stochastic trend
f2 <- sin(seq(0, 3 * pi, length.out = TT)) # a slower cycle
tau <- c(rep(0, T0), EFF * seq_len(TT - T0)) # the true effect path
make_panel <- function(lam_treated, seed = 20260720) {
set.seed(seed)
lam <- cbind(runif(J, 0, 1), runif(J, 0, 1)) # donor loadings in the unit square
alpha <- rnorm(J, 0, 3)
Y <- matrix(NA_real_, nrow = J + 1, ncol = TT)
for (i in 1:J) {
Y[i + 1, ] <- 120 + alpha[i] +
8 * lam[i, 1] * f1 + 8 * lam[i, 2] * f2 + rnorm(TT, 0, 1)
}
Y[1, ] <- 120 + 8 * lam_treated[1] * f1 + 8 * lam_treated[2] * f2 +
rnorm(TT, 0, 1) + tau
data.frame(
id = rep(1:(J + 1), each = TT),
name = rep(paste0("u", 1:(J + 1)), each = TT),
time = rep(1:TT, times = J + 1),
Y = as.vector(t(Y))
)
}
# Freeze the two scenarios and the ground truth to disk for the other languages.
dir_sim <- "synthetic-control_files/data"
write.csv(make_panel(c(0.5, 0.5)), file.path(dir_sim, "sim_good.csv"), row.names = FALSE)
write.csv(make_panel(c(2.6, 2.6)), file.path(dir_sim, "sim_bad.csv"), row.names = FALSE)
write.csv(data.frame(time = 1:TT, f1 = f1, f2 = f2, tau = tau),
file.path(dir_sim, "sim_truth.csv"), row.names = FALSE)
cat("True mean post-treatment effect:", mean(tau[(T0 + 1):TT]))True mean post-treatment effect: -9.6One helper does all the fitting. It matches on every pre-treatment outcome, which is the standard outcome-only synthetic control, and returns the pieces every diagnostic needs.
fit_scm <- function(df, treated, donors, T0, TT) {
sp <- lapply(1:T0, function(t) list("Y", t, "mean"))
prep <- dataprep(
foo = df, dependent = "Y",
unit.variable = "id", unit.names.variable = "name", time.variable = "time",
special.predictors = sp,
treatment.identifier = treated, controls.identifier = donors,
time.predictors.prior = 1:T0, time.optimize.ssr = 1:T0, time.plot = 1:TT
)
invisible(capture.output(res <- synth(prep)))
synth_path <- as.vector(prep$Y0plot %*% res$solution.w)
observed <- as.vector(prep$Y1plot)
gap <- observed - synth_path
list(
observed = observed, synth = synth_path, gap = gap,
w = setNames(as.vector(res$solution.w), colnames(prep$Y0plot)),
rmspe_pre = sqrt(mean(gap[1:T0]^2)),
att_post = mean(gap[(T0 + 1):TT])
)
}import pandas as pd
import numpy as np
from pysyncon import Dataprep, Synth
J, T0, TT = 25, 20, 35 # same constants as the simulation design
def fit_scm(df, treated, donors, T0, TT):
"""Outcome-only synthetic control, matching every pre-period outcome."""
dp = Dataprep(
foo=df, dependent="Y", unit_variable="name", time_variable="time",
treatment_identifier=treated, controls_identifier=donors,
predictors=[], predictors_op="mean",
special_predictors=[("Y", [t], "mean") for t in range(1, T0 + 1)],
time_predictors_prior=range(1, T0 + 1), time_optimize_ssr=range(1, T0 + 1),
)
s = Synth()
s.fit(dp)
wide = df.pivot(index="time", columns="name", values="Y").sort_index()
synth = wide[donors].values @ np.asarray(s.W)
obs = wide[treated].values
gap = obs - synth
return dict(observed=obs, synth=synth, gap=gap,
w=pd.Series(np.asarray(s.W), index=donors),
rmspe_pre=float(np.sqrt(np.mean(gap[:T0] ** 2))),
att_post=float(np.mean(gap[T0:])))* Stata runs `synth` per fit. Because each code block here is a fresh Stata
* session, the analyses below call `synth` inline rather than sharing a helper;
* the reusable shape of that call, for an outcome-only fit on a tsset panel
* with outcome Y over periods 1..T0, is:
*
* local pre ""
* forvalues t = 1/`T0' {
* local pre "`pre' Y(`t')"
* }
* synth Y `pre', trunit(<id>) trperiod(`=`T0'+1') ///
* mspeperiod(1(1)`T0') resultsperiod(1(1)`TT') keep(<file>, replace)
display "See the analysis blocks below for the inline synth calls."See the analysis blocks below for the inline synth calls.Two Fits, One Truth
panel_good <- read.csv("synthetic-control_files/data/sim_good.csv") # inside the hull
panel_bad <- read.csv("synthetic-control_files/data/sim_bad.csv") # far outside it
fit_good <- fit_scm(panel_good, 1, 2:(J + 1), T0, TT)
fit_bad <- fit_scm(panel_bad, 1, 2:(J + 1), T0, TT)
data.frame(
Scenario = c("Good support", "Poor support"),
`Pre-RMSPE` = round(c(fit_good$rmspe_pre, fit_bad$rmspe_pre), 2),
`Estimated ATT` = round(c(fit_good$att_post, fit_bad$att_post), 2),
`True ATT` = round(mean(tau[(T0 + 1):TT]), 2),
check.names = FALSE
) Scenario Pre-RMSPE Estimated ATT True ATT
1 Good support 0.67 -9.73 -9.6
2 Poor support 10.97 -1.48 -9.6panel_good = pd.read_csv("synthetic-control_files/data/sim_good.csv")
panel_bad = pd.read_csv("synthetic-control_files/data/sim_bad.csv")
donors = [f"u{i}" for i in range(2, J + 2)]
fit_good = fit_scm(panel_good, "u1", donors, T0, TT)
fit_bad = fit_scm(panel_bad, "u1", donors, T0, TT)
truth = pd.read_csv("synthetic-control_files/data/sim_truth.csv")
true_att = truth.loc[truth.time > T0, "tau"].mean()
print(pd.DataFrame({
"Scenario": ["Good support", "Poor support"],
"Pre-RMSPE": [round(fit_good["rmspe_pre"], 2), round(fit_bad["rmspe_pre"], 2)],
"Estimated ATT": [round(fit_good["att_post"], 2), round(fit_bad["att_post"], 2)],
"True ATT": round(true_att, 2),
}).to_string(index=False)) Scenario Pre-RMSPE Estimated ATT True ATT
Good support 0.67 -9.73 -9.6
Poor support 11.03 -2.47 -9.6tempname R
postfile `R' str12 scenario prermspe estatt using "synthetic-control_files/stata/_fitboth", replace
foreach s in good bad {
import delimited "synthetic-control_files/data/sim_`s'.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly synth y `pre', trunit(1) trperiod(21) ///
mspeperiod(1(1)20) resultsperiod(1(1)35) keep("synthetic-control_files/stata/_scm_`s'", replace)
preserve
use "synthetic-control_files/stata/_scm_`s'", clear
gen gap = _Y_treated - _Y_synthetic
gen gap2 = gap^2 if _time <= 20
quietly summarize gap2
scalar pr = sqrt(r(mean))
quietly summarize gap if _time > 20
post `R' ("`s'") (pr) (r(mean))
restore
}
postclose `R'
use "synthetic-control_files/stata/_fitboth", clear
list scenario prermspe estatt, noobs 9. preserve
10. use "synthetic-control_files/stata/_scm_`s'", clear
11. gen gap = _Y_treated - _Y_synthetic
12. gen gap2 = gap^2 if _time <= 20
13. quietly summarize gap2
14. scalar pr = sqrt(r(mean))
15. quietly summarize gap if _time > 20
16. post `R' ("`s'") (pr) (r(mean))
17. restore
18. }
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
+---------------------------------+
| scenario prermspe estatt |
|---------------------------------|
| good .6693835 -9.730545 |
| bad 10.97087 -1.478045 |
+---------------------------------+
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.32e-08 | -3.96e+02 | -6.46e+02 | 0.201 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.62e-11 | -3.96e+02 | -8.51e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.68e-03 | 3.31e-13 | -3.96e+02 | -4.05e+02 | 1.633 | 0.9950 | 2.07e-06
3 | 5.64e-05 | 1.11e-14 | -3.96e+02 | -3.96e+02 | 3.106 | 0.9664 | 1.74e-06
4 | 6.30e-06 | 1.24e-15 | -3.96e+02 | -3.96e+02 | 4.026 | 0.8884 | 6.53e-07
5 | 1.83e-06 | 3.61e-16 | -3.96e+02 | -3.96e+02 | 4.560 | 0.7090 | 5.53e-07
6 | 3.44e-07 | 6.78e-17 | -3.96e+02 | -3.96e+02 | 5.235 | 0.8123 | 6.59e-08
7 | 7.80e-09 | 1.57e-18 | -3.96e+02 | -3.96e+02 | 6.565 | 0.9773 | 1.97e-10
8 | 1.62e-10 | 5.22e-20 | -3.96e+02 | -3.96e+02 | 8.034 | 0.9792 | 7.76e-12
9 | 9.57e-13 | 6.42e-20 | -3.96e+02 | -3.96e+02 | 10.040 | 0.9941 | 2.11e-14
10 | 4.90e-15 | 6.69e-20 | -3.96e+02 | -3.96e+02 | 12.338 | 0.9950 | 1.50e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.43e-07 | -2.85e+02 | -5.34e+02 | 0.061 | 0.0000 | 1.00e+02
1 | 1.59e+00 | 3.38e-09 | -2.84e+02 | -7.34e+02 | 0.000 | 0.9763 | 4.65e-04
2 | 2.47e-02 | 5.23e-11 | -2.82e+02 | -3.18e+02 | 0.907 | 0.9845 | 2.64e-05
3 | 1.45e-03 | 3.08e-12 | -2.82e+02 | -2.85e+02 | 2.117 | 0.9411 | 2.12e-05
4 | 4.11e-04 | 8.71e-13 | -2.83e+02 | -2.84e+02 | 2.615 | 0.7172 | 3.02e-05
5 | 8.50e-05 | 1.80e-13 | -2.83e+02 | -2.83e+02 | 3.294 | 0.7933 | 5.07e-06
6 | 7.87e-07 | 1.67e-15 | -2.83e+02 | -2.83e+02 | 4.879 | 0.9907 | 2.05e-09
7 | 1.13e-08 | 2.39e-17 | -2.83e+02 | -2.83e+02 | 6.324 | 0.9857 | 1.28e-10
8 | 9.72e-11 | 2.38e-19 | -2.83e+02 | -2.83e+02 | 7.260 | 0.9914 | 1.66e-12
9 | 6.11e-13 | 1.03e-19 | -2.83e+02 | -2.83e+02 | 8.109 | 0.9937 | 1.02e-13
10 | 3.04e-15 | 7.93e-20 | -2.83e+02 | -2.83e+02 | 8.956 | 0.9949 | 9.36e-15
11 | 4.16e-16 | 9.21e-20 | -2.83e+02 | -2.83e+02 | 9.817 | 0.9950 | 1.30e-15
12 | 5.70e-16 | 8.45e-20 | -2.83e+02 | -2.83e+02 | 10.770 | 0.9950 | 1.79e-16
13 | 4.01e-16 | 9.27e-20 | -2.83e+02 | -2.83e+02 | 12.263 | 0.9950 | 2.00e-17
----------------------------------------------------------------------------------
optimization convergedThe numbers already tell the story. With good support the estimate lands within a rounding error of the truth. With poor support it is severely attenuated: the method recovers roughly a sixth of the true effect, so an analyst would conclude the intervention did almost nothing when in fact it did a great deal. And the pre-treatment RMSPE is roughly sixteen times larger, which is the tell.
path_df <- bind_rows(
data.frame(time = 1:TT, Observed = fit_good$observed,
Synthetic = fit_good$synth, panel = "Good support"),
data.frame(time = 1:TT, Observed = fit_bad$observed,
Synthetic = fit_bad$synth, panel = "Poor support")
) |>
pivot_longer(c(Observed, Synthetic), names_to = "series", values_to = "Y")
ggplot(path_df, aes(time, Y, colour = series, linetype = series)) +
geom_vline(xintercept = T0, linetype = "dotted", colour = "grey40") +
geom_line(linewidth = 0.8) +
facet_wrap(~panel) +
scale_colour_manual(values = c(Observed = PAL[["treated"]],
Synthetic = PAL[["synth"]])) +
labs(x = "Period", y = "Outcome", colour = NULL, linetype = NULL,
subtitle = "Dotted line marks the treatment date")
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5))
for ax, fit, title in [(axes[0], fit_good, "Good support"),
(axes[1], fit_bad, "Poor support")]:
t = range(1, TT + 1)
ax.plot(t, fit["observed"], color="#B22222", label="Observed")
ax.plot(t, fit["synth"], color="#1F4E79", linestyle="--", label="Synthetic")
ax.axvline(T0, color="grey", linestyle=":")
ax.set_title(title)
ax.set_xlabel("Period")
axes[0].set_ylabel("Outcome")
axes[0].legend(frameon=False)
fig.suptitle("Dotted line marks the treatment date")
plt.tight_layout()
plt.show()
* Uses the fit files saved by the fit-both block (_scm_good, _scm_bad).
use "synthetic-control_files/stata/_scm_good", clear
twoway (line _Y_treated _time, lcolor("178 34 34")) ///
(line _Y_synthetic _time, lcolor("31 78 121") lpattern(dash)), ///
xline(20, lpattern(dot)) title("Good support") ///
legend(order(1 "Observed" 2 "Synthetic")) name(g_good, replace) nodraw
use "synthetic-control_files/stata/_scm_bad", clear
twoway (line _Y_treated _time, lcolor("178 34 34")) ///
(line _Y_synthetic _time, lcolor("31 78 121") lpattern(dash)), ///
xline(20, lpattern(dot)) title("Poor support") ///
legend(off) name(g_bad, replace) nodraw
graph combine g_good g_bad, ycommon
quietly graph export "sc_figs/fig-paths-stata.svg", replace width(1600)
The left panel is what a usable design looks like: the two lines are indistinguishable before period 20, then separate. The right panel is what failure looks like, and it fails visibly and before the treatment date. Anyone can see it without knowing the true effect.
gap_df <- bind_rows(
data.frame(time = 1:TT, gap = fit_good$gap, panel = "Good support"),
data.frame(time = 1:TT, gap = fit_bad$gap, panel = "Poor support")
)
truth_df <- data.frame(time = 1:TT, gap = tau)
ggplot(gap_df, aes(time, gap)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = T0, linetype = "dotted", colour = "grey40") +
geom_line(data = truth_df, aes(time, gap), colour = "grey45",
linetype = "dashed", linewidth = 0.7) +
geom_line(colour = PAL[["treated"]], linewidth = 0.9) +
facet_wrap(~panel) +
labs(x = "Period", y = "Treated - synthetic",
subtitle = "Dashed grey line is the true effect")
tau = truth["tau"].values
fig, axes = plt.subplots(1, 2, figsize=(8, 4.5), sharey=True)
for ax, fit, title in [(axes[0], fit_good, "Good support"),
(axes[1], fit_bad, "Poor support")]:
t = range(1, TT + 1)
ax.axhline(0, color="grey")
ax.axvline(T0, color="grey", linestyle=":")
ax.plot(t, tau, color="grey", linestyle="--", label="True effect")
ax.plot(t, fit["gap"], color="#B22222", label="Estimated gap")
ax.set_title(title)
ax.set_xlabel("Period")
axes[0].set_ylabel("Treated - synthetic")
axes[0].legend(frameon=False)
fig.suptitle("Dashed grey line is the true effect")
plt.tight_layout()
plt.show()
foreach s in good bad {
use "synthetic-control_files/stata/_scm_`s'", clear
gen gap = _Y_treated - _Y_synthetic
gen tau = cond(_time > 20, -1.2 * (_time - 20), 0) // deterministic truth
twoway (line tau _time, lcolor(gs8) lpattern(dash)) ///
(line gap _time, lcolor("178 34 34")), ///
yline(0) xline(20, lpattern(dot)) title("`s' support") ///
legend(off) name(gap_`s', replace) nodraw
}
graph combine gap_good gap_bad, ycommon
quietly graph export "sc_figs/fig-gaps-stata.svg", replace width(1600) 2. use "synthetic-control_files/stata/_scm_`s'", clear
3. gen gap = _Y_treated - _Y_synthetic
4. gen tau = cond(_time > 20, -1.2 * (_time - 20), 0) // deterministic truth
5. twoway (line tau _time, lcolor(gs8) lpattern(dash)) ///
> (line gap _time, lcolor("178 34 34")), ///
> yline(0) xline(20, lpattern(dot)) title("`s' support") ///
> legend(off) name(gap_`s', replace) nodraw
6. }
In the good-support panel the estimated gap hugs zero before treatment and then tracks the true ramp. In the poor-support panel the gap is already swinging by ten units before anything happened, which is the diagnostic saying, in advance, that the post-treatment gap means nothing.
Placebo Inference
Now the permutation logic. We reassign the treatment to each donor in turn, refit, and compare the distribution of placebo gaps to the treated unit’s.
placebo_gaps <- function(panel, T0, TT, J) {
out <- lapply(2:(J + 1), function(u) {
donors <- setdiff(2:(J + 1), u)
f <- fit_scm(panel, u, donors, T0, TT)
data.frame(time = 1:TT, gap = f$gap, unit = paste0("u", u),
rmspe_pre = f$rmspe_pre)
})
bind_rows(out)
}
pl_good <- placebo_gaps(panel_good, T0, TT, J)
# Screen out placebos that never fit their own pre-period well.
cutoff <- 5 * fit_good$rmspe_pre
pl_keep <- pl_good |> filter(rmspe_pre <= cutoff)
n_dropped <- length(unique(pl_good$unit)) - length(unique(pl_keep$unit))
cat("Placebo units retained:", length(unique(pl_keep$unit)),
"| dropped for poor pre-fit:", n_dropped, "\n")Placebo units retained: 25 | dropped for poor pre-fit: 0 units = [f"u{i}" for i in range(2, J + 2)]
pl_records = []
for u in units:
others = [d for d in units if d != u]
f = fit_scm(panel_good, u, others, T0, TT)
pl_records.append({"unit": u, "gap": f["gap"], "rmspe_pre": f["rmspe_pre"]})
cutoff = 5 * fit_good["rmspe_pre"]
pl_keep = [p for p in pl_records if p["rmspe_pre"] <= cutoff]
print(f"Placebo units retained: {len(pl_keep)} | "
f"dropped for poor pre-fit: {len(pl_records) - len(pl_keep)}")Placebo units retained: 25 | dropped for poor pre-fit: 0* Refit with each donor (units 2..26) cast as the fake-treated unit, appending
* each one's gap series to _placebo.dta so a single failed fit is skipped, not fatal.
local first = 1
forvalues u = 2/26 {
import delimited "synthetic-control_files/data/sim_good.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
capture quietly synth y `pre', trunit(`u') trperiod(21) ///
mspeperiod(1(1)20) resultsperiod(1(1)35) keep("synthetic-control_files/stata/_pl_tmp", replace)
if _rc continue
use "synthetic-control_files/stata/_pl_tmp", clear
gen unit = `u'
gen gap = _Y_treated - _Y_synthetic
gen g2 = gap^2 if _time <= 20
quietly summarize g2
gen rmspe = sqrt(r(mean))
rename _time time
keep unit time gap rmspe
if `first' == 0 {
quietly append using "synthetic-control_files/stata/_placebo"
}
quietly save "synthetic-control_files/stata/_placebo", replace
local first = 0
}
* Screen placebos at 5x the treated unit's good-support pre-RMSPE (about 0.67).
use "synthetic-control_files/stata/_placebo", clear
egen tag = tag(unit)
count if tag & rmspe <= 5 * 0.67
display "Placebo units retained: " r(N) 9. if _rc continue
10. use "synthetic-control_files/stata/_pl_tmp", clear
11. gen unit = `u'
12. gen gap = _Y_treated - _Y_synthetic
13. gen g2 = gap^2 if _time <= 20
14. quietly summarize g2
15. gen rmspe = sqrt(r(mean))
16. rename _time time
17. keep unit time gap rmspe
18. if `first' == 0 {
19. quietly append using "synthetic-control_files/stata/_placebo"
20. }
21. quietly save "synthetic-control_files/stata/_placebo", replace
22. local first = 0
23. }
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(15 missing values generated)
25
Placebo units retained: 25
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.85e-08 | -4.04e+02 | -6.54e+02 | 0.209 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 9.26e-11 | -4.04e+02 | -8.59e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.83e-03 | 5.03e-13 | -4.04e+02 | -4.13e+02 | 1.632 | 0.9946 | 2.44e-06
3 | 7.04e-05 | 1.94e-14 | -4.04e+02 | -4.04e+02 | 3.045 | 0.9614 | 2.35e-06
4 | 1.17e-05 | 3.21e-15 | -4.04e+02 | -4.04e+02 | 3.699 | 0.8345 | 1.70e-06
5 | 3.01e-07 | 8.28e-17 | -4.04e+02 | -4.04e+02 | 4.668 | 0.9742 | 8.96e-09
6 | 6.00e-08 | 1.65e-17 | -4.04e+02 | -4.04e+02 | 5.311 | 0.8006 | 5.93e-08
7 | 3.09e-10 | 1.10e-19 | -4.04e+02 | -4.04e+02 | 6.477 | 0.9949 | 8.65e-12
8 | 4.86e-12 | 5.55e-20 | -4.04e+02 | -4.04e+02 | 7.828 | 0.9843 | 5.54e-12
9 | 3.62e-14 | 5.27e-20 | -4.04e+02 | -4.04e+02 | 9.428 | 0.9926 | 5.52e-14
10 | 7.47e-16 | 4.42e-20 | -4.04e+02 | -4.04e+02 | 11.682 | 0.9950 | 6.24e-16
11 | 6.75e-16 | 5.03e-20 | -4.04e+02 | -4.04e+02 | 13.995 | 0.9950 | 3.48e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.99e-08 | -3.80e+02 | -6.30e+02 | 0.183 | 0.0000 | 1.00e+02
1 | 5.37e-01 | 1.59e-10 | -3.80e+02 | -8.34e+02 | 0.000 | 0.9920 | 5.27e-05
2 | 6.30e-03 | 1.86e-12 | -3.80e+02 | -3.96e+02 | 1.372 | 0.9883 | 1.20e-05
3 | 5.63e-04 | 1.66e-13 | -3.80e+02 | -3.81e+02 | 2.418 | 0.9107 | 2.20e-05
4 | 3.54e-05 | 1.05e-14 | -3.80e+02 | -3.80e+02 | 3.253 | 0.9371 | 9.64e-07
5 | 9.23e-06 | 2.73e-15 | -3.80e+02 | -3.80e+02 | 3.824 | 0.7395 | 2.50e-06
6 | 2.38e-06 | 7.04e-16 | -3.80e+02 | -3.80e+02 | 4.376 | 0.7420 | 6.60e-07
7 | 1.52e-07 | 4.49e-17 | -3.80e+02 | -3.80e+02 | 5.225 | 0.9362 | 1.09e-08
8 | 2.71e-08 | 8.01e-18 | -3.80e+02 | -3.80e+02 | 5.966 | 0.8217 | 1.23e-08
9 | 4.99e-10 | 1.56e-19 | -3.80e+02 | -3.80e+02 | 7.516 | 0.9816 | 2.32e-11
10 | 2.58e-12 | 2.35e-20 | -3.80e+02 | -3.80e+02 | 9.735 | 0.9948 | 5.15e-14
11 | 1.30e-14 | 2.54e-20 | -3.80e+02 | -3.80e+02 | 12.033 | 0.9950 | 2.90e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.50e-08 | -3.80e+02 | -6.30e+02 | 0.183 | 0.0000 | 1.00e+02
1 | 6.79e-01 | 1.51e-10 | -3.80e+02 | -8.33e+02 | 0.000 | 0.9899 | 8.43e-05
2 | 4.93e-03 | 1.10e-12 | -3.79e+02 | -3.96e+02 | 1.358 | 0.9927 | 4.76e-06
3 | 3.26e-04 | 7.24e-14 | -3.79e+02 | -3.80e+02 | 2.536 | 0.9340 | 1.23e-05
4 | 6.58e-05 | 1.46e-14 | -3.80e+02 | -3.80e+02 | 3.150 | 0.7978 | 7.77e-06
5 | 2.49e-06 | 5.53e-16 | -3.80e+02 | -3.80e+02 | 4.081 | 0.9622 | 6.41e-08
6 | 2.95e-07 | 6.56e-17 | -3.80e+02 | -3.80e+02 | 4.831 | 0.8814 | 7.50e-08
7 | 5.69e-08 | 1.27e-17 | -3.80e+02 | -3.80e+02 | 5.464 | 0.8073 | 3.57e-08
8 | 5.67e-10 | 1.43e-19 | -3.80e+02 | -3.80e+02 | 7.188 | 0.9900 | 2.15e-11
9 | 2.92e-12 | 6.05e-20 | -3.80e+02 | -3.80e+02 | 9.446 | 0.9948 | 1.09e-13
10 | 1.46e-14 | 5.52e-20 | -3.80e+02 | -3.80e+02 | 11.743 | 0.9950 | 5.64e-16
11 | 4.52e-16 | 4.44e-20 | -3.80e+02 | -3.80e+02 | 14.062 | 0.9950 | 2.84e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.05e-08 | -4.03e+02 | -6.53e+02 | 0.209 | 0.0000 | 1.00e+02
1 | 3.56e-01 | 1.08e-10 | -4.03e+02 | -8.58e+02 | 0.000 | 0.9947 | 2.31e-05
2 | 2.42e-03 | 7.39e-13 | -4.03e+02 | -4.14e+02 | 1.586 | 0.9932 | 3.87e-06
3 | 1.20e-04 | 3.66e-14 | -4.03e+02 | -4.04e+02 | 2.891 | 0.9505 | 4.31e-06
4 | 1.47e-05 | 4.49e-15 | -4.03e+02 | -4.03e+02 | 3.626 | 0.8773 | 1.32e-06
5 | 1.66e-06 | 5.06e-16 | -4.03e+02 | -4.03e+02 | 4.406 | 0.8871 | 2.05e-07
6 | 3.66e-08 | 1.12e-17 | -4.03e+02 | -4.03e+02 | 5.565 | 0.9779 | 1.28e-09
7 | 3.26e-10 | 1.17e-19 | -4.03e+02 | -4.03e+02 | 7.328 | 0.9911 | 1.44e-11
8 | 1.64e-12 | 7.66e-20 | -4.03e+02 | -4.03e+02 | 9.585 | 0.9950 | 7.97e-14
9 | 8.18e-15 | 7.46e-20 | -4.03e+02 | -4.03e+02 | 11.883 | 0.9950 | 4.35e-16
10 | 3.88e-16 | 6.59e-20 | -4.03e+02 | -4.03e+02 | 14.180 | 0.9950 | 2.19e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.83e-08 | -4.19e+02 | -6.69e+02 | 0.226 | 0.0000 | 1.00e+02
1 | 4.84e-01 | 1.31e-10 | -4.19e+02 | -8.74e+02 | 0.000 | 0.9928 | 4.28e-05
2 | 5.15e-03 | 1.40e-12 | -4.19e+02 | -4.34e+02 | 1.455 | 0.9894 | 9.74e-06
3 | 4.08e-04 | 1.11e-13 | -4.19e+02 | -4.20e+02 | 2.554 | 0.9207 | 1.57e-05
4 | 6.69e-05 | 1.82e-14 | -4.19e+02 | -4.19e+02 | 3.229 | 0.8361 | 5.38e-06
5 | 2.73e-06 | 7.42e-16 | -4.19e+02 | -4.19e+02 | 4.207 | 0.9592 | 6.90e-08
6 | 1.20e-07 | 3.25e-17 | -4.19e+02 | -4.19e+02 | 5.009 | 0.9562 | 8.34e-09
7 | 5.52e-09 | 1.52e-18 | -4.19e+02 | -4.19e+02 | 6.094 | 0.9539 | 1.46e-09
8 | 4.00e-10 | 1.19e-19 | -4.19e+02 | -4.19e+02 | 7.202 | 0.9275 | 2.98e-10
9 | 2.39e-12 | 5.56e-20 | -4.19e+02 | -4.19e+02 | 9.105 | 0.9940 | 1.56e-13
10 | 1.19e-14 | 5.50e-20 | -4.19e+02 | -4.19e+02 | 11.389 | 0.9950 | 1.37e-15
11 | 5.62e-16 | 5.56e-20 | -4.19e+02 | -4.19e+02 | 13.699 | 0.9950 | 7.10e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.00e-08 | -4.03e+02 | -6.53e+02 | 0.209 | 0.0000 | 1.00e+02
1 | 3.44e-01 | 1.02e-10 | -4.03e+02 | -8.58e+02 | 0.000 | 0.9949 | 2.16e-05
2 | 1.88e-03 | 5.58e-13 | -4.03e+02 | -4.13e+02 | 1.623 | 0.9945 | 2.47e-06
3 | 6.67e-05 | 1.98e-14 | -4.03e+02 | -4.03e+02 | 3.071 | 0.9645 | 2.03e-06
4 | 1.48e-05 | 4.39e-15 | -4.03e+02 | -4.03e+02 | 3.670 | 0.7784 | 2.90e-06
5 | 1.54e-06 | 4.59e-16 | -4.03e+02 | -4.03e+02 | 4.490 | 0.8955 | 1.59e-07
6 | 1.55e-07 | 4.60e-17 | -4.03e+02 | -4.03e+02 | 5.421 | 0.8998 | 2.21e-08
7 | 2.51e-09 | 7.54e-19 | -4.03e+02 | -4.03e+02 | 7.107 | 0.9838 | 6.70e-11
8 | 1.47e-11 | 6.73e-20 | -4.03e+02 | -4.03e+02 | 9.170 | 0.9941 | 1.79e-13
9 | 7.36e-14 | 5.76e-20 | -4.03e+02 | -4.03e+02 | 11.468 | 0.9950 | 1.13e-15
10 | 5.90e-16 | 5.16e-20 | -4.03e+02 | -4.03e+02 | 13.748 | 0.9950 | 5.69e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.28e-08 | -4.18e+02 | -6.68e+02 | 0.225 | 0.0000 | 1.00e+02
1 | 5.55e-01 | 1.88e-10 | -4.18e+02 | -8.72e+02 | 0.000 | 0.9917 | 5.63e-05
2 | 7.06e-03 | 2.40e-12 | -4.18e+02 | -4.35e+02 | 1.392 | 0.9873 | 1.42e-05
3 | 5.74e-04 | 1.95e-13 | -4.18e+02 | -4.20e+02 | 2.480 | 0.9187 | 1.91e-05
4 | 1.10e-04 | 3.72e-14 | -4.18e+02 | -4.19e+02 | 3.107 | 0.8091 | 8.67e-06
5 | 2.94e-05 | 9.98e-15 | -4.18e+02 | -4.18e+02 | 3.652 | 0.7318 | 4.10e-06
6 | 2.25e-07 | 7.64e-17 | -4.18e+02 | -4.18e+02 | 4.766 | 0.9923 | 9.06e-10
7 | 1.88e-08 | 6.38e-18 | -4.18e+02 | -4.18e+02 | 5.764 | 0.9164 | 8.44e-09
8 | 9.59e-11 | 7.40e-20 | -4.18e+02 | -4.18e+02 | 7.770 | 0.9949 | 3.11e-12
9 | 4.92e-13 | 4.23e-20 | -4.18e+02 | -4.18e+02 | 10.027 | 0.9949 | 3.10e-14
10 | 2.55e-15 | 4.89e-20 | -4.18e+02 | -4.18e+02 | 12.324 | 0.9950 | 1.63e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.94e-08 | -3.85e+02 | -6.35e+02 | 0.189 | 0.0000 | 1.00e+02
1 | 9.55e-01 | 4.18e-10 | -3.85e+02 | -8.38e+02 | 0.000 | 0.9858 | 1.67e-04
2 | 1.11e-02 | 4.86e-12 | -3.85e+02 | -4.08e+02 | 1.220 | 0.9884 | 1.30e-05
3 | 8.50e-04 | 3.72e-13 | -3.85e+02 | -3.86e+02 | 2.332 | 0.9235 | 2.33e-05
4 | 3.23e-04 | 1.41e-13 | -3.85e+02 | -3.86e+02 | 2.719 | 0.6206 | 4.59e-05
5 | 9.11e-05 | 3.99e-14 | -3.85e+02 | -3.85e+02 | 3.248 | 0.7177 | 1.02e-05
6 | 2.23e-06 | 9.76e-16 | -3.85e+02 | -3.85e+02 | 4.481 | 0.9755 | 2.17e-08
7 | 1.28e-08 | 5.62e-18 | -3.85e+02 | -3.85e+02 | 6.403 | 0.9943 | 6.97e-11
8 | 6.40e-11 | 6.74e-20 | -3.85e+02 | -3.85e+02 | 8.702 | 0.9950 | 6.31e-13
9 | 3.20e-13 | 6.20e-20 | -3.85e+02 | -3.85e+02 | 11.000 | 0.9950 | 3.17e-15
10 | 1.72e-15 | 5.81e-20 | -3.85e+02 | -3.85e+02 | 13.299 | 0.9950 | 1.60e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.44e-08 | -4.25e+02 | -6.75e+02 | 0.232 | 0.0000 | 1.00e+02
1 | 5.87e-01 | 1.26e-10 | -4.25e+02 | -8.79e+02 | 0.000 | 0.9913 | 6.30e-05
2 | 5.71e-03 | 1.23e-12 | -4.25e+02 | -4.41e+02 | 1.420 | 0.9903 | 8.36e-06
3 | 4.45e-04 | 9.55e-14 | -4.25e+02 | -4.26e+02 | 2.527 | 0.9221 | 1.67e-05
4 | 2.22e-05 | 4.76e-15 | -4.25e+02 | -4.25e+02 | 3.363 | 0.9501 | 5.26e-07
5 | 4.27e-06 | 9.17e-16 | -4.25e+02 | -4.25e+02 | 4.046 | 0.8075 | 1.17e-06
6 | 1.09e-06 | 2.34e-16 | -4.25e+02 | -4.25e+02 | 4.608 | 0.7449 | 4.32e-07
7 | 2.17e-07 | 4.65e-17 | -4.25e+02 | -4.25e+02 | 5.287 | 0.8011 | 7.14e-08
8 | 1.05e-07 | 2.25e-17 | -4.25e+02 | -4.25e+02 | 5.612 | 0.5176 | 9.26e-08
9 | 5.77e-10 | 1.48e-19 | -4.25e+02 | -4.25e+02 | 7.331 | 0.9945 | 5.24e-12
10 | 2.90e-12 | 6.26e-20 | -4.25e+02 | -4.25e+02 | 9.622 | 0.9950 | 8.35e-14
11 | 1.45e-14 | 5.33e-20 | -4.25e+02 | -4.25e+02 | 11.920 | 0.9950 | 4.21e-16
12 | 3.74e-16 | 4.51e-20 | -4.25e+02 | -4.25e+02 | 14.252 | 0.9950 | 2.12e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.37e-08 | -4.30e+02 | -6.80e+02 | 0.236 | 0.0000 | 1.00e+02
1 | 5.51e-01 | 1.12e-10 | -4.30e+02 | -8.84e+02 | 0.000 | 0.9918 | 5.55e-05
2 | 8.21e-03 | 1.67e-12 | -4.29e+02 | -4.47e+02 | 1.382 | 0.9851 | 1.95e-05
3 | 8.14e-04 | 1.66e-13 | -4.29e+02 | -4.31e+02 | 2.383 | 0.9009 | 3.00e-05
4 | 1.94e-05 | 3.95e-15 | -4.30e+02 | -4.30e+02 | 3.269 | 0.9762 | 1.68e-07
5 | 3.07e-06 | 6.26e-16 | -4.30e+02 | -4.30e+02 | 4.035 | 0.8415 | 9.92e-07
6 | 4.51e-07 | 9.18e-17 | -4.30e+02 | -4.30e+02 | 4.683 | 0.8533 | 1.45e-07
7 | 1.43e-08 | 2.92e-18 | -4.30e+02 | -4.30e+02 | 6.010 | 0.9682 | 1.50e-09
8 | 8.13e-11 | 6.33e-20 | -4.30e+02 | -4.30e+02 | 8.143 | 0.9943 | 2.25e-12
9 | 4.06e-13 | 5.16e-20 | -4.30e+02 | -4.30e+02 | 10.440 | 0.9950 | 1.28e-14
10 | 2.10e-15 | 5.82e-20 | -4.30e+02 | -4.30e+02 | 12.737 | 0.9950 | 6.46e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.48e-08 | -3.81e+02 | -6.31e+02 | 0.185 | 0.0000 | 1.00e+02
1 | 6.03e-01 | 1.33e-10 | -3.81e+02 | -8.35e+02 | 0.000 | 0.9910 | 6.66e-05
2 | 4.43e-03 | 9.73e-13 | -3.81e+02 | -3.97e+02 | 1.395 | 0.9927 | 4.77e-06
3 | 2.82e-04 | 6.19e-14 | -3.81e+02 | -3.82e+02 | 2.589 | 0.9363 | 1.06e-05
4 | 3.89e-05 | 8.55e-15 | -3.81e+02 | -3.82e+02 | 3.317 | 0.8620 | 3.18e-06
5 | 6.50e-07 | 1.43e-16 | -3.81e+02 | -3.81e+02 | 4.325 | 0.9833 | 8.54e-09
6 | 1.60e-07 | 3.52e-17 | -3.81e+02 | -3.81e+02 | 4.902 | 0.7539 | 1.90e-07
7 | 1.30e-08 | 2.86e-18 | -3.81e+02 | -3.81e+02 | 5.697 | 0.9190 | 5.28e-09
8 | 3.59e-10 | 7.78e-20 | -3.81e+02 | -3.81e+02 | 7.117 | 0.9723 | 9.78e-11
9 | 4.51e-12 | 4.22e-20 | -3.81e+02 | -3.81e+02 | 8.620 | 0.9874 | 7.67e-13
10 | 3.21e-14 | 3.61e-20 | -3.81e+02 | -3.81e+02 | 9.900 | 0.9929 | 7.64e-15
11 | 7.03e-16 | 4.26e-20 | -3.81e+02 | -3.81e+02 | 11.964 | 0.9950 | 1.99e-16
12 | 4.30e-16 | 5.30e-20 | -3.81e+02 | -3.81e+02 | 14.284 | 0.9950 | 1.72e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.29e-08 | -4.37e+02 | -6.86e+02 | 0.243 | 0.0000 | 1.00e+02
1 | 9.68e-01 | 3.30e-10 | -4.36e+02 | -8.89e+02 | 0.000 | 0.9856 | 1.72e-04
2 | 2.44e-02 | 8.32e-12 | -4.36e+02 | -4.65e+02 | 1.173 | 0.9748 | 6.15e-05
3 | 2.56e-03 | 8.71e-13 | -4.36e+02 | -4.39e+02 | 2.143 | 0.8953 | 5.63e-05
4 | 4.73e-04 | 1.61e-13 | -4.36e+02 | -4.37e+02 | 2.706 | 0.8150 | 1.85e-05
5 | 2.32e-05 | 7.90e-15 | -4.37e+02 | -4.37e+02 | 3.825 | 0.9510 | 3.46e-07
6 | 4.95e-06 | 1.69e-15 | -4.37e+02 | -4.37e+02 | 4.468 | 0.7865 | 5.13e-07
7 | 2.83e-08 | 9.62e-18 | -4.37e+02 | -4.37e+02 | 6.001 | 0.9943 | 8.04e-11
8 | 1.41e-10 | 7.97e-20 | -4.37e+02 | -4.37e+02 | 8.294 | 0.9950 | 1.81e-12
9 | 7.07e-13 | 4.96e-20 | -4.37e+02 | -4.37e+02 | 10.591 | 0.9950 | 9.21e-15
10 | 3.54e-15 | 4.64e-20 | -4.37e+02 | -4.37e+02 | 12.893 | 0.9950 | 4.64e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.94e-08 | -3.99e+02 | -6.49e+02 | 0.205 | 0.0000 | 1.00e+02
1 | 4.59e-01 | 1.33e-10 | -3.99e+02 | -8.54e+02 | 0.000 | 0.9932 | 3.85e-05
2 | 3.13e-03 | 9.05e-13 | -3.99e+02 | -4.12e+02 | 1.504 | 0.9932 | 3.97e-06
3 | 1.33e-04 | 3.84e-14 | -3.99e+02 | -4.00e+02 | 2.875 | 0.9576 | 3.79e-06
4 | 4.20e-06 | 1.21e-15 | -3.99e+02 | -3.99e+02 | 3.680 | 0.9683 | 8.89e-08
5 | 1.10e-06 | 3.18e-16 | -3.99e+02 | -3.99e+02 | 4.239 | 0.7378 | 9.96e-07
6 | 3.44e-07 | 9.96e-17 | -3.99e+02 | -3.99e+02 | 4.759 | 0.6873 | 3.95e-07
7 | 6.15e-08 | 1.78e-17 | -3.99e+02 | -3.99e+02 | 5.497 | 0.8215 | 3.79e-08
8 | 3.07e-10 | 1.10e-19 | -3.99e+02 | -3.99e+02 | 7.524 | 0.9950 | 5.27e-12
9 | 1.54e-12 | 4.72e-20 | -3.99e+02 | -3.99e+02 | 9.823 | 0.9950 | 4.95e-14
10 | 7.69e-15 | 4.27e-20 | -3.99e+02 | -3.99e+02 | 12.121 | 0.9950 | 2.49e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.26e-08 | -4.29e+02 | -6.79e+02 | 0.236 | 0.0000 | 1.00e+02
1 | 8.61e-01 | 2.90e-10 | -4.29e+02 | -8.82e+02 | 0.000 | 0.9872 | 1.36e-04
2 | 2.06e-02 | 6.92e-12 | -4.29e+02 | -4.56e+02 | 1.202 | 0.9761 | 5.39e-05
3 | 2.30e-03 | 7.75e-13 | -4.29e+02 | -4.32e+02 | 2.145 | 0.8881 | 5.91e-05
4 | 1.26e-04 | 4.25e-14 | -4.29e+02 | -4.30e+02 | 2.922 | 0.9452 | 1.55e-06
5 | 5.03e-06 | 1.69e-15 | -4.29e+02 | -4.29e+02 | 4.104 | 0.9602 | 1.36e-07
6 | 1.57e-06 | 5.28e-16 | -4.29e+02 | -4.29e+02 | 4.499 | 0.6876 | 5.79e-07
7 | 2.52e-07 | 8.45e-17 | -4.29e+02 | -4.29e+02 | 5.232 | 0.8398 | 5.96e-08
8 | 4.31e-08 | 1.45e-17 | -4.29e+02 | -4.29e+02 | 6.010 | 0.8288 | 1.26e-08
9 | 2.20e-10 | 8.96e-20 | -4.29e+02 | -4.29e+02 | 8.186 | 0.9949 | 1.81e-12
10 | 1.10e-12 | 5.28e-20 | -4.29e+02 | -4.29e+02 | 10.484 | 0.9950 | 1.16e-14
11 | 5.47e-15 | 6.18e-20 | -4.29e+02 | -4.29e+02 | 12.783 | 0.9950 | 5.84e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.74e-08 | -3.76e+02 | -6.26e+02 | 0.179 | 0.0000 | 1.00e+02
1 | 7.09e-01 | 1.83e-10 | -3.76e+02 | -8.30e+02 | 0.000 | 0.9895 | 9.20e-05
2 | 7.59e-03 | 1.96e-12 | -3.76e+02 | -3.94e+02 | 1.303 | 0.9893 | 1.04e-05
3 | 5.89e-04 | 1.52e-13 | -3.76e+02 | -3.77e+02 | 2.411 | 0.9223 | 1.92e-05
4 | 6.92e-05 | 1.79e-14 | -3.76e+02 | -3.76e+02 | 3.131 | 0.8826 | 3.41e-06
5 | 1.29e-06 | 3.33e-16 | -3.76e+02 | -3.76e+02 | 4.072 | 0.9814 | 1.61e-08
6 | 3.83e-07 | 9.90e-17 | -3.76e+02 | -3.76e+02 | 4.591 | 0.7028 | 4.93e-07
7 | 6.89e-08 | 1.78e-17 | -3.76e+02 | -3.76e+02 | 5.303 | 0.8200 | 5.35e-08
8 | 6.68e-10 | 1.64e-19 | -3.76e+02 | -3.76e+02 | 6.981 | 0.9903 | 2.92e-11
9 | 3.35e-12 | 4.61e-20 | -3.76e+02 | -3.76e+02 | 9.270 | 0.9950 | 1.64e-13
10 | 1.67e-14 | 3.63e-20 | -3.76e+02 | -3.76e+02 | 11.568 | 0.9950 | 8.37e-16
11 | 4.71e-16 | 4.07e-20 | -3.76e+02 | -3.76e+02 | 13.892 | 0.9950 | 4.22e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.93e-08 | -4.43e+02 | -6.93e+02 | 0.249 | 0.0000 | 1.00e+02
1 | 7.76e-01 | 2.23e-10 | -4.42e+02 | -8.96e+02 | 0.000 | 0.9885 | 1.10e-04
2 | 1.48e-02 | 4.27e-12 | -4.42e+02 | -4.65e+02 | 1.272 | 0.9809 | 3.39e-05
3 | 1.21e-03 | 3.49e-13 | -4.42e+02 | -4.44e+02 | 2.353 | 0.9182 | 2.73e-05
4 | 2.57e-04 | 7.37e-14 | -4.42e+02 | -4.43e+02 | 2.860 | 0.7888 | 1.51e-05
5 | 3.30e-05 | 9.47e-15 | -4.42e+02 | -4.42e+02 | 3.709 | 0.8716 | 1.71e-06
6 | 6.57e-07 | 1.89e-16 | -4.42e+02 | -4.42e+02 | 5.189 | 0.9801 | 5.70e-09
7 | 3.33e-09 | 9.50e-19 | -4.42e+02 | -4.42e+02 | 7.460 | 0.9949 | 1.22e-11
8 | 1.66e-11 | 6.10e-20 | -4.42e+02 | -4.42e+02 | 9.757 | 0.9950 | 6.37e-14
9 | 8.33e-14 | 4.19e-20 | -4.42e+02 | -4.42e+02 | 12.055 | 0.9950 | 3.21e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.73e-08 | -3.63e+02 | -6.13e+02 | 0.163 | 0.0000 | 1.00e+02
1 | 9.10e-01 | 2.34e-10 | -3.63e+02 | -8.16e+02 | 0.000 | 0.9865 | 1.52e-04
2 | 2.00e-02 | 5.15e-12 | -3.62e+02 | -3.89e+02 | 1.130 | 0.9780 | 4.62e-05
3 | 2.38e-03 | 6.13e-13 | -3.62e+02 | -3.66e+02 | 2.047 | 0.8810 | 6.68e-05
4 | 9.72e-05 | 2.50e-14 | -3.63e+02 | -3.64e+02 | 2.834 | 0.9592 | 9.09e-07
5 | 2.08e-05 | 5.34e-15 | -3.63e+02 | -3.63e+02 | 3.479 | 0.7864 | 4.19e-06
6 | 1.21e-07 | 3.11e-17 | -3.63e+02 | -3.63e+02 | 4.567 | 0.9942 | 6.77e-10
7 | 3.37e-09 | 8.57e-19 | -3.63e+02 | -3.63e+02 | 5.589 | 0.9721 | 1.28e-09
8 | 1.14e-09 | 2.91e-19 | -3.63e+02 | -3.63e+02 | 6.085 | 0.6607 | 1.90e-08
9 | 5.72e-12 | 4.25e-20 | -3.63e+02 | -3.63e+02 | 7.694 | 0.9950 | 1.24e-12
10 | 7.24e-14 | 5.15e-20 | -3.63e+02 | -3.63e+02 | 9.270 | 0.9874 | 1.95e-13
11 | 6.90e-16 | 4.37e-20 | -3.63e+02 | -3.63e+02 | 11.486 | 0.9949 | 8.48e-16
12 | 4.85e-16 | 4.46e-20 | -3.63e+02 | -3.63e+02 | 13.789 | 0.9950 | 4.92e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.10e-08 | -3.87e+02 | -6.37e+02 | 0.191 | 0.0000 | 1.00e+02
1 | 7.36e-01 | 2.30e-10 | -3.87e+02 | -8.41e+02 | 0.000 | 0.9891 | 9.92e-05
2 | 5.48e-03 | 1.71e-12 | -3.87e+02 | -4.04e+02 | 1.339 | 0.9926 | 5.07e-06
3 | 2.75e-04 | 8.61e-14 | -3.87e+02 | -3.88e+02 | 2.636 | 0.9497 | 7.58e-06
4 | 5.27e-05 | 1.65e-14 | -3.87e+02 | -3.87e+02 | 3.208 | 0.8087 | 5.63e-06
5 | 1.05e-05 | 3.28e-15 | -3.87e+02 | -3.87e+02 | 3.878 | 0.8009 | 1.63e-06
6 | 1.81e-06 | 5.67e-16 | -3.87e+02 | -3.87e+02 | 4.605 | 0.8272 | 2.62e-07
7 | 2.82e-08 | 8.81e-18 | -3.87e+02 | -3.87e+02 | 6.176 | 0.9844 | 3.86e-10
8 | 4.64e-10 | 1.69e-19 | -3.87e+02 | -3.87e+02 | 7.689 | 0.9835 | 1.16e-11
9 | 4.44e-12 | 5.59e-20 | -3.87e+02 | -3.87e+02 | 8.872 | 0.9904 | 1.21e-13
10 | 2.28e-14 | 6.64e-20 | -3.87e+02 | -3.87e+02 | 10.163 | 0.9949 | 2.28e-15
11 | 5.04e-16 | 5.86e-20 | -3.87e+02 | -3.87e+02 | 12.406 | 0.9950 | 1.10e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.38e-08 | -4.06e+02 | -6.56e+02 | 0.212 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.92e-11 | -4.06e+02 | -8.61e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.69e-03 | 3.47e-13 | -4.06e+02 | -4.15e+02 | 1.644 | 0.9950 | 2.08e-06
3 | 6.56e-05 | 1.35e-14 | -4.06e+02 | -4.06e+02 | 3.053 | 0.9611 | 2.34e-06
4 | 6.43e-06 | 1.32e-15 | -4.06e+02 | -4.06e+02 | 3.973 | 0.9019 | 5.83e-07
5 | 1.39e-06 | 2.87e-16 | -4.06e+02 | -4.06e+02 | 4.469 | 0.7832 | 3.51e-07
6 | 2.46e-07 | 5.07e-17 | -4.06e+02 | -4.06e+02 | 5.196 | 0.8234 | 7.36e-08
7 | 3.40e-09 | 6.80e-19 | -4.06e+02 | -4.06e+02 | 6.866 | 0.9862 | 8.20e-11
8 | 2.36e-11 | 4.98e-20 | -4.06e+02 | -4.06e+02 | 8.847 | 0.9931 | 4.41e-13
9 | 1.18e-13 | 6.07e-20 | -4.06e+02 | -4.06e+02 | 11.136 | 0.9950 | 2.40e-15
10 | 7.34e-16 | 5.69e-20 | -4.06e+02 | -4.06e+02 | 13.430 | 0.9950 | 1.23e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.01e-08 | -3.91e+02 | -6.41e+02 | 0.195 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 1.00e-10 | -3.91e+02 | -8.46e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 2.32e-03 | 6.92e-13 | -3.91e+02 | -4.01e+02 | 1.588 | 0.9931 | 3.93e-06
3 | 1.62e-04 | 4.83e-14 | -3.91e+02 | -3.91e+02 | 2.743 | 0.9302 | 8.29e-06
4 | 3.25e-05 | 9.69e-15 | -3.91e+02 | -3.91e+02 | 3.384 | 0.7992 | 4.90e-06
5 | 8.12e-07 | 2.42e-16 | -3.91e+02 | -3.91e+02 | 4.314 | 0.9750 | 1.68e-08
6 | 1.80e-07 | 5.37e-17 | -3.91e+02 | -3.91e+02 | 4.930 | 0.7782 | 1.61e-07
7 | 7.04e-09 | 2.11e-18 | -3.91e+02 | -3.91e+02 | 6.178 | 0.9609 | 1.17e-09
8 | 3.77e-11 | 3.70e-20 | -3.91e+02 | -3.91e+02 | 8.344 | 0.9946 | 1.23e-12
9 | 1.89e-13 | 3.69e-20 | -3.91e+02 | -3.91e+02 | 10.643 | 0.9950 | 7.33e-15
10 | 1.01e-15 | 3.33e-20 | -3.91e+02 | -3.91e+02 | 12.941 | 0.9950 | 3.69e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.34e-08 | -3.97e+02 | -6.47e+02 | 0.202 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.68e-11 | -3.97e+02 | -8.52e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.68e-03 | 3.34e-13 | -3.97e+02 | -4.06e+02 | 1.634 | 0.9950 | 2.07e-06
3 | 4.19e-05 | 8.32e-15 | -3.97e+02 | -3.97e+02 | 3.236 | 0.9751 | 9.56e-07
4 | 1.07e-05 | 2.12e-15 | -3.97e+02 | -3.97e+02 | 3.817 | 0.7450 | 2.60e-06
5 | 3.47e-06 | 6.89e-16 | -3.97e+02 | -3.97e+02 | 4.305 | 0.6753 | 1.12e-06
6 | 5.25e-07 | 1.04e-16 | -3.97e+02 | -3.97e+02 | 5.099 | 0.8488 | 7.65e-08
7 | 3.53e-08 | 7.03e-18 | -3.97e+02 | -3.97e+02 | 6.084 | 0.9326 | 2.41e-09
8 | 8.18e-10 | 1.56e-19 | -3.97e+02 | -3.97e+02 | 7.569 | 0.9769 | 2.91e-11
9 | 4.39e-12 | 7.11e-20 | -3.97e+02 | -3.97e+02 | 9.812 | 0.9946 | 5.10e-14
10 | 2.19e-14 | 7.00e-20 | -3.97e+02 | -3.97e+02 | 12.110 | 0.9950 | 2.53e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.27e-08 | -4.17e+02 | -6.67e+02 | 0.224 | 0.0000 | 1.00e+02
1 | 3.68e-01 | 6.95e-11 | -4.17e+02 | -8.72e+02 | 0.000 | 0.9945 | 2.48e-05
2 | 2.48e-03 | 4.67e-13 | -4.17e+02 | -4.28e+02 | 1.592 | 0.9933 | 3.78e-06
3 | 1.60e-04 | 3.01e-14 | -4.17e+02 | -4.18e+02 | 2.782 | 0.9355 | 7.48e-06
4 | 2.74e-06 | 5.18e-16 | -4.17e+02 | -4.17e+02 | 3.917 | 0.9828 | 3.39e-08
5 | 5.89e-07 | 1.11e-16 | -4.17e+02 | -4.17e+02 | 4.537 | 0.7855 | 4.01e-07
6 | 1.56e-07 | 2.95e-17 | -4.17e+02 | -4.17e+02 | 5.096 | 0.7344 | 1.49e-07
7 | 1.33e-08 | 2.51e-18 | -4.17e+02 | -4.17e+02 | 6.006 | 0.9149 | 4.08e-09
8 | 2.26e-10 | 7.87e-20 | -4.17e+02 | -4.17e+02 | 7.576 | 0.9830 | 1.97e-11
9 | 1.67e-12 | 6.04e-20 | -4.17e+02 | -4.17e+02 | 9.283 | 0.9926 | 9.99e-14
10 | 8.54e-15 | 6.03e-20 | -4.17e+02 | -4.17e+02 | 11.430 | 0.9949 | 9.21e-16
11 | 4.74e-16 | 6.16e-20 | -4.17e+02 | -4.17e+02 | 13.736 | 0.9950 | 6.42e-18
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.30e-08 | -4.19e+02 | -6.69e+02 | 0.226 | 0.0000 | 1.00e+02
1 | 7.87e-01 | 2.70e-10 | -4.19e+02 | -8.73e+02 | 0.000 | 0.9883 | 1.13e-04
2 | 1.38e-02 | 4.71e-12 | -4.19e+02 | -4.42e+02 | 1.262 | 0.9825 | 2.83e-05
3 | 1.25e-03 | 4.27e-13 | -4.19e+02 | -4.21e+02 | 2.299 | 0.9093 | 3.26e-05
4 | 7.10e-05 | 2.43e-14 | -4.19e+02 | -4.20e+02 | 3.117 | 0.9431 | 1.15e-06
5 | 1.41e-05 | 4.81e-15 | -4.19e+02 | -4.19e+02 | 3.802 | 0.8021 | 2.15e-06
6 | 4.08e-06 | 1.40e-15 | -4.19e+02 | -4.19e+02 | 4.281 | 0.7095 | 9.76e-07
7 | 1.23e-06 | 4.22e-16 | -4.19e+02 | -4.19e+02 | 4.798 | 0.6982 | 3.50e-07
8 | 7.76e-09 | 2.64e-18 | -4.19e+02 | -4.19e+02 | 6.417 | 0.9937 | 4.40e-11
9 | 4.59e-11 | 5.23e-20 | -4.19e+02 | -4.19e+02 | 8.556 | 0.9941 | 9.31e-13
10 | 2.29e-13 | 5.48e-20 | -4.19e+02 | -4.19e+02 | 10.854 | 0.9950 | 4.83e-15
11 | 1.42e-15 | 4.80e-20 | -4.19e+02 | -4.19e+02 | 13.148 | 0.9950 | 2.43e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.92e-08 | -3.60e+02 | -6.10e+02 | 0.160 | 0.0000 | 1.00e+02
1 | 8.62e-01 | 2.46e-10 | -3.60e+02 | -8.13e+02 | 0.000 | 0.9872 | 1.36e-04
2 | 1.98e-02 | 5.66e-12 | -3.59e+02 | -3.86e+02 | 1.131 | 0.9770 | 4.99e-05
3 | 2.31e-03 | 6.60e-13 | -3.59e+02 | -3.62e+02 | 2.057 | 0.8833 | 6.35e-05
4 | 1.16e-04 | 3.32e-14 | -3.60e+02 | -3.60e+02 | 2.812 | 0.9497 | 1.34e-06
5 | 5.98e-06 | 1.71e-15 | -3.60e+02 | -3.60e+02 | 3.934 | 0.9487 | 2.45e-07
6 | 9.91e-07 | 2.83e-16 | -3.60e+02 | -3.60e+02 | 4.580 | 0.8341 | 1.97e-07
7 | 2.26e-08 | 6.44e-18 | -3.60e+02 | -3.60e+02 | 6.031 | 0.9772 | 8.18e-10
8 | 1.26e-10 | 6.10e-20 | -3.60e+02 | -3.60e+02 | 8.121 | 0.9944 | 1.74e-12
9 | 6.31e-13 | 4.64e-20 | -3.60e+02 | -3.60e+02 | 10.419 | 0.9950 | 1.13e-14
10 | 3.25e-15 | 4.68e-20 | -3.60e+02 | -3.60e+02 | 12.720 | 0.9950 | 5.68e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.00e-08 | -4.04e+02 | -6.54e+02 | 0.209 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 9.99e-11 | -4.04e+02 | -8.59e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 2.30e-03 | 6.82e-13 | -4.04e+02 | -4.14e+02 | 1.604 | 0.9932 | 3.86e-06
3 | 1.17e-04 | 3.46e-14 | -4.04e+02 | -4.04e+02 | 2.897 | 0.9492 | 4.36e-06
4 | 2.47e-05 | 7.32e-15 | -4.04e+02 | -4.04e+02 | 3.499 | 0.7886 | 3.94e-06
5 | 2.64e-06 | 7.83e-16 | -4.04e+02 | -4.04e+02 | 4.299 | 0.8931 | 2.47e-07
6 | 1.37e-07 | 4.08e-17 | -4.04e+02 | -4.04e+02 | 5.158 | 0.9479 | 9.19e-09
7 | 4.14e-09 | 1.22e-18 | -4.04e+02 | -4.04e+02 | 6.476 | 0.9698 | 4.25e-10
8 | 3.82e-11 | 5.88e-20 | -4.04e+02 | -4.04e+02 | 7.927 | 0.9908 | 1.91e-12
9 | 2.30e-13 | 4.70e-20 | -4.04e+02 | -4.04e+02 | 9.131 | 0.9940 | 2.89e-14
10 | 1.28e-15 | 4.64e-20 | -4.04e+02 | -4.04e+02 | 10.449 | 0.9950 | 1.24e-15
11 | 4.08e-16 | 3.94e-20 | -4.04e+02 | -4.04e+02 | 12.697 | 0.9950 | 5.95e-17
----------------------------------------------------------------------------------
optimization convergedggplot() +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = T0, linetype = "dotted", colour = "grey40") +
geom_line(data = pl_keep, aes(time, gap, group = unit),
colour = PAL[["donor"]], alpha = 0.7, linewidth = 0.5) +
geom_line(data = data.frame(time = 1:TT, gap = fit_good$gap),
aes(time, gap), colour = PAL[["treated"]], linewidth = 1.1) +
labs(x = "Period", y = "Treated - synthetic",
subtitle = "Good-support scenario; grey lines are fake-treated donors")
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(T0, color="grey", linestyle=":")
for p in pl_keep:
ax.plot(range(1, TT + 1), p["gap"], color="0.7", alpha=0.7, lw=0.5)
ax.plot(range(1, TT + 1), fit_good["gap"], color="#B22222", lw=1.5)
ax.set_xlabel("Period")
ax.set_ylabel("Treated - synthetic")
ax.set_title("Good-support scenario; grey lines are fake-treated donors")
plt.tight_layout()
plt.show()
* Treated unit's own gap (fit here so the block is self-contained).
import delimited "synthetic-control_files/data/sim_good.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly synth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20) ///
resultsperiod(1(1)35) keep("synthetic-control_files/stata/_tr", replace)
use "synthetic-control_files/stata/_tr", clear
gen treated_gap = _Y_treated - _Y_synthetic
rename _time time
keep time treated_gap
tempfile trg
save `trg'
* Placebo spaghetti (from placebo block), screened, with the treated line on top.
use "synthetic-control_files/stata/_placebo", clear
keep if rmspe <= 5 * 0.67
sort unit time // connect(ascending) draws one line per unit
merge m:1 time using `trg', nogenerate
twoway (line gap time, connect(ascending) lcolor(gs11) lwidth(vthin)) ///
(line treated_gap time, lcolor("178 34 34") lwidth(medthick)), ///
yline(0) xline(20, lpattern(dot)) legend(off) ///
xtitle("Period") ytitle("Treated - synthetic")
quietly graph export "sc_figs/fig-placebo-stata.svg", replace width(1600)(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
file /var/folders/2f/5dgnsms91455tbbpyrkgb15c9rnnyq/T//S_21672.000001 saved as .dta format
(0 observations deleted)
Result Number of obs
-----------------------------------------
Not matched 0
Matched 875
-----------------------------------------
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.32e-08 | -3.96e+02 | -6.46e+02 | 0.201 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.62e-11 | -3.96e+02 | -8.51e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.68e-03 | 3.31e-13 | -3.96e+02 | -4.05e+02 | 1.633 | 0.9950 | 2.07e-06
3 | 5.64e-05 | 1.11e-14 | -3.96e+02 | -3.96e+02 | 3.106 | 0.9664 | 1.74e-06
4 | 6.30e-06 | 1.24e-15 | -3.96e+02 | -3.96e+02 | 4.026 | 0.8884 | 6.53e-07
5 | 1.83e-06 | 3.61e-16 | -3.96e+02 | -3.96e+02 | 4.560 | 0.7090 | 5.53e-07
6 | 3.44e-07 | 6.78e-17 | -3.96e+02 | -3.96e+02 | 5.235 | 0.8123 | 6.59e-08
7 | 7.80e-09 | 1.57e-18 | -3.96e+02 | -3.96e+02 | 6.565 | 0.9773 | 1.97e-10
8 | 1.62e-10 | 5.22e-20 | -3.96e+02 | -3.96e+02 | 8.034 | 0.9792 | 7.76e-12
9 | 9.57e-13 | 6.42e-20 | -3.96e+02 | -3.96e+02 | 10.040 | 0.9941 | 2.11e-14
10 | 4.90e-15 | 6.69e-20 | -3.96e+02 | -3.96e+02 | 12.338 | 0.9950 | 1.50e-16
----------------------------------------------------------------------------------
optimization converged
The permutation p-value follows directly. We use the mean post-treatment gap as the statistic and compare in the direction of the hypothesized negative effect.
stat <- function(g, T0, TT) mean(g[(T0 + 1):TT])
placebo_stats <- pl_keep |>
group_by(unit) |>
summarise(S = stat(gap, T0, TT), .groups = "drop")
S_treated <- stat(fit_good$gap, T0, TT)
p_val <- (1 + sum(placebo_stats$S <= S_treated)) / (1 + nrow(placebo_stats))
cat("Treated statistic:", round(S_treated, 2), "\n")Treated statistic: -9.73 cat("Placebo units :", nrow(placebo_stats), "\n")Placebo units : 25 cat("One-sided p-value:", round(p_val, 4), "\n")One-sided p-value: 0.0385 def stat(gap):
return gap[T0:].mean()
placebo_S = [stat(p["gap"]) for p in pl_keep]
S_treated = stat(fit_good["gap"])
p_val = (1 + sum(s <= S_treated for s in placebo_S)) / (1 + len(placebo_S))
print(f"Treated statistic: {round(S_treated, 2)}")Treated statistic: -9.73print(f"Placebo units : {len(placebo_S)}")Placebo units : 25print(f"One-sided p-value: {round(p_val, 4)}")One-sided p-value: 0.0385* Treated statistic (mean post-period gap) from the good-support fit.
import delimited "synthetic-control_files/data/sim_good.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly synth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20) ///
resultsperiod(1(1)35) keep("synthetic-control_files/stata/_tr", replace)
use "synthetic-control_files/stata/_tr", clear
gen gap = _Y_treated - _Y_synthetic
quietly summarize gap if _time > 20
scalar S_tr = r(mean)
* Placebo statistics, screened, one per unit.
use "synthetic-control_files/stata/_placebo", clear
keep if rmspe <= 5 * 0.67
collapse (mean) S = gap if time > 20, by(unit)
count if S <= S_tr
scalar p = (1 + r(N)) / (1 + _N)
display "Treated statistic: " %5.2f S_tr
display "Placebo units : " _N
display "One-sided p-value: " %6.4f p(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(0 observations deleted)
0
Treated statistic: -9.73
Placebo units : 25
One-sided p-value: 0.0385
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.32e-08 | -3.96e+02 | -6.46e+02 | 0.201 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.62e-11 | -3.96e+02 | -8.51e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.68e-03 | 3.31e-13 | -3.96e+02 | -4.05e+02 | 1.633 | 0.9950 | 2.07e-06
3 | 5.64e-05 | 1.11e-14 | -3.96e+02 | -3.96e+02 | 3.106 | 0.9664 | 1.74e-06
4 | 6.30e-06 | 1.24e-15 | -3.96e+02 | -3.96e+02 | 4.026 | 0.8884 | 6.53e-07
5 | 1.83e-06 | 3.61e-16 | -3.96e+02 | -3.96e+02 | 4.560 | 0.7090 | 5.53e-07
6 | 3.44e-07 | 6.78e-17 | -3.96e+02 | -3.96e+02 | 5.235 | 0.8123 | 6.59e-08
7 | 7.80e-09 | 1.57e-18 | -3.96e+02 | -3.96e+02 | 6.565 | 0.9773 | 1.97e-10
8 | 1.62e-10 | 5.22e-20 | -3.96e+02 | -3.96e+02 | 8.034 | 0.9792 | 7.76e-12
9 | 9.57e-13 | 6.42e-20 | -3.96e+02 | -3.96e+02 | 10.040 | 0.9941 | 2.11e-14
10 | 4.90e-15 | 6.69e-20 | -3.96e+02 | -3.96e+02 | 12.338 | 0.9950 | 1.50e-16
----------------------------------------------------------------------------------
optimization convergedLeave-One-Out
A conclusion that rests on a single donor is fragile. We refit while dropping each donor that carries meaningful weight.
top_donors <- names(sort(fit_good$w, decreasing = TRUE))[1:4]
top_ids <- as.integer(sub("^u", "", top_donors))
loo <- lapply(top_ids, function(drop_id) {
donors <- setdiff(2:(J + 1), drop_id)
f <- fit_scm(panel_good, 1, donors, T0, TT)
data.frame(time = 1:TT, gap = f$gap, dropped = paste0("without u", drop_id))
}) |> bind_rows()
ggplot() +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = T0, linetype = "dotted", colour = "grey40") +
geom_line(data = loo, aes(time, gap, group = dropped, colour = dropped),
linewidth = 0.6) +
geom_line(data = data.frame(time = 1:TT, gap = fit_good$gap),
aes(time, gap), colour = PAL[["treated"]], linewidth = 1.2) +
labs(x = "Period", y = "Treated - synthetic", colour = "Leave-one-out",
subtitle = "Thick red line is the full-donor-pool estimate")
top = fit_good["w"].sort_values(ascending=False).index[:4].tolist()
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(T0, color="grey", linestyle=":")
for d in top:
others = [u for u in donors if u != d]
f = fit_scm(panel_good, "u1", others, T0, TT)
ax.plot(range(1, TT + 1), f["gap"], lw=0.8, label=f"without {d}")
ax.plot(range(1, TT + 1), fit_good["gap"], color="#B22222", lw=2)
ax.legend(frameon=False, title="Leave-one-out")
ax.set_xlabel("Period")
ax.set_ylabel("Treated - synthetic")
plt.tight_layout()
plt.show()
* Identify the top-weight donors from a full fit, then refit dropping each.
import delimited "synthetic-control_files/data/sim_good.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly synth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20)
matrix W = e(W_weights)
* (first-draft) drop a fixed set of high-weight donor ids and overlay refits
graph drop _all
local k = 0
foreach d in 2 3 4 5 {
local ++k
import delimited "synthetic-control_files/data/sim_good.csv", clear
drop if id == `d'
tsset id time
quietly synth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20) ///
resultsperiod(1(1)35) keep("synthetic-control_files/stata/_loo_`d'", replace)
use "synthetic-control_files/stata/_loo_`d'", clear
gen gap = _Y_treated - _Y_synthetic
twoway line gap _time, name(loo`k', replace) nodraw ///
title("without u`d'") yline(0) xline(20, lpattern(dot)) legend(off)
}
graph combine loo1 loo2 loo3 loo4, ycommon
quietly graph export "sc_figs/fig-loo-stata.svg", replace width(1600)(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
7. use "synthetic-control_files/stata/_loo_`d'", clear
8. gen gap = _Y_treated - _Y_synthetic
9. twoway line gap _time, name(loo`k', replace) nodraw ///
> title("without u`d'") yline(0) xline(20, lpattern(dot)) legend(off)
10. }
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
(35 observations deleted)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
(35 observations deleted)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
(35 observations deleted)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(encoding automatically selected: UTF-8)
(4 vars, 910 obs)
(35 observations deleted)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.32e-08 | -3.96e+02 | -6.46e+02 | 0.201 | 0.0000 | 1.00e+02
1 | 3.36e-01 | 6.62e-11 | -3.96e+02 | -8.51e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.68e-03 | 3.31e-13 | -3.96e+02 | -4.05e+02 | 1.633 | 0.9950 | 2.07e-06
3 | 5.64e-05 | 1.11e-14 | -3.96e+02 | -3.96e+02 | 3.106 | 0.9664 | 1.74e-06
4 | 6.30e-06 | 1.24e-15 | -3.96e+02 | -3.96e+02 | 4.026 | 0.8884 | 6.53e-07
5 | 1.83e-06 | 3.61e-16 | -3.96e+02 | -3.96e+02 | 4.560 | 0.7090 | 5.53e-07
6 | 3.44e-07 | 6.78e-17 | -3.96e+02 | -3.96e+02 | 5.235 | 0.8123 | 6.59e-08
7 | 7.80e-09 | 1.57e-18 | -3.96e+02 | -3.96e+02 | 6.565 | 0.9773 | 1.97e-10
8 | 1.62e-10 | 5.22e-20 | -3.96e+02 | -3.96e+02 | 8.034 | 0.9792 | 7.76e-12
9 | 9.57e-13 | 6.42e-20 | -3.96e+02 | -3.96e+02 | 10.040 | 0.9941 | 2.11e-14
10 | 4.90e-15 | 6.69e-20 | -3.96e+02 | -3.96e+02 | 12.338 | 0.9950 | 1.50e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.59e+01 | 1.48e-08 | -3.84e+02 | -6.24e+02 | 0.206 | 0.0000 | 1.00e+02
1 | 3.29e-01 | 7.38e-11 | -3.84e+02 | -8.21e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.65e-03 | 3.69e-13 | -3.84e+02 | -3.93e+02 | 1.637 | 0.9950 | 2.07e-06
3 | 5.31e-05 | 1.19e-14 | -3.84e+02 | -3.84e+02 | 3.128 | 0.9678 | 1.60e-06
4 | 5.37e-06 | 1.20e-15 | -3.84e+02 | -3.84e+02 | 4.087 | 0.8988 | 5.16e-07
5 | 1.63e-06 | 3.64e-16 | -3.84e+02 | -3.84e+02 | 4.604 | 0.6972 | 5.26e-07
6 | 2.97e-07 | 6.65e-17 | -3.84e+02 | -3.84e+02 | 5.292 | 0.8176 | 5.68e-08
7 | 1.30e-08 | 2.93e-18 | -3.84e+02 | -3.84e+02 | 6.475 | 0.9560 | 6.60e-10
8 | 2.52e-10 | 7.59e-20 | -3.84e+02 | -3.84e+02 | 7.990 | 0.9807 | 8.35e-12
9 | 1.52e-12 | 6.33e-20 | -3.84e+02 | -3.84e+02 | 10.014 | 0.9940 | 2.47e-14
10 | 7.60e-15 | 9.35e-20 | -3.84e+02 | -3.84e+02 | 12.311 | 0.9950 | 1.61e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.59e+01 | 1.43e-08 | -3.95e+02 | -6.35e+02 | 0.218 | 0.0000 | 1.00e+02
1 | 3.29e-01 | 7.13e-11 | -3.95e+02 | -8.32e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.65e-03 | 3.56e-13 | -3.95e+02 | -4.04e+02 | 1.649 | 0.9950 | 2.07e-06
3 | 6.58e-05 | 1.42e-14 | -3.95e+02 | -3.95e+02 | 3.046 | 0.9600 | 2.47e-06
4 | 4.36e-06 | 9.44e-16 | -3.95e+02 | -3.95e+02 | 4.143 | 0.9337 | 2.73e-07
5 | 1.11e-06 | 2.40e-16 | -3.95e+02 | -3.95e+02 | 4.710 | 0.7461 | 3.31e-07
6 | 1.56e-07 | 3.36e-17 | -3.95e+02 | -3.95e+02 | 5.482 | 0.8596 | 2.69e-08
7 | 1.58e-08 | 3.40e-18 | -3.95e+02 | -3.95e+02 | 6.447 | 0.8987 | 2.35e-09
8 | 1.38e-10 | 5.27e-20 | -3.95e+02 | -3.95e+02 | 8.258 | 0.9912 | 1.87e-12
9 | 1.04e-12 | 4.96e-20 | -3.95e+02 | -3.95e+02 | 9.839 | 0.9924 | 2.16e-14
10 | 5.24e-15 | 5.15e-20 | -3.95e+02 | -3.95e+02 | 12.007 | 0.9950 | 2.47e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.59e+01 | 1.42e-08 | -3.93e+02 | -6.33e+02 | 0.215 | 0.0000 | 1.00e+02
1 | 3.29e-01 | 7.09e-11 | -3.93e+02 | -8.30e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.65e-03 | 3.54e-13 | -3.93e+02 | -4.02e+02 | 1.646 | 0.9950 | 2.07e-06
3 | 6.18e-05 | 1.33e-14 | -3.93e+02 | -3.93e+02 | 3.071 | 0.9625 | 2.18e-06
4 | 5.48e-06 | 1.18e-15 | -3.93e+02 | -3.93e+02 | 4.076 | 0.9113 | 4.61e-07
5 | 1.67e-06 | 3.60e-16 | -3.93e+02 | -3.93e+02 | 4.589 | 0.6953 | 5.59e-07
6 | 3.03e-07 | 6.53e-17 | -3.93e+02 | -3.93e+02 | 5.280 | 0.8184 | 5.96e-08
7 | 5.48e-09 | 1.17e-18 | -3.93e+02 | -3.93e+02 | 6.700 | 0.9819 | 1.17e-10
8 | 6.23e-11 | 8.24e-20 | -3.93e+02 | -3.93e+02 | 8.203 | 0.9886 | 1.75e-12
9 | 4.22e-13 | 6.32e-20 | -3.93e+02 | -3.93e+02 | 9.961 | 0.9932 | 1.95e-14
10 | 2.16e-15 | 4.65e-20 | -3.93e+02 | -3.93e+02 | 12.259 | 0.9950 | 1.86e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.59e+01 | 1.47e-08 | -3.86e+02 | -6.26e+02 | 0.208 | 0.0000 | 1.00e+02
1 | 3.29e-01 | 7.34e-11 | -3.86e+02 | -8.23e+02 | 0.000 | 0.9950 | 2.07e-05
2 | 1.65e-03 | 3.67e-13 | -3.86e+02 | -3.95e+02 | 1.639 | 0.9950 | 2.07e-06
3 | 5.40e-05 | 1.20e-14 | -3.86e+02 | -3.87e+02 | 3.122 | 0.9672 | 1.66e-06
4 | 5.17e-06 | 1.15e-15 | -3.86e+02 | -3.86e+02 | 4.101 | 0.9042 | 4.70e-07
5 | 1.60e-06 | 3.57e-16 | -3.86e+02 | -3.86e+02 | 4.608 | 0.6908 | 5.36e-07
6 | 2.89e-07 | 6.45e-17 | -3.86e+02 | -3.86e+02 | 5.301 | 0.8193 | 5.55e-08
7 | 1.05e-08 | 2.35e-18 | -3.86e+02 | -3.86e+02 | 6.531 | 0.9635 | 4.46e-10
8 | 1.95e-10 | 7.29e-20 | -3.86e+02 | -3.86e+02 | 8.040 | 0.9815 | 6.77e-12
9 | 1.18e-12 | 7.00e-20 | -3.86e+02 | -3.86e+02 | 10.044 | 0.9940 | 2.23e-14
10 | 5.94e-15 | 5.59e-20 | -3.86e+02 | -3.86e+02 | 12.342 | 0.9950 | 1.51e-16
----------------------------------------------------------------------------------
optimization converged
What the Diagnostics Caught
| Diagnostic | Good support | Poor support |
|---|---|---|
| Pre-treatment RMSPE | Small relative to outcome scale | An order of magnitude larger, an immediate warning |
| Path plot | Lines overlap before treatment | Lines never align, visible without the truth |
| Gap plot | Flat at zero pre-treatment | Large pre-treatment swings |
| Estimated vs. true ATT | Close | Attenuated to roughly a sixth of the truth |
The point is not that SCM worked in one case and failed in the other. It is that the failure was detectable from the pre-treatment period alone, using only information available to any reader of the figures. That is the property that makes the method defensible, and it is why the diagnostics are not optional decoration.
Classical SCM with Synth: Proposition 99
Now the real data. The tidysynth package ships the Abadie, Diamond, and Hainmueller panel: 39 states, 1970 to 2000, annual per-capita cigarette sales (Abadie, Diamond, and Hainmueller 2010).1
Data Requirements
The Synth package implements the classical estimator used in the foundational applications (Hainmueller and Diamond 2026). It expects long panel data, one row per unit-period.
| Variable | Role |
|---|---|
| Unit identifier | Numeric ID used by the optimizer |
| Unit name | Readable label for tables and plots |
| Time identifier | A consistently coded period variable |
| Outcome | The variable whose untreated path you want to estimate |
| Predictors | Pre-treatment covariates and outcome summaries |
| Treatment date | A design choice, not a column Synth infers |
Donors must be untreated throughout the analysis window and plausibly unaffected by the intervention. The published donor pool already excludes states that adopted their own large tobacco-control programs.
# Local cache of tidysynth::smoking (see synthetic-control_files/data/PROVENANCE.md).
smoking <- read.csv("synthetic-control_files/data/smoking.csv")
smoking_df <- smoking |>
mutate(state_id = as.numeric(factor(state))) |>
as.data.frame()
ca_id <- unique(smoking_df$state_id[smoking_df$state == "California"])
donor_ids <- sort(setdiff(unique(smoking_df$state_id), ca_id))
cat("States:", length(unique(smoking_df$state)),
"| California id:", ca_id, "| donors:", length(donor_ids),
"| years:", paste(range(smoking_df$year), collapse = "-"), "\n")States: 39 | California id: 3 | donors: 38 | years: 1970-2000 import pandas as pd
smoking = pd.read_csv("synthetic-control_files/data/smoking.csv")
# as.numeric(factor(state)): 1-based alphabetical codes (California -> 3)
smoking["state_id"] = pd.Categorical(smoking["state"]).codes + 1
ca_id = int(smoking.loc[smoking.state == "California", "state_id"].iloc[0])
donor_ids = sorted(set(smoking.state_id) - {ca_id})
print(f"States: {smoking.state.nunique()} | California id: {ca_id} | "
f"donors: {len(donor_ids)} | "
f"years: {smoking.year.min()}-{smoking.year.max()}")States: 39 | California id: 3 | donors: 38 | years: 1970-2000import delimited "synthetic-control_files/data/smoking.csv", clear
* R writes missing as "NA", so columns with gaps import as strings.
destring lnincome beer age15to24, replace force
encode state, gen(state_id) // alphabetical, California -> 3
quietly summarize state_id if state == "California"
local ca = r(mean)
quietly levelsof state_id, local(ids)
local nstates : word count `ids'
local ndonors = `nstates' - 1
quietly summarize year
display "States: " `nstates' " | California id: " `ca' ///
" | donors: " `ndonors' " | years: " r(min) "-" r(max)(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
States: 39 | California id: 3 | donors: 38 | years: 1970-2000Fitting the Synthetic California
The specification follows the published analysis: income, retail price, and the share aged 15 to 24 averaged over 1980 to 1988, beer consumption over 1984 to 1988 (the years it is observed), and cigarette sales in 1975, 1980, and 1988.
prep_ca <- dataprep(
foo = smoking_df,
predictors = c("lnincome", "retprice", "age15to24"),
predictors.op = "mean",
dependent = "cigsale",
unit.variable = "state_id",
unit.names.variable = "state",
time.variable = "year",
special.predictors = list(
list("beer", 1984:1988, "mean"),
list("cigsale", 1975, "mean"),
list("cigsale", 1980, "mean"),
list("cigsale", 1988, "mean")
),
treatment.identifier = ca_id,
controls.identifier = donor_ids,
time.predictors.prior = 1980:1988,
time.optimize.ssr = 1970:1988,
time.plot = 1970:2000
)
invisible(capture.output(fit_ca <- synth(prep_ca)))
tabs_ca <- synth.tab(dataprep.res = prep_ca, synth.res = fit_ca)from pysyncon import Dataprep, Synth
donors = [s for s in smoking.state.unique() if s != "California"]
dp_ca = Dataprep(
foo=smoking, dependent="cigsale",
unit_variable="state", time_variable="year",
treatment_identifier="California", controls_identifier=donors,
predictors=["lnincome", "retprice", "age15to24"], predictors_op="mean",
time_predictors_prior=range(1980, 1989),
special_predictors=[
("beer", range(1984, 1989), "mean"),
("cigsale", [1975], "mean"),
("cigsale", [1980], "mean"),
("cigsale", [1988], "mean"),
],
time_optimize_ssr=range(1970, 1989),
)
synth_ca = Synth()
synth_ca.fit(dp_ca) # silent fit; weights and gaps read out belowimport delimited "synthetic-control_files/data/smoking.csv", clear
destring lnincome beer age15to24, replace force
encode state, gen(state_id)
tsset state_id year
quietly summarize state_id if state == "California"
* synth prints unit weights, predictor balance, and RMSPE in one table.
synth cigsale lnincome(1980(1)1988) retprice(1980(1)1988) age15to24(1980(1)1988) ///
beer(1984(1)1988) cigsale(1975) cigsale(1980) cigsale(1988), ///
trunit(`r(mean)') trperiod(1989) mspeperiod(1970(1)1988) ///
resultsperiod(1970(1)2000)(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
Panel variable: state_id (strongly balanced)
Time variable: year, 1970 to 2000
Delta: 1 unit
--------------------------------------------------------------------------------------------------------------------------------------------
Synthetic Control Method for Comparative Case Studies
--------------------------------------------------------------------------------------------------------------------------------------------
First Step: Data Setup
--------------------------------------------------------------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------------------------------------------------------
Data Setup successful
--------------------------------------------------------------------------------------------------------------------------------------------
Treated Unit: California
Control Units: Alabama, Arkansas, Colorado, Connecticut, Delaware, Georgia, Idaho, Illinois, Indiana, Iowa, Kansas, Kentucky,
Louisiana, Maine, Minnesota, Mississippi, Missouri, Montana, Nebraska, Nevada, New Hampshire, New Mexico,
North Carolina, North Dakota, Ohio, Oklahoma, Pennsylvania, Rhode Island, South Carolina, South Dakota,
Tennessee, Texas, Utah, Vermont, Virginia, West Virginia, Wisconsin, Wyoming
--------------------------------------------------------------------------------------------------------------------------------------------
Dependent Variable: cigsale
MSPE minimized for periods: 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988
Results obtained for periods: 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991
1992 1993 1994 1995 1996 1997 1998 1999 2000
--------------------------------------------------------------------------------------------------------------------------------------------
Predictors: lnincome(1980(1)1988) retprice(1980(1)1988) age15to24(1980(1)1988) beer(1984(1)1988) cigsale(1975)
cigsale(1980) cigsale(1988)
--------------------------------------------------------------------------------------------------------------------------------------------
Unless period is specified
predictors are averaged over: 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988
--------------------------------------------------------------------------------------------------------------------------------------------
Second Step: Run Optimization
--------------------------------------------------------------------------------------------------------------------------------------------
--------------------------------------------------------------------------------------------------------------------------------------------
Optimization done
--------------------------------------------------------------------------------------------------------------------------------------------
Third Step: Obtain Results
--------------------------------------------------------------------------------------------------------------------------------------------
Loss: Root Mean Squared Prediction Error
--------------------
RMSPE | 2.03785
--------------------
--------------------------------------------------------------------------------------------------------------------------------------------
Unit Weights:
----------------------------
Co_No | Unit_Weight
---------------+------------
Alabama | 0
Arkansas | 0
Colorado | .356
Connecticut | .083
Delaware | 0
Georgia | 0
Idaho | 0
Illinois | 0
Indiana | 0
Iowa | 0
Kansas | 0
Kentucky | 0
Louisiana | 0
Maine | 0
Minnesota | 0
Mississippi | 0
Missouri | 0
Montana | 0
Nebraska | 0
Nevada | .217
New Hampshire | 0
New Mexico | 0
North Carolina | 0
North Dakota | 0
Ohio | 0
Oklahoma | 0
Pennsylvania | 0
Rhode Island | 0
South Carolina | 0
South Dakota | 0
Tennessee | 0
Texas | 0
Utah | .344
Vermont | 0
Virginia | 0
West Virginia | 0
Wisconsin | 0
Wyoming | 0
----------------------------
--------------------------------------------------------------------------------------------------------------------------------------------
Predictor Balance:
-----------------------------------------------------
| Treated Synthetic
-------------------------------+---------------------
lnincome(1980(1)1988) | 10.07656 9.904773
retprice(1980(1)1988) | 89.42222 89.00128
age15to24(1980(1)1988) | .1735324 .175389
beer(1984(1)1988) | 24.28 23.26454
cigsale(1975) | 127.1 126.3862
cigsale(1980) | 120.2 120.7221
cigsale(1988) | 90.1 92.0883
-----------------------------------------------------
--------------------------------------------------------------------------------------------------------------------------------------------
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 8.29e+01 | 6.80e-06 | -1.36e+01 | -3.94e+02 | 0.000 | 0.0000 | 1.00e+02
1 | 5.05e-01 | 4.14e-08 | -1.36e+01 | -7.03e+02 | 0.000 | 0.9939 | 3.06e-05
2 | 2.87e-03 | 2.35e-10 | -1.36e+01 | -2.90e+01 | 0.000 | 0.9943 | 2.72e-06
3 | 1.63e-04 | 1.33e-11 | -1.36e+01 | -1.45e+01 | 1.220 | 0.9433 | 5.50e-06
4 | 1.60e-05 | 1.32e-12 | -1.36e+01 | -1.38e+01 | 2.024 | 0.9013 | 9.47e-07
5 | 8.35e-06 | 6.85e-13 | -1.36e+01 | -1.37e+01 | 2.305 | 0.4792 | 4.50e-06
6 | 2.37e-06 | 1.94e-13 | -1.36e+01 | -1.37e+01 | 2.833 | 0.7163 | 6.69e-07
7 | 7.27e-07 | 5.96e-14 | -1.36e+01 | -1.37e+01 | 3.334 | 0.6933 | 2.33e-07
8 | 1.61e-07 | 1.32e-14 | -1.36e+01 | -1.36e+01 | 3.924 | 0.7790 | 3.75e-08
9 | 3.88e-09 | 3.18e-16 | -1.36e+01 | -1.36e+01 | 5.302 | 0.9758 | 1.11e-10
10 | 1.99e-11 | 2.13e-18 | -1.36e+01 | -1.36e+01 | 7.563 | 0.9949 | 2.09e-13
11 | 9.92e-14 | 1.22e-18 | -1.36e+01 | -1.36e+01 | 9.861 | 0.9950 | 1.09e-15
12 | 8.56e-16 | 1.74e-18 | -1.36e+01 | -1.36e+01 | 12.158 | 0.9950 | 5.49e-18
----------------------------------------------------------------------------------
optimization convergedRead the weights before the gap plot.
w_ca <- tabs_ca$tab.w
w_ca[w_ca$w.weights > 0.001, ] |>
arrange(desc(w.weights)) w.weights unit.names unit.numbers
34 0.343 Utah 34
21 0.236 Nevada 21
19 0.182 Montana 19
4 0.175 Colorado 4
5 0.062 Connecticut 5w_py = synth_ca.weights().sort_values(ascending=False)
print(w_py[w_py > 0.001].round(4).to_string())Colorado 0.625
Connecticut 0.273
Utah 0.101import delimited "synthetic-control_files/data/smoking.csv", clear
destring lnincome beer age15to24, replace force
encode state, gen(state_id)
tsset state_id year
quietly summarize state_id if state == "California"
quietly synth cigsale lnincome(1980(1)1988) retprice(1980(1)1988) ///
age15to24(1980(1)1988) beer(1984(1)1988) cigsale(1975) cigsale(1980) ///
cigsale(1988), trunit(`r(mean)') trperiod(1989) mspeperiod(1970(1)1988)
matrix W = e(W_weights) // column 1 = weight, column 2 = unit id
matrix list W, format(%6.4f)(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
Panel variable: state_id (strongly balanced)
Time variable: year, 1970 to 2000
Delta: 1 unit
W[38,2]
_Co_Number _W_Weight
Alabama 1.0000 0.0000
Arkansas 2.0000 0.0000
Colorado 4.0000 0.3560
Connecticut 5.0000 0.0830
Delaware 6.0000 0.0000
Georgia 7.0000 0.0000
Idaho 8.0000 0.0000
Illinois 9.0000 0.0000
Indiana 10.0000 0.0000
Iowa 11.0000 0.0000
Kansas 12.0000 0.0000
Kentucky 13.0000 0.0000
Louisiana 14.0000 0.0000
Maine 15.0000 0.0000
Minnesota 16.0000 0.0000
Mississippi 17.0000 0.0000
Missouri 18.0000 0.0000
Montana 19.0000 0.0000
Nebraska 20.0000 0.0000
Nevada 21.0000 0.2170
NewHampshire 22.0000 0.0000
NewMexico 23.0000 0.0000
NorthCarol~a 24.0000 0.0000
NorthDakota 25.0000 0.0000
Ohio 26.0000 0.0000
Oklahoma 27.0000 0.0000
Pennsylvania 28.0000 0.0000
RhodeIsland 29.0000 0.0000
SouthCarol~a 30.0000 0.0000
SouthDakota 31.0000 0.0000
Tennessee 32.0000 0.0000
Texas 33.0000 0.0000
Utah 34.0000 0.3440
Vermont 35.0000 0.0000
Virginia 36.0000 0.0000
WestVirginia 37.0000 0.0000
Wisconsin 38.0000 0.0000
Wyoming 39.0000 0.0000
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 8.29e+01 | 6.80e-06 | -1.36e+01 | -3.94e+02 | 0.000 | 0.0000 | 1.00e+02
1 | 5.05e-01 | 4.14e-08 | -1.36e+01 | -7.03e+02 | 0.000 | 0.9939 | 3.06e-05
2 | 2.87e-03 | 2.35e-10 | -1.36e+01 | -2.90e+01 | 0.000 | 0.9943 | 2.72e-06
3 | 1.63e-04 | 1.33e-11 | -1.36e+01 | -1.45e+01 | 1.220 | 0.9433 | 5.50e-06
4 | 1.60e-05 | 1.32e-12 | -1.36e+01 | -1.38e+01 | 2.024 | 0.9013 | 9.47e-07
5 | 8.35e-06 | 6.85e-13 | -1.36e+01 | -1.37e+01 | 2.305 | 0.4792 | 4.50e-06
6 | 2.37e-06 | 1.94e-13 | -1.36e+01 | -1.37e+01 | 2.833 | 0.7163 | 6.69e-07
7 | 7.27e-07 | 5.96e-14 | -1.36e+01 | -1.37e+01 | 3.334 | 0.6933 | 2.33e-07
8 | 1.61e-07 | 1.32e-14 | -1.36e+01 | -1.36e+01 | 3.924 | 0.7790 | 3.75e-08
9 | 3.88e-09 | 3.18e-16 | -1.36e+01 | -1.36e+01 | 5.302 | 0.9758 | 1.11e-10
10 | 1.99e-11 | 2.13e-18 | -1.36e+01 | -1.36e+01 | 7.563 | 0.9949 | 2.09e-13
11 | 9.92e-14 | 1.22e-18 | -1.36e+01 | -1.36e+01 | 9.861 | 0.9950 | 1.09e-15
12 | 8.56e-16 | 1.74e-18 | -1.36e+01 | -1.36e+01 | 12.158 | 0.9950 | 5.49e-18
----------------------------------------------------------------------------------
optimization converged
The three implementations agree on the effect, not the weights
All three fit the same specification, and all three imply a similar Proposition 99 effect (mean post-1988 gap near \(-18\) to \(-20\) packs, and roughly \(-25\) by 2000). But the donor weights differ: R leans on Utah and Nevada, pysyncon on Colorado and Connecticut, Stata’s synth on a Colorado/Utah/Nevada blend. This is not a bug; it is the \(V/W\) non-uniqueness discussed above (Klößner et al. 2018). The estimand is robust across solvers; the exact synthetic unit is not. Report which implementation and settings produced your weights.
Synthetic California is a blend of five states, with Utah and Nevada carrying the most weight. Whether that is substantively sensible is a question for a reader who knows something about tobacco markets, which is exactly the sort of judgment the method makes possible.
tabs_ca$tab.pred Treated Synthetic Sample Mean
lnincome 10.077 9.859 9.829
retprice 89.422 89.318 87.266
age15to24 0.174 0.174 0.173
special.beer.1984.1988 24.280 24.095 23.655
special.cigsale.1975 127.100 126.897 136.932
special.cigsale.1980 120.200 120.248 138.089
special.cigsale.1988 90.100 91.432 113.824# pysyncon prints predictor balance (treated vs synthetic vs donor sample) via summary().
synth_ca.summary() V treated synthetic sample mean
lnincome 0.077 10.077 10.009 9.829
retprice 0.120 89.422 89.016 87.266
age15to24 0.218 0.174 0.172 0.173
special.1.beer 0.017 24.280 22.694 23.655
special.2.cigsale 0.150 127.100 119.724 136.932
special.3.cigsale 0.246 120.200 121.756 138.089
special.4.cigsale 0.172 90.100 93.382 113.824import delimited "synthetic-control_files/data/smoking.csv", clear
destring lnincome beer age15to24, replace force
encode state, gen(state_id)
tsset state_id year
quietly summarize state_id if state == "California"
quietly synth cigsale lnincome(1980(1)1988) retprice(1980(1)1988) ///
age15to24(1980(1)1988) beer(1984(1)1988) cigsale(1975) cigsale(1980) ///
cigsale(1988), trunit(`r(mean)') trperiod(1989) mspeperiod(1970(1)1988)
matrix list e(X_balance) // treated, synthetic, and donor-average predictors(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
Panel variable: state_id (strongly balanced)
Time variable: year, 1970 to 2000
Delta: 1 unit
e(X_balance)[7,2]
Treated Synthetic
lninco~1988) 10.076559 9.9047732
retpri~1988) 89.422223 89.001278
age15t~1988) .17353238 .17538896
beer(1~1988) 24.28 23.26454
cigsale(1975) 127.1 126.3862
cigsale(1980) 120.2 120.7221
cigsale(1988) 90.099998 92.088298
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 8.29e+01 | 6.80e-06 | -1.36e+01 | -3.94e+02 | 0.000 | 0.0000 | 1.00e+02
1 | 5.05e-01 | 4.14e-08 | -1.36e+01 | -7.03e+02 | 0.000 | 0.9939 | 3.06e-05
2 | 2.87e-03 | 2.35e-10 | -1.36e+01 | -2.90e+01 | 0.000 | 0.9943 | 2.72e-06
3 | 1.63e-04 | 1.33e-11 | -1.36e+01 | -1.45e+01 | 1.220 | 0.9433 | 5.50e-06
4 | 1.60e-05 | 1.32e-12 | -1.36e+01 | -1.38e+01 | 2.024 | 0.9013 | 9.47e-07
5 | 8.35e-06 | 6.85e-13 | -1.36e+01 | -1.37e+01 | 2.305 | 0.4792 | 4.50e-06
6 | 2.37e-06 | 1.94e-13 | -1.36e+01 | -1.37e+01 | 2.833 | 0.7163 | 6.69e-07
7 | 7.27e-07 | 5.96e-14 | -1.36e+01 | -1.37e+01 | 3.334 | 0.6933 | 2.33e-07
8 | 1.61e-07 | 1.32e-14 | -1.36e+01 | -1.36e+01 | 3.924 | 0.7790 | 3.75e-08
9 | 3.88e-09 | 3.18e-16 | -1.36e+01 | -1.36e+01 | 5.302 | 0.9758 | 1.11e-10
10 | 1.99e-11 | 2.13e-18 | -1.36e+01 | -1.36e+01 | 7.563 | 0.9949 | 2.09e-13
11 | 9.92e-14 | 1.22e-18 | -1.36e+01 | -1.36e+01 | 9.861 | 0.9950 | 1.09e-15
12 | 8.56e-16 | 1.74e-18 | -1.36e+01 | -1.36e+01 | 12.158 | 0.9950 | 5.49e-18
----------------------------------------------------------------------------------
optimization convergedPredictor balance improves markedly on the raw donor average, most visibly for the pre-treatment sales figures.
years_ca <- as.numeric(rownames(prep_ca$Y1plot))
synth_ca <- as.vector(prep_ca$Y0plot %*% fit_ca$solution.w)
obs_ca <- as.vector(prep_ca$Y1plot)
gap_ca <- obs_ca - synth_ca
data.frame(year = years_ca, California = obs_ca,
`Synthetic California` = synth_ca, check.names = FALSE) |>
pivot_longer(-year, names_to = "series", values_to = "cigsale") |>
ggplot(aes(year, cigsale, colour = series, linetype = series)) +
geom_vline(xintercept = 1988, linetype = "dotted", colour = "grey40") +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = setNames(c(PAL[["treated"]], PAL[["synth"]]),
c("California", "Synthetic California"))) +
labs(x = NULL, y = "Cigarette sales (packs per capita)",
colour = NULL, linetype = NULL,
subtitle = "Dotted line marks Proposition 99 (1988)")
wide_ca = smoking.pivot(index="year", columns="state", values="cigsale").sort_index()
synth_path = wide_ca[donors].values @ np.asarray(synth_ca.W)
obs_path = wide_ca["California"].values
years = wide_ca.index.values
gap_ca = obs_path - synth_path
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(years, obs_path, color="#B22222", label="California")
ax.plot(years, synth_path, color="#1F4E79", linestyle="--", label="Synthetic California")
ax.axvline(1988, color="grey", linestyle=":")
ax.set_ylabel("Cigarette sales (packs per capita)")
ax.set_title("Dotted line marks Proposition 99 (1988)")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/smoking.csv", clear
destring lnincome beer age15to24, replace force
encode state, gen(state_id)
tsset state_id year
quietly summarize state_id if state == "California"
quietly synth cigsale lnincome(1980(1)1988) retprice(1980(1)1988) ///
age15to24(1980(1)1988) beer(1984(1)1988) cigsale(1975) cigsale(1980) ///
cigsale(1988), trunit(`r(mean)') trperiod(1989) mspeperiod(1970(1)1988) ///
resultsperiod(1970(1)2000) keep("synthetic-control_files/stata/_ca_out", replace)
use "synthetic-control_files/stata/_ca_out", clear
twoway (line _Y_treated _time, lcolor("178 34 34")) ///
(line _Y_synthetic _time, lcolor("31 78 121") lpattern(dash)), ///
xline(1988, lpattern(dot)) ytitle("Cigarette sales (packs per capita)") ///
legend(order(1 "California" 2 "Synthetic California")) xtitle("")
quietly graph export "sc_figs/prop99-paths-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
Panel variable: state_id (strongly balanced)
Time variable: year, 1970 to 2000
Delta: 1 unit
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 8.29e+01 | 6.80e-06 | -1.36e+01 | -3.94e+02 | 0.000 | 0.0000 | 1.00e+02
1 | 5.05e-01 | 4.14e-08 | -1.36e+01 | -7.03e+02 | 0.000 | 0.9939 | 3.06e-05
2 | 2.87e-03 | 2.35e-10 | -1.36e+01 | -2.90e+01 | 0.000 | 0.9943 | 2.72e-06
3 | 1.63e-04 | 1.33e-11 | -1.36e+01 | -1.45e+01 | 1.220 | 0.9433 | 5.50e-06
4 | 1.60e-05 | 1.32e-12 | -1.36e+01 | -1.38e+01 | 2.024 | 0.9013 | 9.47e-07
5 | 8.35e-06 | 6.85e-13 | -1.36e+01 | -1.37e+01 | 2.305 | 0.4792 | 4.50e-06
6 | 2.37e-06 | 1.94e-13 | -1.36e+01 | -1.37e+01 | 2.833 | 0.7163 | 6.69e-07
7 | 7.27e-07 | 5.96e-14 | -1.36e+01 | -1.37e+01 | 3.334 | 0.6933 | 2.33e-07
8 | 1.61e-07 | 1.32e-14 | -1.36e+01 | -1.36e+01 | 3.924 | 0.7790 | 3.75e-08
9 | 3.88e-09 | 3.18e-16 | -1.36e+01 | -1.36e+01 | 5.302 | 0.9758 | 1.11e-10
10 | 1.99e-11 | 2.13e-18 | -1.36e+01 | -1.36e+01 | 7.563 | 0.9949 | 2.09e-13
11 | 9.92e-14 | 1.22e-18 | -1.36e+01 | -1.36e+01 | 9.861 | 0.9950 | 1.09e-15
12 | 8.56e-16 | 1.74e-18 | -1.36e+01 | -1.36e+01 | 12.158 | 0.9950 | 5.49e-18
----------------------------------------------------------------------------------
optimization converged
pre_ca <- years_ca <= 1988
data.frame(year = years_ca, gap = gap_ca) |>
ggplot(aes(year, gap)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = 1988, linetype = "dotted", colour = "grey40") +
geom_line(colour = PAL[["treated"]], linewidth = 0.9) +
labs(x = NULL, y = "Gap in packs per capita")
cat("Pre-treatment RMSPE:", round(sqrt(mean(gap_ca[pre_ca]^2)), 3), "\n")Pre-treatment RMSPE: 1.791 cat("Mean post-1988 gap :", round(mean(gap_ca[!pre_ca]), 2), "packs\n")Mean post-1988 gap : -18.72 packscat("Gap in 2000 :", round(gap_ca[years_ca == 2000], 2), "packs\n")Gap in 2000 : -25.43 packspre = years <= 1988
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(1988, color="grey", linestyle=":")
ax.plot(years, gap_ca, color="#B22222")
ax.set_ylabel("Gap in packs per capita")
plt.tight_layout()
plt.show()
print(f"Pre-treatment RMSPE: {round(np.sqrt(np.mean(gap_ca[pre]**2)), 3)}")Pre-treatment RMSPE: 4.628print(f"Mean post-1988 gap : {round(gap_ca[~pre].mean(), 2)} packs")Mean post-1988 gap : -19.92 packsprint(f"Gap in 2000 : {round(gap_ca[years == 2000][0], 2)} packs")Gap in 2000 : -27.69 packsuse "synthetic-control_files/stata/_ca_out", clear // saved by the paths block
gen gap = _Y_treated - _Y_synthetic
twoway line gap _time, yline(0) xline(1988, lpattern(dot)) ///
lcolor("178 34 34") ytitle("Gap in packs per capita") xtitle("")
gen g2 = gap^2 if _time <= 1988
quietly summarize g2
display "Pre-treatment RMSPE: " %6.3f sqrt(r(mean))
quietly summarize gap if _time > 1988
display "Mean post-1988 gap : " %6.2f r(mean) " packs"
quietly summarize gap if _time == 2000
display "Gap in 2000 : " %6.2f r(mean) " packs"
quietly graph export "sc_figs/prop99-gap-stata.svg", replace width(1600)(7 missing values generated)
(19 missing values generated)
Pre-treatment RMSPE: 2.037
Mean post-1988 gap : -17.96 packs
Gap in 2000 : -24.54 packs
The pre-treatment RMSPE is small on a scale where sales run above 100 packs per capita, and the estimated effect grows over time, reaching roughly 25 fewer packs per capita by 2000. That growth pattern is itself substantively meaningful: a tax-plus-campaign intervention would be expected to accumulate, not to jump once and flatten.
Synth also ships quick base-R versions of these two figures, path.plot() and gaps.plot(), which take the same synth.res and dataprep.res objects.
What to Check, and in What Order
- Pre-treatment path: does the synthetic control track the treated unit before treatment?
- Donor weights: which units are carrying the comparison?
- Predictor balance: does the synthetic control improve on the raw donor average?
- Post-treatment gap: how large and persistent is the divergence?
- Placebos: does the treated gap look unusual relative to fake-treated donors?
- Sensitivity: does the conclusion survive reasonable donor exclusions and leave-one-out checks?
If item one fails, item four is not worth reading.
The Optimizer’s Answer Is Not Always Unique
synth() solves the \(V\)-weighted minimization by numerical optimization, and the nested optimization used to choose \(V\) is not guaranteed to have a unique solution. Klößner, Kaul, Pfeifer, and Schieler show that reproductions of the same application can land on different predictor weights, and therefore different donor weights, across software and settings that all claim to implement the same method (Klößner et al. 2018). The weights above are close to, but not identical to, the published figures for exactly this reason. Report the software, version, and optimization settings, and revisit the fit if a small settings change moves the conclusion.
Synth is not the only implementation. tidysynth reimplements the same method with a pipeline interface that some readers find easier to adapt (Dunford 2025).
Generalizing Synthetic Control
Classical SCM assumes a specific and rather restrictive world: one treated unit, one treatment date, treatment that never switches off. Most panel data are not that tidy.
Where the Classical Design Runs Out
| Complication | Why classical SCM strains |
|---|---|
| Several treated units | You could fit a separate synthetic control per treated unit, but there is no single estimand, no shared uncertainty statement, and donor pools overlap. |
| Staggered adoption | Units treated at different dates means “the pre-period” is unit-specific, and a unit treated late is a valid donor early and an invalid one later. |
| Treatment reversal | If treatment switches off, post-treatment periods are not uniformly treated, and the “after” window stops being well defined. |
| Missing observations | The optimizer needs a balanced matrix of pre-treatment outcomes; holes break it. |
| Covariate-dependent trends | Convex weights on outcome paths may not absorb a trend driven by a covariate. |
The unifying fix is to stop thinking about weighting donors and start thinking about imputing the missing potential outcome from a model fitted on the untreated observations.
The Factor-Model View
Write untreated outcomes as an interactive fixed-effects model:
\[Y_{it}(0) = \mu + \alpha_i + \xi_t + \lambda_i' f_t + \varepsilon_{it},\]
where \(\alpha_i\) is a unit effect, \(\xi_t\) a common time effect, \(f_t\) a vector of \(r\) latent time-varying factors, and \(\lambda_i\) the unit’s loadings on them. Additive two-way fixed effects are the special case \(r = 0\); parallel trends is what that assumption buys you.
The generalized synthetic control estimator proceeds in three steps (Xu 2017):
- Fit the factor model on untreated observations only, estimating \(f_t\) and each unit’s \(\lambda_i\).
- For a treated unit, estimate its loadings \(\lambda_i\) from its pre-treatment periods.
- Impute \(\widehat{Y}_{it}(0)\) for treated periods and take the difference from what was observed.
This generalizes classical SCM rather than simply competing with it, though the two are not the same estimator and can disagree, as we show below. Classical SCM matches the pre-treatment path in order to implicitly match \(\lambda_i\); the factor model estimates \(\lambda_i\) directly. Because step 3 is an imputation, everything that broke above now works: multiple treated units simply give more imputations to average, staggered adoption means each unit has its own set of imputed periods, and reversal is handled because treated status is defined observation by observation rather than by a single date.
The fect package implements this family of counterfactual estimators, including the generalized synthetic control procedure (Xu et al. 2026).
| Setting | Starting point in fect |
|---|---|
| Few treated units; block or staggered adoption; no reversal; time components learned from never-treated donors | method = "gsynth" |
| Counterfactual imputation using not-yet-treated observations, including treatment reversal | method = "ife" |
| Standard additive fixed-effects counterfactual | method = "fe" |
| Matrix-completion formulation | method = "mc" |
The choice changes both the source of identifying variation and the uncertainty procedure.
Recovering a Known Effect
fect ships a simulated panel generated from exactly this factor structure, including the true factors, the true loadings, and the true treatment effect, so we can check the estimator against ground truth.2
The generalized SCM estimator across languages
R uses fect; Stata uses fect (same authors; R’s method="gsynth" is Stata’s method("ife")). The official Python port fect_python is an rpy2 wrapper around R’s fect, but its import conflicts with this environment’s rpy2 build, so the Python tabs below use a compact native NumPy implementation of the interactive fixed-effects estimator instead. It recovers the same known effect (about \(5.5\) against a truth of \(5.46\)); it is not a full reimplementation of fect (no cross-validation search or bootstrap inference), which is why the R and Stata tabs remain the reference.
# Local cache of fect::simgsynth (see synthetic-control_files/data/PROVENANCE.md).
simgsynth <- read.csv("synthetic-control_files/data/simgsynth.csv")
cat("Units:", length(unique(simgsynth$id)),
"| periods:", length(unique(simgsynth$time)),
"| treated observations:", sum(simgsynth$D == 1), "\n")Units: 50 | periods: 30 | treated observations: 50 cat("True ATT:", round(mean(simgsynth$eff[simgsynth$D == 1]), 3), "\n")True ATT: 5.462 simgsynth = pd.read_csv("synthetic-control_files/data/simgsynth.csv")
print("Units:", simgsynth["id"].nunique(),
"| periods:", simgsynth["time"].nunique(),
"| treated observations:", int((simgsynth["D"] == 1).sum()))Units: 50 | periods: 30 | treated observations: 50print("True ATT:", round(simgsynth.loc[simgsynth["D"] == 1, "eff"].mean(), 3))True ATT: 5.462import delimited "synthetic-control_files/data/simgsynth.csv", clear
quietly levelsof id, local(ids)
quietly levelsof time, local(tms)
quietly count if d == 1
display "Units: " `:word count `ids'' " | periods: " `:word count `tms'' ///
" | treated observations: " r(N)
quietly summarize eff if d == 1
display "True ATT: " %6.3f r(mean)(encoding automatically selected: ISO-8859-1)
(15 vars, 1,500 obs)
Units: 50 | periods: 30 | treated observations: 50
True ATT: 5.462We fit the estimator on the simulated panel. The factor model is estimated on the never treated units, each treated unit’s untreated outcome is imputed from its pre-treatment factor loadings, and the difference between what we observe and that imputation is the estimated effect. The fit itself prints nothing; we read its numbers in the next step. The R and Stata tabs call fect; the Python tab uses the compact native estimator defined above.
fit_gs <- fect(
Y ~ D + X1 + X2,
data = simgsynth,
index = c("id", "time"),
method = "gsynth",
force = "two-way",
CV = TRUE,
r = c(0, 5),
se = TRUE,
nboots = 200,
parallel = FALSE
)# Compact native interactive-fixed-effects (gsynth) estimator. Fit the factor
# model on the never-treated units, impute each treated unit's untreated outcome
# from its pre-treatment loadings, and bootstrap control units for confidence
# bands. Handles a common date or staggered adoption, with or without covariates.
def _core(Y, D, Xs, ctrl, r, iters, tol):
Yc = Y[ctrl]; Xc = [X[ctrl] for X in Xs]; T = Y.shape[1]
beta = np.zeros(len(Xs)); mu = 0.0; xi = np.zeros(T); F = np.zeros((T, r)); prev = np.inf
for _ in range(iters): # alternate: covariates / FE / factors
W = Yc - sum(b * X for b, X in zip(beta, Xc))
mu = W.mean(); alpha = (W - mu).mean(1); xi = (W - mu - alpha[:, None]).mean(0)
R = W - mu - alpha[:, None] - xi[None, :]
_, S, Vt = np.linalg.svd(R, full_matrices=False)
F = Vt[:r].T * S[:r]
fac = R @ F @ np.linalg.inv(F.T @ F) @ F.T
Ytil = Yc - mu - alpha[:, None] - xi[None, :] - fac
if Xs:
A = np.column_stack([X.ravel() for X in Xc])
beta, *_ = np.linalg.lstsq(A, Ytil.ravel(), rcond=None)
obj = ((Ytil - fac) ** 2).sum()
if abs(prev - obj) < tol:
break
prev = obj
return beta, mu, xi, F
def gsynth_fit(df, y="Y", d="D", unit="id", time="time", x=("X1", "X2"),
r=2, iters=200, tol=1e-7, nboot=200, seed=1):
units = np.sort(df[unit].unique()); times = np.sort(df[time].unique())
N, T = len(units), len(times)
ui = {u: i for i, u in enumerate(units)}; ti = {t: j for j, t in enumerate(times)}
Y = np.zeros((N, T)); D = np.zeros((N, T), int); Xs = [np.zeros((N, T)) for _ in x]
for _, row in df.iterrows():
i, j = ui[row[unit]], ti[row[time]]; Y[i, j] = row[y]; D[i, j] = int(row[d])
for k, xv in enumerate(x):
Xs[k][i, j] = row[xv]
ctrl = np.where(D.sum(1) == 0)[0]; trt = np.where(D.sum(1) > 0)[0]
def effects(ctrl_): # returns ATT, per-period eff, factors, event study
beta, mu, xi, F = _core(Y, D, Xs, ctrl_, r, iters, tol)
per = [[] for _ in range(T)]; evt = {}; num = den = 0.0
for i in trt:
pre = np.where(D[i] == 0)[0]; post = np.where(D[i] == 1)[0]; a = int(post.min())
yi = Y[i] - sum(b * X[i] for b, X in zip(beta, Xs)) - mu - xi
Z = np.column_stack([np.ones(len(pre)), F[pre]])
coef, *_ = np.linalg.lstsq(Z, yi[pre], rcond=None)
y0 = sum(b * X[i] for b, X in zip(beta, Xs)) + mu + xi + coef[0] + F @ coef[1:]
e = Y[i] - y0
for j in range(T):
per[j].append(e[j]); evt.setdefault(j - a, []).append(e[j])
num += e[post].sum(); den += len(post)
eff = np.array([np.mean(v) if v else np.nan for v in per])
ek = np.array(sorted(evt)); em = np.array([np.mean(evt[k]) for k in ek])
return num / den, eff, F, ek, em
att, eff, F, ek, em = effects(ctrl)
rng = np.random.default_rng(seed); bcal, bevt = [], []
for _ in range(nboot): # resample control units for bands
try:
_, e2, _, k2, m2 = effects(rng.choice(ctrl, len(ctrl), replace=True))
bcal.append(e2); bevt.append(np.interp(ek, k2, m2))
except Exception:
pass
lo, hi = np.nanpercentile(bcal, [2.5, 97.5], axis=0)
elo, ehi = np.nanpercentile(bevt, [2.5, 97.5], axis=0)
return dict(att=att, eff=eff, eff_lo=lo, eff_hi=hi, F=F, times=times,
evt_k=ek, evt_m=em, evt_lo=elo, evt_hi=ehi,
T0=int(df.loc[df[d] == 1, time].min()) - 1, r=r)
simgsynth = pd.read_csv("synthetic-control_files/data/simgsynth.csv")
fit_gs = gsynth_fit(simgsynth, r=2) # r selected by CV in the R/Stata tabsimport delimited "synthetic-control_files/data/simgsynth.csv", clear
* Stata fect (Liu-Wang-Xu). R's method="gsynth" is interactive fixed effects,
* which Stata fect calls "ife"; r() is a single integer (the true value is 2).
fect y, treat(d) unit(id) time(time) cov(x1 x2) method("ife") r(2)(encoding automatically selected: ISO-8859-1)
(15 vars, 1,500 obs)
Balanced Panel DataNow we read the fit. The estimator reports how many latent factors it used and the average effect across all treated observations. Because this panel was simulated with a known effect, we can print the truth alongside and check that the estimate lands close to it.
cat("Factors selected by cross-validation:", fit_gs$r.cv, "(true value: 2)\n")Factors selected by cross-validation: 2 (true value: 2)cat("Estimated ATT:", round(fit_gs$est.avg[1], 3),
" SE:", round(fit_gs$est.avg[2], 3), "\n")Estimated ATT: 5.543 SE: 0.228 cat("True ATT :", round(mean(simgsynth$eff[simgsynth$D == 1]), 3), "\n")True ATT : 5.462 print("Factors used:", fit_gs["r"], "(true value: 2)")Factors used: 2 (true value: 2)print("Estimated ATT:", round(fit_gs["att"], 3))Estimated ATT: 5.543print("True ATT :", round(simgsynth.loc[simgsynth.D == 1, "eff"].mean(), 3))True ATT : 5.462import delimited "synthetic-control_files/data/simgsynth.csv", clear
quietly fect y, treat(d) unit(id) time(time) cov(x1 x2) method("ife") r(2)
matrix A = e(ATT) // columns: ATT, N (treated observations)
display "Estimated ATT: " %6.3f A[1,1](encoding automatically selected: ISO-8859-1)
(15 vars, 1,500 obs)
Estimated ATT: 5.579Cross-validation recovers the correct number of latent factors, and the estimate lands within a standard error of the truth. We use 200 bootstrap replications here to keep the guide quick to render; report more in real work.
The average hides the shape of the response. The next plot shows the effect period by period with its bootstrap interval. Read it like a pre-treatment fit check: the estimates should sit near zero before the intervention, with intervals covering zero, and then lift away from it once treatment begins.
plot(fit_gs, main = "")
periods = fit_gs["times"]
eff = fit_gs["eff"]
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(fit_gs["T0"] + 0.5, color="grey", linestyle=":")
ax.errorbar(periods, eff,
yerr=[eff - fit_gs["eff_lo"], fit_gs["eff_hi"] - eff],
fmt="o", color="#B22222", ecolor="0.6", capsize=2, markersize=3)
ax.set_xlabel("Period")
ax.set_ylabel("Estimated effect")
plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/simgsynth.csv", clear
quietly fect y, treat(d) unit(id) time(time) cov(x1 x2) method("ife") r(2)
* e(ATTs) holds the event-study: column 1 is event time, column 2 the effect.
matrix E = e(ATTs)
clear
svmat E
twoway line E2 E1, yline(0) xline(0, lpattern(dot)) ///
xtitle("Periods since treatment") ytitle("Estimated effect")
quietly graph export "sc_figs/fig-gsynth-gap-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-1)
(15 vars, 1,500 obs)
number of observations will be reset to 30
Press any key to continue, or Break to abort
Number of observations (_N) was 0, now 30.
The estimator also exposes the machinery behind the imputation, which is worth inspecting rather than trusting.
plot(fit_gs, type = "factors", main = "")
fig, ax = plt.subplots(figsize=(8, 4.5))
for k in range(fit_gs["F"].shape[1]):
ax.plot(fit_gs["times"], fit_gs["F"][:, k], marker="o", markersize=3,
label=f"Factor {k + 1}")
ax.axhline(0, color="grey")
ax.set_xlabel("Period")
ax.set_ylabel("Latent factor")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()
* Stata fect does not expose a built-in latent-factor plot; the factors are an
* internal object. This view is R-specific (gsynth/fect in R).
display "Latent-factor plot is R-specific; see the R tab."Latent-factor plot is R-specific; see the R tab.A Staggered-Adoption Example
fect also ships the election-day registration panel used in the original paper: state-level voter turnout, with states adopting reform in different years and never reversing.3
The first step loads the panel and, for every state that ever adopts election day registration, records the first year it did. That gives a quick census of how many states are treated and when the treatment starts spreading.
# Local cache of fect::turnout (see synthetic-control_files/data/PROVENANCE.md).
turnout <- read.csv("synthetic-control_files/data/turnout.csv")
adopt <- turnout |>
filter(policy_edr == 1) |>
group_by(abb) |>
summarise(first_year = min(year), .groups = "drop") |>
arrange(first_year)
cat("States:", length(unique(turnout$abb)),
"| ever treated:", nrow(adopt), "\n")States: 47 | ever treated: 9 print(as.data.frame(adopt)) abb first_year
1 ME 1976
2 MN 1976
3 WI 1976
4 ID 1996
5 NH 1996
6 WY 1996
7 IA 2008
8 MT 2008
9 CT 2012turnout = pd.read_csv("synthetic-control_files/data/turnout.csv")
adopt = (turnout[turnout.policy_edr == 1]
.groupby("abb", as_index=False)["year"].min()
.rename(columns={"year": "first_year"})
.sort_values("first_year"))
print("States:", turnout.abb.nunique(), "| ever treated:", len(adopt))States: 47 | ever treated: 9print(adopt.to_string(index=False))abb first_year
ME 1976
MN 1976
WI 1976
ID 1996
NH 1996
WY 1996
IA 2008
MT 2008
CT 2012import delimited "synthetic-control_files/data/turnout.csv", clear
quietly levelsof abb, local(st)
preserve
keep if policy_edr == 1
collapse (min) first_year = year, by(abb)
sort first_year
display "States: " `:word count `st'' " | ever treated: " _N
list abb first_year, noobs
restore(encoding automatically selected: ISO-8859-1)
(6 vars, 1,128 obs)
(1,078 observations deleted)
States: 47 | ever treated: 9
+----------------+
| abb first_~r |
|----------------|
| WI 1976 |
| ME 1976 |
| MN 1976 |
| ID 1996 |
| WY 1996 |
|----------------|
| NH 1996 |
| IA 2008 |
| MT 2008 |
| CT 2012 |
+----------------+Nine states adopt across four decades, which is precisely the staggered structure classical SCM cannot handle in one estimand.
The grid below makes the staggering visible. Each cell marks whether a state has election day registration in a given year, and the states are ordered by how long they spend treated, so the block of treated cells grows toward the top.
turnout |>
mutate(status = factor(policy_edr, levels = c(0, 1),
labels = c("Untreated", "Treated"))) |>
ggplot(aes(year, reorder(abb, policy_edr, sum), fill = status)) +
geom_tile(colour = "white", linewidth = 0.1) +
scale_fill_manual(values = c(Untreated = "grey85",
Treated = PAL[["treated"]])) +
labs(x = NULL, y = NULL, fill = NULL) +
theme(axis.text.y = element_text(size = 6))
grid = turnout.pivot(index="abb", columns="year", values="policy_edr")
grid = grid.loc[grid.sum(axis=1).sort_values().index] # order by exposure
fig, ax = plt.subplots(figsize=(8, 5))
ax.pcolormesh(grid.columns, range(len(grid.index)), grid.values,
cmap=plt.cm.colors.ListedColormap(["grey", "#B22222"]))
ax.set_yticks(range(len(grid.index)))
ax.set_yticklabels(grid.index, fontsize=6)
ax.set_xlabel("")
plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/turnout.csv", clear
encode abb, gen(state)
* first-draft heatmap via scatter of treated cells (a heatplot package gives a
* cleaner grid: ssc install heatplot).
twoway (scatter state year if policy_edr == 0, msymbol(square) mcolor(gs13) msize(small)) ///
(scatter state year if policy_edr == 1, msymbol(square) mcolor("178 34 34") msize(small)), ///
legend(order(1 "Untreated" 2 "Treated")) ytitle("") xtitle("") ylabel(none)
quietly graph export "sc_figs/fig-turnout-status-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-1)
(6 vars, 1,128 obs)
With the design mapped out, we fit the generalized synthetic control estimator. Instead of weighting donor states, it estimates a factor model on the states that are not yet treated and uses it to impute what each treated state’s turnout would have been without the reform. A single model then covers every adoption date at once. The fit is silent here; we read its numbers in the next step.
fit_turnout <- fect(
turnout ~ policy_edr,
data = turnout,
index = c("abb", "year"),
method = "gsynth",
force = "two-way",
CV = TRUE,
r = c(0, 5),
se = TRUE,
nboots = 200,
parallel = FALSE
)# Same native estimator, now with staggered adoption and no covariates.
turnout = pd.read_csv("synthetic-control_files/data/turnout.csv")
fit_turnout = gsynth_fit(turnout, y="turnout", d="policy_edr",
unit="abb", time="year", x=(), r=2, nboot=200)import delimited "synthetic-control_files/data/turnout.csv", clear
encode abb, gen(state)
* R's method="gsynth" is interactive fixed effects; in Stata fect that is "ife".
* Stata fect takes a single integer r(); r=2 is the count CV selects here.
fect turnout, treat(policy_edr) unit(state) time(year) method("ife") r(2)(encoding automatically selected: ISO-8859-1)
(6 vars, 1,128 obs)
Balanced Panel DataThe fit reports a single average effect across all treated state years, together with the number of latent factors the procedure kept. Because each Stata block runs in its own session, the Stata tab refits the model before reading the result.
cat("Factors selected:", fit_turnout$r.cv, "\n")Factors selected: 2 cat("Estimated ATT:", round(fit_turnout$est.avg[1], 3),
" SE:", round(fit_turnout$est.avg[2], 3), "\n")Estimated ATT: 5.13 SE: 1.856 print("Factors used:", fit_turnout["r"])Factors used: 2print("Estimated ATT:", round(fit_turnout["att"], 3))Estimated ATT: 5.13import delimited "synthetic-control_files/data/turnout.csv", clear
encode abb, gen(state)
quietly fect turnout, treat(policy_edr) unit(state) time(year) method("ife") r(2)
matrix A = e(ATT) // columns: ATT, N (treated observations)
display "Estimated ATT: " %6.3f A[1,1](encoding automatically selected: ISO-8859-1)
(6 vars, 1,128 obs)
Estimated ATT: 5.203Averages hide timing. Aligning every state to its own adoption year and plotting the effect period by period shows whether turnout responds at once or gradually, and, just as important, whether anything is already moving before the reform. Flat pre-adoption estimates are the staggered analogue of a clean pre-treatment fit.
plot(fit_turnout, main = "")
k, m = fit_turnout["evt_k"], fit_turnout["evt_m"]
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(0, color="grey", linestyle=":")
ax.errorbar(k, m, yerr=[m - fit_turnout["evt_lo"], fit_turnout["evt_hi"] - m],
fmt="o", color="#B22222", ecolor="0.6", capsize=2, markersize=3)
ax.set_xlabel("Periods since adoption")
ax.set_ylabel("Effect on turnout")
plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/turnout.csv", clear
encode abb, gen(state)
quietly fect turnout, treat(policy_edr) unit(state) time(year) method("ife") r(2)
* e(ATTs) holds the event-study: column 1 is event time, column 2 the effect.
matrix E = e(ATTs)
clear
svmat E
twoway line E2 E1, yline(0) xline(0, lpattern(dot)) ///
xtitle("Periods since adoption") ytitle("Effect on turnout")
quietly graph export "sc_figs/fig-turnout-gap-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-1)
(6 vars, 1,128 obs)
number of observations will be reset to 33
Press any key to continue, or Break to abort
Number of observations (_N) was 0, now 33.
Because effects are aligned to each state’s own adoption year, period zero means something different in calendar time for each unit. The pre-treatment periods now serve the role the pre-treatment path served in classical SCM: if estimated effects are already nonzero before adoption, the imputation model is not capturing the untreated trajectory.
Does the General Tool Replace the Specific One?
Gsynth handles everything classical SCM handles, plus cases classical SCM cannot express. A natural conclusion is that you should simply always use it. That conclusion is wrong, and Proposition 99 shows why.
The two are not the same estimator. Classical SCM puts convex weights on donor units, so the counterfactual is a blend of real states and cannot leave their range. Gsynth fits a factor model to the untreated observations, estimates California’s loadings from its pre-period, and imputes the counterfactual. Nothing constrains that imputation to resemble any real state. Run both on the same data and they disagree.
The disagreement is not a small numerical wobble. It depends almost entirely on one tuning choice: how many latent factors the model is allowed.
smoking_fect <- smoking_df |>
mutate(treat = as.integer(state == "California" & year >= 1989))
r_sweep <- lapply(0:5, function(rr) {
set.seed(20260720)
invisible(capture.output(
f <- fect(cigsale ~ treat, data = smoking_fect,
index = c("state_id", "year"), method = "gsynth",
force = "two-way", CV = FALSE, r = rr,
se = TRUE, nboots = 200, parallel = FALSE)
))
pos <- which(f$id == ca_id)
g <- f$Y.dat[, pos] - f$Y.ct[, pos]
data.frame(Factors = rr,
ATT = round(f$est.avg[1], 2),
SE = round(f$est.avg[2], 2),
`Pre-RMSE` = round(sqrt(mean(g[f$rawtime <= 1988]^2)), 2),
check.names = FALSE)
}) |> bind_rows()
r_sweep$`Classical SCM ATT` <- round(mean(gap_ca[!pre_ca]), 2)
r_sweep Factors ATT SE Pre-RMSE Classical SCM ATT
1 0 -27.35 2.66 7.16 -18.72
2 1 -13.89 3.39 4.06 -18.72
3 2 -0.40 9.94 1.36 -18.72
4 3 -1.87 8.22 1.35 -18.72
5 4 -1.90 8.13 1.35 -18.72
6 5 -0.45 8.03 1.35 -18.72smoking_fect = smoking.assign(
treat=((smoking.state == "California") & (smoking.year >= 1989)).astype(int)
)
rows = []
for rr in range(0, 6):
f = Fect(data=smoking_fect, Y="cigsale", D="treat",
index=["state_id", "year"], method="gsynth", force="two-way",
CV=False, r=(rr, rr), se=True, nboots=200, parallel=False)
rows.append({"Factors": rr, "ATT": round(f.est_avg[0], 2),
"SE": round(f.est_avg[1], 2)})
print(pd.DataFrame(rows).to_string(index=False))import delimited "synthetic-control_files/data/smoking.csv", clear
destring lnincome beer age15to24, replace force
encode state, gen(state_id)
gen treat = (state == "California" & year >= 1989)
tempname R
postfile `R' int factors double att using "synthetic-control_files/stata/_rsweep", replace
forvalues rr = 0/5 {
if `rr' == 0 {
quietly fect cigsale, treat(treat) unit(state_id) time(year) method("fe")
}
else {
quietly fect cigsale, treat(treat) unit(state_id) time(year) method("ife") r(`rr')
}
matrix A = e(ATT)
post `R' (`rr') (A[1,1])
}
postclose `R'
use "synthetic-control_files/stata/_rsweep", clear
list, noobs(encoding automatically selected: ISO-8859-2)
(7 vars, 1,209 obs)
lnincome: contains nonnumeric characters; replaced as double
(195 missing values generated)
beer: contains nonnumeric characters; replaced as double
(663 missing values generated)
age15to24: contains nonnumeric characters; replaced as double
(390 missing values generated)
+----------------------+
| factors att |
|----------------------|
| 0 -27.349111 |
| 1 -13.870756 |
| 2 -.72578989 |
| 3 -23.052441 |
| 4 -22.224885 |
|----------------------|
| 5 -25.717255 |
+----------------------+Read the first two columns together. As the model is given more factors, its pre-treatment fit improves steadily, from an RMSE above 7 down to about 1.35, which is better than classical SCM manages. And as the fit improves, the estimated effect collapses from \(-27\) toward zero, with standard errors wide enough to cover almost anything.
That is overfitting, and it is worth being precise about the mechanism. California is one unit with nineteen pre-treatment periods. Four latent factors give the model enough freedom to reproduce almost any pre-1988 path, and that same freedom carries into the post-1988 period, where it absorbs the treatment effect instead of revealing it. Cross-validation does not save you here: before we fixed a seed, repeated runs of the same specification selected different factor counts and returned estimates ranging from \(-1.9\) to \(-13.9\). Constraining the loadings with loading.bound = "simplex" changes the answer again rather than settling it.
set.seed(20260720)
invisible(capture.output(
fit_gs_ca <- fect(cigsale ~ treat, data = smoking_fect,
index = c("state_id", "year"), method = "gsynth",
force = "two-way", CV = TRUE, r = c(0, 5),
se = TRUE, nboots = 200, parallel = FALSE)
))
pos_ca <- which(fit_gs_ca$id == ca_id)
data.frame(
year = years_ca,
`California (observed)` = obs_ca,
`Synthetic control (convex)` = synth_ca,
`Gsynth imputation` = fit_gs_ca$Y.ct[, pos_ca],
check.names = FALSE
) |>
pivot_longer(-year, names_to = "series", values_to = "cigsale") |>
ggplot(aes(year, cigsale, colour = series, linetype = series)) +
geom_vline(xintercept = 1988, linetype = "dotted", colour = "grey40") +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = c(`California (observed)` = "grey25",
`Synthetic control (convex)` = PAL[["synth"]],
`Gsynth imputation` = PAL[["treated"]])) +
labs(x = NULL, y = "Cigarette sales (packs per capita)",
colour = NULL, linetype = NULL,
subtitle = "Dotted line marks Proposition 99 (1988)")
# convex synthetic path (from pysyncon, above) vs the gsynth imputation
fit_gs_ca = Fect(data=smoking_fect, Y="cigsale", D="treat",
index=["state_id", "year"], method="gsynth", force="two-way",
CV=True, r=(0, 5), se=True, nboots=200, parallel=False)
gsynth_ct = fit_gs_ca.Y_ct_for("California") # imputed counterfactual
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot(years, obs_path, color="0.25", label="California (observed)")
ax.plot(years, synth_path, color="#1F4E79", label="Synthetic control (convex)")
ax.plot(years, gsynth_ct, color="#B22222", label="Gsynth imputation")
ax.axvline(1988, color="grey", linestyle=":")
ax.legend(frameon=False)
plt.show()* Overlaying the convex synth path (from the Prop 99 fit, _ca_out.dta) and the
* gsynth imputation requires exporting both counterfactuals; this comparison is
* cleanest in R. First-draft placeholder.
display "See the R tab: overlay of convex SCM vs gsynth imputation for California."See the R tab: overlay of convex SCM vs gsynth imputation for California.The lesson generalizes beyond this dataset, and it qualifies advice given earlier in this guide. Pre-treatment fit is necessary but not sufficient. An estimator flexible enough to fit any pre-treatment path tells you nothing by fitting yours. What made the pre-treatment fit meaningful in the classical case was that convexity made a good fit hard to achieve: only a genuinely similar blend of real states could produce one. Remove the constraint and a good fit stops being evidence. The proportion-outcome simulation made the same point from a different direction, where the unconstrained estimator also fit better and also produced counterfactuals that could not be true.
So the choice is not “general tool versus specific tool.” It is a question about your design:
| Your design | Use |
|---|---|
| One treated unit, clean treatment date, adequate pre-period | Classical SCM. Convexity is doing real work, and the weights are reportable. |
| Multiple treated units, staggered adoption, reversal, or missing data | Gsynth or another fect estimator. Classical SCM cannot express the design at all. |
| One treated unit, but poor donor support | Augmented SCM, with the extrapolation stated openly. |
The turnout example earlier is where gsynth earns its keep: nine states adopting across four decades is a design classical SCM cannot represent in a single estimand. Proposition 99 is the opposite case, and there the extra flexibility is a liability rather than an upgrade.
Choices You Must Report
Report the treatment coding, treatment timing, panel indices, method, fixed effects, candidate factor range, cross-validation settings, inference procedure, bootstrap count, missing-data handling, and package versions. Also state whether treatment can switch on and off.
Two fect specifics are easy to get wrong. loading.bound = "simplex", which keeps treated-unit loadings inside the donor support and is closest to the classical SCM intuition, and vartype = "parametric" both require time.component.from = "nevertreated"; without it the call errors or silently falls back. And fect objects have no summary() method, so use print().
Augmented Synthetic Controls
Classical SCM works best when a convex combination of donors reproduces the treated unit. When the treated unit sits outside the donor convex hull, as in the poor-support scenario above, no set of nonnegative weights summing to one can close the gap.
Augmented SCM adds an outcome-model correction to the classical estimate (Ben-Michael, Feller, and Rothstein 2021):
\[\widehat{Y}_{1t}^{AugSCM}(0) = \widehat{Y}_{1t}^{SCM}(0) + \widehat{\text{bias correction}},\]
where the correction, commonly from a ridge regression, estimates the remaining pre-treatment imbalance and adjusts for it. The logic is bias-corrected matching: build the most transparent synthetic control available, measure what it still gets wrong before treatment, and correct by that amount.
| What augmentation buys | What it costs |
|---|---|
| Better fit when exact classical fit is infeasible | More model dependence |
| Explicit bias correction | Possible negative donor weights |
| A controlled form of extrapolation | Weaker “weighted average of real donors” intuition |
How Much Does It Help?
We can measure this directly. Using the same factor-model simulation, we vary how far outside the donor hull the treated unit sits and compare classical SCM against ridge-augmented SCM, with the true effect known.
Simulation design (click to expand)
# One donor pool, treated loadings 1.1 / 1.5 / 2.6 (just outside to far outside
# the hull). Written to disk so all three languages read identical panels.
build_panel <- function(lam_val, seed = 20260720) {
set.seed(seed)
lam <- cbind(runif(J, 0, 1), runif(J, 0, 1))
alpha <- rnorm(J, 0, 3)
Y <- matrix(NA_real_, nrow = J + 1, ncol = TT)
for (i in 1:J) {
Y[i + 1, ] <- 120 + alpha[i] +
8 * lam[i, 1] * f1 + 8 * lam[i, 2] * f2 + rnorm(TT, 0, 1)
}
Y[1, ] <- 120 + 8 * lam_val * f1 + 8 * lam_val * f2 + rnorm(TT, 0, 1) + tau
data.frame(id = rep(1:(J + 1), each = TT),
time = rep(1:TT, times = J + 1),
Y = as.vector(t(Y))) |>
mutate(trt = as.integer(id == 1 & time > T0))
}
for (v in c(1.1, 1.5, 2.6)) {
tag <- sub("[.]", "p", sprintf("%.1f", v))
write.csv(build_panel(v),
file.path(dir_sim, paste0("sim_aug_", tag, ".csv")), row.names = FALSE)
}post_att <- function(a) {
s <- summary(a)$att
mean(s$Estimate[s$Time > T0], na.rm = TRUE)
}
sweep_res <- lapply(c(1.1, 1.5, 2.6), function(v) {
tag <- sub("[.]", "p", sprintf("%.1f", v))
pan <- read.csv(paste0("synthetic-control_files/data/sim_aug_", tag, ".csv"))
a_scm <- augsynth::augsynth(Y ~ trt, unit = id, time = time, data = pan,
progfunc = "none", scm = TRUE)
a_ridge <- augsynth::augsynth(Y ~ trt, unit = id, time = time, data = pan,
progfunc = "ridge", scm = TRUE)
data.frame(
`Treated loading` = v,
`SCM ATT` = round(post_att(a_scm), 2),
`SCM imbalance` = round(a_scm$l2_imbalance, 2),
`Augmented ATT` = round(post_att(a_ridge), 2),
`Aug. imbalance` = round(a_ridge$l2_imbalance, 2),
`Min weight` = round(min(a_ridge$weights), 3),
check.names = FALSE
)
}) |> bind_rows()
sweep_res$`True ATT` <- round(mean(tau[(T0 + 1):TT]), 2)
sweep_res Treated loading SCM ATT SCM imbalance Augmented ATT Aug. imbalance Min weight
1 1.1 -9.24 8.20 -9.92 3.03 -0.196
2 1.5 -7.30 17.98 -9.27 6.56 -0.251
3 2.6 -1.48 49.06 -9.03 8.80 -0.700
True ATT
1 -9.6
2 -9.6
3 -9.6from pysyncon import AugSynth
rows = []
for v in [1.1, 1.5, 2.6]:
tag = f"{v:.1f}".replace(".", "p")
pan = pd.read_csv(f"synthetic-control_files/data/sim_aug_{tag}.csv")
donors = [i for i in pan.id.unique() if i != 1]
dp = Dataprep(
foo=pan, dependent="Y", unit_variable="id", time_variable="time",
treatment_identifier=1, controls_identifier=donors,
predictors=[], predictors_op="mean",
special_predictors=[("Y", [t], "mean") for t in range(1, T0 + 1)],
time_predictors_prior=range(1, T0 + 1), time_optimize_ssr=range(1, T0 + 1),
)
aug = AugSynth()
aug.fit(dp)
rows.append({"Treated loading": v, "Augmented ATT": round(aug.att(range(T0 + 1, TT + 1))["att"], 2)})
print(pd.DataFrame(rows).to_string(index=False)) Treated loading Augmented ATT
1.1 -9.79
1.5 -9.64
2.6 -9.63* Stata augmented SCM via allsynth (Wiltshire; needs distinct, elasticregress).
* allsynth's kept file has `gap` (classical) and `gap_bc` (ridge-augmented).
tempname A
postfile `A' double loading double scm_att double aug_att using "synthetic-control_files/stata/_augsweep", replace
foreach lv in 1.1 1.5 2.6 {
local tag = subinstr("`lv'", ".", "p", .)
import delimited "synthetic-control_files/data/sim_aug_`tag'.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly allsynth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20) ///
bcorrect(merge ridge) keep("synthetic-control_files/stata/_a", replace)
preserve
use "synthetic-control_files/stata/_a", clear
quietly summarize gap if _time > 20
scalar sa = r(mean)
quietly summarize gap_bc if _time > 20
scalar aa = r(mean)
restore
post `A' (`lv') (sa) (aa)
}
postclose `A'
use "synthetic-control_files/stata/_augsweep", clear
format scm_att aug_att %6.2f
list, noobs 10. preserve
11. use "synthetic-control_files/stata/_a", clear
12. quietly summarize gap if _time > 20
13. scalar sa = r(mean)
14. quietly summarize gap_bc if _time > 20
15. scalar aa = r(mean)
16. restore
17. post `A' (`lv') (sa) (aa)
18. }
(encoding automatically selected: ISO-8859-1)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(encoding automatically selected: ISO-8859-1)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
(encoding automatically selected: ISO-8859-1)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
+-----------------------------+
| loading scm_att aug_att |
|-----------------------------|
| 1.1 -9.24 -9.29 |
| 1.5 -7.30 -7.75 |
| 2.6 -1.48 -3.55 |
+-----------------------------+
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 2.57e-08 | -3.84e+02 | -6.34e+02 | 0.188 | 0.0000 | 1.00e+02
1 | 6.51e-01 | 2.48e-10 | -3.84e+02 | -8.38e+02 | 0.000 | 0.9903 | 7.75e-05
2 | 3.98e-03 | 1.52e-12 | -3.84e+02 | -3.99e+02 | 1.390 | 0.9939 | 3.34e-06
3 | 1.55e-04 | 5.90e-14 | -3.84e+02 | -3.84e+02 | 2.798 | 0.9611 | 3.99e-06
4 | 2.55e-05 | 9.74e-15 | -3.84e+02 | -3.84e+02 | 3.425 | 0.8347 | 2.85e-06
5 | 4.39e-07 | 1.67e-16 | -3.84e+02 | -3.84e+02 | 4.501 | 0.9828 | 7.07e-09
6 | 2.12e-08 | 8.13e-18 | -3.84e+02 | -3.84e+02 | 5.604 | 0.9516 | 4.75e-09
7 | 6.74e-10 | 2.75e-19 | -3.84e+02 | -3.84e+02 | 6.978 | 0.9683 | 1.60e-10
8 | 8.46e-12 | 7.42e-20 | -3.84e+02 | -3.84e+02 | 8.289 | 0.9874 | 1.06e-12
9 | 6.41e-14 | 7.36e-20 | -3.84e+02 | -3.84e+02 | 9.334 | 0.9924 | 1.88e-14
10 | 5.53e-16 | 7.00e-20 | -3.84e+02 | -3.84e+02 | 10.778 | 0.9950 | 7.38e-16
11 | 4.37e-16 | 7.53e-20 | -3.84e+02 | -3.84e+02 | 13.072 | 0.9950 | 2.66e-17
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 4.49e-08 | -3.64e+02 | -6.14e+02 | 0.165 | 0.0000 | 1.00e+02
1 | 9.72e-01 | 6.48e-10 | -3.64e+02 | -8.16e+02 | 0.000 | 0.9856 | 1.73e-04
2 | 9.26e-03 | 6.18e-12 | -3.63e+02 | -3.86e+02 | 1.208 | 0.9905 | 8.77e-06
3 | 4.50e-04 | 3.00e-13 | -3.63e+02 | -3.65e+02 | 2.516 | 0.9514 | 9.07e-06
4 | 6.51e-05 | 4.34e-14 | -3.64e+02 | -3.64e+02 | 3.162 | 0.8555 | 3.94e-06
5 | 7.41e-06 | 4.94e-15 | -3.64e+02 | -3.64e+02 | 3.973 | 0.8861 | 5.50e-07
6 | 1.23e-07 | 8.21e-17 | -3.64e+02 | -3.64e+02 | 5.557 | 0.9834 | 1.77e-09
7 | 1.47e-09 | 9.64e-19 | -3.64e+02 | -3.64e+02 | 6.780 | 0.9881 | 2.39e-11
8 | 1.02e-11 | 5.94e-20 | -3.64e+02 | -3.64e+02 | 7.864 | 0.9931 | 4.82e-13
9 | 5.10e-14 | 6.07e-20 | -3.64e+02 | -3.64e+02 | 9.752 | 0.9950 | 2.07e-14
10 | 5.09e-16 | 6.93e-20 | -3.64e+02 | -3.64e+02 | 12.051 | 0.9950 | 2.67e-16
----------------------------------------------------------------------------------
optimization converged
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 1.43e-07 | -2.85e+02 | -5.34e+02 | 0.061 | 0.0000 | 1.00e+02
1 | 1.59e+00 | 3.38e-09 | -2.84e+02 | -7.34e+02 | 0.000 | 0.9763 | 4.65e-04
2 | 2.47e-02 | 5.23e-11 | -2.82e+02 | -3.18e+02 | 0.907 | 0.9845 | 2.64e-05
3 | 1.45e-03 | 3.08e-12 | -2.82e+02 | -2.85e+02 | 2.117 | 0.9411 | 2.12e-05
4 | 4.11e-04 | 8.71e-13 | -2.83e+02 | -2.84e+02 | 2.615 | 0.7172 | 3.02e-05
5 | 8.50e-05 | 1.80e-13 | -2.83e+02 | -2.83e+02 | 3.294 | 0.7933 | 5.07e-06
6 | 7.87e-07 | 1.67e-15 | -2.83e+02 | -2.83e+02 | 4.879 | 0.9907 | 2.05e-09
7 | 1.13e-08 | 2.39e-17 | -2.83e+02 | -2.83e+02 | 6.324 | 0.9857 | 1.28e-10
8 | 9.72e-11 | 2.38e-19 | -2.83e+02 | -2.83e+02 | 7.260 | 0.9914 | 1.66e-12
9 | 6.11e-13 | 1.03e-19 | -2.83e+02 | -2.83e+02 | 8.109 | 0.9937 | 1.02e-13
10 | 3.04e-15 | 7.93e-20 | -2.83e+02 | -2.83e+02 | 8.956 | 0.9949 | 9.36e-15
11 | 4.16e-16 | 9.21e-20 | -2.83e+02 | -2.83e+02 | 9.817 | 0.9950 | 1.30e-15
12 | 5.70e-16 | 8.45e-20 | -2.83e+02 | -2.83e+02 | 10.770 | 0.9950 | 1.79e-16
13 | 4.01e-16 | 9.27e-20 | -2.83e+02 | -2.83e+02 | 12.263 | 0.9950 | 2.00e-17
----------------------------------------------------------------------------------
optimization convergedDonor loadings run from 0 to 1, so a treated loading of 1.1 is just outside the hull and 2.6 is far outside. Three things are visible.
First, classical SCM degrades sharply with distance from the donor support, falling from \(-9.24\) to \(-1.48\) against a true effect of \(-9.6\). This is the same failure the diagnostics section documented, now traced across its full range. Second, the augmented estimate stays close to the truth throughout, and its pre-treatment imbalance is consistently a fraction of the unaugmented one. Third, that accuracy is bought with extrapolation: the most negative donor weight grows from \(-0.20\) to \(-0.70\) as the correction works harder, and the synthetic unit becomes correspondingly less interpretable as a blend of real donors.
One caveat about reading this table too optimistically. These data were generated from a factor model, and a ridge correction on pre-treatment outcomes is well suited to exactly that structure, so this is close to a best case for augmentation. It demonstrates that the bias correction does what it claims when the outcome model is roughly right. It is not a promise that augmentation recovers the truth when the model is wrong, and it is not a licence to treat a badly chosen donor pool as fixable.
pan_mid <- read.csv("synthetic-control_files/data/sim_aug_1p5.csv")
a_scm <- augsynth::augsynth(Y ~ trt, unit = id, time = time, data = pan_mid,
progfunc = "none", scm = TRUE)
a_ridge <- augsynth::augsynth(Y ~ trt, unit = id, time = time, data = pan_mid,
progfunc = "ridge", scm = TRUE)
bind_rows(
summary(a_scm)$att |> mutate(Method = "Classical SCM"),
summary(a_ridge)$att |> mutate(Method = "Ridge-augmented")
) |>
filter(!is.na(Estimate)) |>
ggplot(aes(Time, Estimate, colour = Method)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_vline(xintercept = T0, linetype = "dotted", colour = "grey40") +
geom_line(data = data.frame(Time = 1:TT, Estimate = tau),
aes(Time, Estimate), colour = "grey45",
linetype = "dashed", linewidth = 0.7, inherit.aes = FALSE) +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = c("Classical SCM" = PAL[["synth"]],
"Ridge-augmented" = PAL[["treated"]])) +
labs(x = "Period", y = "Estimated effect", colour = NULL,
subtitle = "Dashed grey line is the true effect")
pan = pd.read_csv("synthetic-control_files/data/sim_aug_1p5.csv")
donors = [i for i in pan.id.unique() if i != 1]
dp = Dataprep(
foo=pan, dependent="Y", unit_variable="id", time_variable="time",
treatment_identifier=1, controls_identifier=donors,
predictors=[], predictors_op="mean",
special_predictors=[("Y", [t], "mean") for t in range(1, T0 + 1)],
time_predictors_prior=range(1, T0 + 1), time_optimize_ssr=range(1, T0 + 1),
)
scm = Synth(); scm.fit(dp)
aug = AugSynth(); aug.fit(dp)
periods = list(range(1, TT + 1))
eff_scm = [scm.att(time_period=[t])["att"] for t in periods]
eff_aug = [aug.att(time_period=[t])["att"] for t in periods]
tau_vec = [(-1.2 * (t - T0) if t > T0 else 0) for t in periods]
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="grey")
ax.axvline(T0, color="grey", linestyle=":")
ax.plot(periods, tau_vec, color="0.45", linestyle="--", label="True effect")
ax.plot(periods, eff_scm, color="#1F4E79", label="Classical SCM")
ax.plot(periods, eff_aug, color="#B22222", label="Ridge-augmented")
ax.set_xlabel("Period")
ax.set_ylabel("Estimated effect")
ax.legend(frameon=False)
plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/sim_aug_1p5.csv", clear
tsset id time
local pre ""
forvalues t = 1/20 {
local pre "`pre' y(`t')"
}
quietly allsynth y `pre', trunit(1) trperiod(21) mspeperiod(1(1)20) ///
bcorrect(merge ridge) keep("synthetic-control_files/stata/_aug15", replace)
use "synthetic-control_files/stata/_aug15", clear
gen tau = cond(_time > 20, -1.2 * (_time - 20), 0) // deterministic truth
twoway (line tau _time, lcolor(gs8) lpattern(dash)) ///
(line gap _time, lcolor("31 78 121")) ///
(line gap_bc _time, lcolor("178 34 34")), ///
yline(0) xline(20, lpattern(dot)) ///
legend(order(1 "True effect" 2 "Classical SCM" 3 "Ridge-augmented")) ///
xtitle("Period") ytitle("Estimated effect")
quietly graph export "sc_figs/fig-augsynth-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-1)
(4 vars, 910 obs)
Panel variable: id (strongly balanced)
Time variable: time, 1 to 35
Delta: 1 unit
counter | pri_inf | dual_inf | pri_obj | dual_obj | sigfig | alpha | nu
----------------------------------------------------------------------------------
0 | 6.73e+01 | 4.49e-08 | -3.64e+02 | -6.14e+02 | 0.165 | 0.0000 | 1.00e+02
1 | 9.72e-01 | 6.48e-10 | -3.64e+02 | -8.16e+02 | 0.000 | 0.9856 | 1.73e-04
2 | 9.26e-03 | 6.18e-12 | -3.63e+02 | -3.86e+02 | 1.208 | 0.9905 | 8.77e-06
3 | 4.50e-04 | 3.00e-13 | -3.63e+02 | -3.65e+02 | 2.516 | 0.9514 | 9.07e-06
4 | 6.51e-05 | 4.34e-14 | -3.64e+02 | -3.64e+02 | 3.162 | 0.8555 | 3.94e-06
5 | 7.41e-06 | 4.94e-15 | -3.64e+02 | -3.64e+02 | 3.973 | 0.8861 | 5.50e-07
6 | 1.23e-07 | 8.21e-17 | -3.64e+02 | -3.64e+02 | 5.557 | 0.9834 | 1.77e-09
7 | 1.47e-09 | 9.64e-19 | -3.64e+02 | -3.64e+02 | 6.780 | 0.9881 | 2.39e-11
8 | 1.02e-11 | 5.94e-20 | -3.64e+02 | -3.64e+02 | 7.864 | 0.9931 | 4.82e-13
9 | 5.10e-14 | 6.07e-20 | -3.64e+02 | -3.64e+02 | 9.752 | 0.9950 | 2.07e-14
10 | 5.09e-16 | 6.93e-20 | -3.64e+02 | -3.64e+02 | 12.051 | 0.9950 | 2.67e-16
----------------------------------------------------------------------------------
optimization converged
Report the unaugmented fit, the augmented fit, the regularization choice, whether weights become negative, and whether conclusions change across plausible tuning choices.
Synthetic Control for Proportion Outcomes
“Proportion synthetic control” is not an estimator name. It is a warning that your outcome is bounded: turnout, vaccination coverage, graduation rates, compliance, or the share of forest area deforested all live in \([0, 1]\).
Why Convexity Matters Here
Recall the constraints \(w_j \geq 0\) and \(\sum_j w_j = 1\). Because the synthetic outcome is a convex combination of donor values, it can never fall outside the range of the donors themselves. If every donor’s deforestation share lies in \([0, 1]\), so does the synthetic counterfactual, automatically.
Estimators that relax convexity lose this guarantee. Augmented SCM permits negative weights; unconstrained regression, matrix completion, and factor models extrapolate freely. That flexibility buys better pre-treatment fit and can produce a counterfactual that is arithmetically impossible.
This is measurable rather than hypothetical. We simulate bounded panels where donors fit the treated unit’s pre-period well, then compare a convex synthetic control against an unconstrained least-squares fit on the same data.
Simulation design (click to expand)
set.seed(31415)
Jp <- 10; T0p <- 18; TTp <- 30
pre_p <- 1:T0p
draw_bounded_panel <- function() {
wave <- sin(seq(0, 2 * pi, length.out = TTp))
common <- 0.9 * wave
div <- c(rep(0, T0p), seq(0, 1, length.out = TTp - T0p)) # diverge after pre-period
P <- matrix(NA_real_, nrow = Jp + 1, ncol = TTp)
for (i in 1:Jp) {
s <- runif(1, -1, 1) * 2.5
eta <- qlogis(runif(1, 0.15, 0.35)) + common +
0.05 * rnorm(1) * wave + s * div + rnorm(TTp, 0, 0.01)
P[i + 1, ] <- plogis(eta)
}
P[1, ] <- plogis(qlogis(runif(1, 0.14, 0.30)) + common + rnorm(TTp, 0, 0.01))
P
}
# 200 draws, stacked long, written to disk for the other languages.
draws <- do.call(rbind, lapply(1:200, function(r) {
P <- draw_bounded_panel()
data.frame(rep = r, unit = rep(1:(Jp + 1), each = TTp),
time = rep(1:TTp, times = Jp + 1), share = as.vector(t(P)))
}))
write.csv(draws, file.path(dir_sim, "sim_prop_draws.csv"), row.names = FALSE)Jp <- 10; T0p <- 18; TTp <- 30; pre_p <- 1:T0p
# Convex SCM weights: gradient steps, clip negatives to zero, renormalize to sum to one.
scm_weights <- function(D, y, pre, iter = 30000, lr = 5e-3) {
Jn <- ncol(D); w <- rep(1 / Jn, Jn)
Dp <- D[pre, , drop = FALSE]; yp <- y[pre]
for (k in 1:iter) {
w <- w - lr * as.vector(t(Dp) %*% (Dp %*% w - yp))
w[w < 0] <- 0
s <- sum(w); if (s > 0) w <- w / s
}
w
}
draws <- read.csv("synthetic-control_files/data/sim_prop_draws.csv")
reps <- t(sapply(1:200, function(r) {
sub <- draws[draws$rep == r, ]
P <- matrix(sub$share, nrow = Jp + 1, byrow = TRUE)
D <- t(P[-1, ]); y <- P[1, ]
uc <- lm(y[pre_p] ~ D[pre_p, ])
cf_uc <- as.vector(cbind(1, D) %*% ifelse(is.na(coef(uc)), 0, coef(uc)))
cf_scm <- as.vector(D %*% scm_weights(D, y, pre_p))
c(min_uc = min(cf_uc), max_uc = max(cf_uc),
min_scm = min(cf_scm), max_scm = max(cf_scm),
rmse_uc = sqrt(mean((cf_uc[pre_p] - y[pre_p])^2)),
rmse_scm = sqrt(mean((cf_scm[pre_p] - y[pre_p])^2)))
}))
data.frame(
Estimator = c("Convex SCM", "Unconstrained"),
`Median pre-period RMSE` = round(c(median(reps[, "rmse_scm"]),
median(reps[, "rmse_uc"])), 5),
`% draws with counterfactual < 0` = round(c(mean(reps[, "min_scm"] < 0),
mean(reps[, "min_uc"] < 0)) * 100, 1),
`% draws with counterfactual > 1` = round(c(mean(reps[, "max_scm"] > 1),
mean(reps[, "max_uc"] > 1)) * 100, 1),
`Worst value` = round(c(min(reps[, "min_scm"]), min(reps[, "min_uc"])), 3),
check.names = FALSE
) Estimator Median pre-period RMSE % draws with counterfactual < 0
1 Convex SCM 0.00353 0.0
2 Unconstrained 0.00126 21.5
% draws with counterfactual > 1 Worst value
1 0.0 0.016
2 2.5 -0.576Jp, T0p, TTp = 10, 18, 30
pre = np.arange(T0p)
def scm_weights(D, y, iters=30000, lr=5e-3):
w = np.full(D.shape[1], 1 / D.shape[1])
Dp, yp = D[:T0p], y[:T0p]
for _ in range(iters):
w = w - lr * (Dp.T @ (Dp @ w - yp))
w[w < 0] = 0
s = w.sum()
if s > 0:
w /= s
return w
draws = pd.read_csv("synthetic-control_files/data/sim_prop_draws.csv")
rec = []
for r in range(1, 201):
P = (draws[draws.rep == r]
.pivot(index="unit", columns="time", values="share").values) # (Jp+1, TTp)
D, y = P[1:].T, P[0]
X = np.column_stack([np.ones(T0p), D[:T0p]])
beta, *_ = np.linalg.lstsq(X, y[:T0p], rcond=None)
cf_uc = np.column_stack([np.ones(TTp), D]) @ beta
cf_scm = D @ scm_weights(D, y)
rec.append((cf_uc.min(), cf_uc.max(), cf_scm.min(), cf_scm.max(),
np.sqrt(((cf_uc[:T0p] - y[:T0p]) ** 2).mean()),
np.sqrt(((cf_scm[:T0p] - y[:T0p]) ** 2).mean())))
reps = np.array(rec)
print(pd.DataFrame({
"Estimator": ["Convex SCM", "Unconstrained"],
"Median pre-period RMSE": [round(np.median(reps[:, 5]), 5), round(np.median(reps[:, 4]), 5)],
"% cf < 0": [round((reps[:, 2] < 0).mean() * 100, 1), round((reps[:, 0] < 0).mean() * 100, 1)],
"% cf > 1": [round((reps[:, 3] > 1).mean() * 100, 1), round((reps[:, 1] > 1).mean() * 100, 1)],
"Worst": [round(reps[:, 2].min(), 3), round(reps[:, 0].min(), 3)],
}).to_string(index=False)) Estimator Median pre-period RMSE % cf < 0 % cf > 1 Worst
Convex SCM 0.00353 0.0 0.0 0.016
Unconstrained 0.00126 21.5 2.5 -0.576import delimited "synthetic-control_files/data/sim_prop_draws.csv", clear
mata:
// convex SCM weights: gradient step, clip negatives, renormalize (simplex heuristic)
real colvector scm_w(real matrix Dm, real colvector y, real scalar T0,
real scalar iters, real scalar lr) {
real scalar J, k, s
real colvector w
real matrix Dp
Dp = Dm[1..T0,]; J = cols(Dm); w = J(J,1,1/J)
for (k=1; k<=iters; k++) {
w = w :- lr * (Dp' * (Dp*w :- y[1..T0]))
w = w :* (w :> 0)
s = colsum(w); if (s>0) w = w / s
}
return(w)
}
real scalar med(real colvector v) {
real colvector s
s = sort(v,1); return((s[100]+s[101])/2)
}
st_view(REP=.,.,"rep"); st_view(UN=.,.,"unit"); st_view(TM=.,.,"time"); st_view(SH=.,.,"share")
T0 = 18
minuc=J(200,1,.); maxuc=J(200,1,.); minsc=J(200,1,.); maxsc=J(200,1,.)
rmuc=J(200,1,.); rmsc=J(200,1,.)
for (r=1; r<=200; r++) {
idx = selectindex(REP:==r)
un = UN[idx]; tm = TM[idx]; sh = SH[idx]
P = J(11,30,.)
for (i=1;i<=rows(idx);i++) P[un[i], tm[i]] = sh[i]
y = P[1,]'; Dm = P[2..11,]'
X = (J(T0,1,1), Dm[1..T0,]); b = invsym(X'X)*X'y[1..T0]
cf_uc = (J(30,1,1), Dm)*b
cf_sc = Dm*scm_w(Dm, y, T0, 3000, 5e-3)
minuc[r]=colmin(cf_uc); maxuc[r]=colmax(cf_uc)
minsc[r]=colmin(cf_sc); maxsc[r]=colmax(cf_sc)
rmuc[r]=sqrt(mean((cf_uc[1..T0]-y[1..T0]):^2))
rmsc[r]=sqrt(mean((cf_sc[1..T0]-y[1..T0]):^2))
}
printf("{txt}Estimator Median pre-RMSE %% cf<0 %% cf>1 Worst\n")
printf("Convex SCM %9.5f %5.1f %5.1f %6.3f\n",
med(rmsc), mean(minsc:<0)*100, mean(maxsc:>1)*100, colmin(minsc))
printf("Unconstrained %9.5f %5.1f %5.1f %6.3f\n",
med(rmuc), mean(minuc:<0)*100, mean(maxuc:>1)*100, colmin(minuc))
end(encoding automatically selected: ISO-8859-2)
(4 vars, 66,000 obs)
------------------------------------------------- mata (type end to exit) ------------------------------------------------------------------
: // convex SCM weights: gradient step, clip negatives, renormalize (simplex heuristic)
: real colvector scm_w(real matrix Dm, real colvector y, real scalar T0,
> real scalar iters, real scalar lr) {
> real scalar J, k, s
> real colvector w
> real matrix Dp
> Dp = Dm[1..T0,]; J = cols(Dm); w = J(J,1,1/J)
> for (k=1; k<=iters; k++) {
> w = w :- lr * (Dp' * (Dp*w :- y[1..T0]))
> w = w :* (w :> 0)
> s = colsum(w); if (s>0) w = w / s
> }
> return(w)
> }
: real scalar med(real colvector v) {
> real colvector s
> s = sort(v,1); return((s[100]+s[101])/2)
> }
: st_view(REP=.,.,"rep"); st_view(UN=.,.,"unit"); st_view(TM=.,.,"time"); st_view(SH=.,.,"share")
: T0 = 18
: minuc=J(200,1,.); maxuc=J(200,1,.); minsc=J(200,1,.); maxsc=J(200,1,.)
: rmuc=J(200,1,.); rmsc=J(200,1,.)
: for (r=1; r<=200; r++) {
> idx = selectindex(REP:==r)
> un = UN[idx]; tm = TM[idx]; sh = SH[idx]
> P = J(11,30,.)
> for (i=1;i<=rows(idx);i++) P[un[i], tm[i]] = sh[i]
> y = P[1,]'; Dm = P[2..11,]'
> X = (J(T0,1,1), Dm[1..T0,]); b = invsym(X'X)*X'y[1..T0]
> cf_uc = (J(30,1,1), Dm)*b
> cf_sc = Dm*scm_w(Dm, y, T0, 3000, 5e-3)
> minuc[r]=colmin(cf_uc); maxuc[r]=colmax(cf_uc)
> minsc[r]=colmin(cf_sc); maxsc[r]=colmax(cf_sc)
> rmuc[r]=sqrt(mean((cf_uc[1..T0]-y[1..T0]):^2))
> rmsc[r]=sqrt(mean((cf_sc[1..T0]-y[1..T0]):^2))
> }
: printf("{txt}Estimator Median pre-RMSE %% cf<0 %% cf>1 Worst\n")
Estimator Median pre-RMSE % cf<0 % cf>1 Worst
: printf("Convex SCM %9.5f %5.1f %5.1f %6.3f\n",
> med(rmsc), mean(minsc:<0)*100, mean(maxsc:>1)*100, colmin(minsc))
Convex SCM 0.00383 0.0 0.0 0.016
: printf("Unconstrained %9.5f %5.1f %5.1f %6.3f\n",
> med(rmuc), mean(minuc:<0)*100, mean(maxuc:>1)*100, colmin(minuc))
Unconstrained 0.00126 21.5 2.5 -0.576
: end
--------------------------------------------------------------------------------------------------------------------------------------------The unconstrained estimator fits the pre-treatment period better than convex SCM, and in a meaningful share of draws it produces a counterfactual share below zero or above one. The convex estimator never does, in any draw. A predicted deforestation share of \(-0.5\) is not a surprising finding; it is a scale problem announcing itself.
# First saved draw whose unconstrained counterfactual leaves [0, 1].
ex <- NULL
for (r in 1:200) {
sub <- draws[draws$rep == r, ]
P <- matrix(sub$share, nrow = Jp + 1, byrow = TRUE)
D <- t(P[-1, ]); y <- P[1, ]
uc <- lm(y[pre_p] ~ D[pre_p, ])
cf_uc <- as.vector(cbind(1, D) %*% ifelse(is.na(coef(uc)), 0, coef(uc)))
if (min(cf_uc) < 0) {
ex <- list(y = y, uc = cf_uc, scm = as.vector(D %*% scm_weights(D, y, pre_p)))
break
}
}
data.frame(time = 1:TTp, Observed = ex$y,
`Convex SCM` = ex$scm, `Unconstrained` = ex$uc,
check.names = FALSE) |>
pivot_longer(-time, names_to = "series", values_to = "share") |>
ggplot(aes(time, share, colour = series, linetype = series)) +
geom_hline(yintercept = 0, colour = "grey30") +
geom_vline(xintercept = T0p, linetype = "dotted", colour = "grey40") +
geom_line(linewidth = 0.9) +
scale_colour_manual(values = c(Observed = "grey30",
`Convex SCM` = PAL[["synth"]],
Unconstrained = PAL[["treated"]])) +
labs(x = "Period", y = "Proportion", colour = NULL, linetype = NULL,
subtitle = "Solid line at zero is the lower bound for a share")
ex = None
for r in range(1, 201):
P = (draws[draws.rep == r]
.pivot(index="unit", columns="time", values="share").values)
D, y = P[1:].T, P[0]
X = np.column_stack([np.ones(T0p), D[:T0p]])
beta, *_ = np.linalg.lstsq(X, y[:T0p], rcond=None)
cf_uc = np.column_stack([np.ones(TTp), D]) @ beta
if cf_uc.min() < 0:
ex = (y, D @ scm_weights(D, y), cf_uc)
break
y, cf_scm, cf_uc = ex
t = range(1, TTp + 1)
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.axhline(0, color="0.3")<matplotlib.lines.Line2D object at 0x16c5cad80>ax.axvline(T0p, color="grey", linestyle=":")<matplotlib.lines.Line2D object at 0x16c5ca540>ax.plot(t, y, color="0.3", label="Observed")[<matplotlib.lines.Line2D object at 0x16c5679e0>]ax.plot(t, cf_scm, color="#1F4E79", label="Convex SCM")[<matplotlib.lines.Line2D object at 0x16c7a2ae0>]ax.plot(t, cf_uc, color="#B22222", label="Unconstrained")[<matplotlib.lines.Line2D object at 0x16c7a0d10>]ax.set_xlabel("Period")Text(0.5, 0, 'Period')ax.set_ylabel("Proportion")Text(0, 0.5, 'Proportion')ax.legend(frameon=False)<matplotlib.legend.Legend object at 0x16c4466f0>plt.tight_layout()
plt.show()
import delimited "synthetic-control_files/data/sim_prop_draws.csv", clear
mata:
real colvector scm_w(real matrix Dm, real colvector y, real scalar T0,
real scalar iters, real scalar lr) {
real scalar J, k, s
real colvector w
real matrix Dp
Dp = Dm[1..T0,]; J = cols(Dm); w = J(J,1,1/J)
for (k=1; k<=iters; k++) {
w = w :- lr * (Dp' * (Dp*w :- y[1..T0]))
w = w :* (w :> 0)
s = colsum(w); if (s>0) w = w / s
}
return(w)
}
st_view(REP=.,.,"rep"); st_view(UN=.,.,"unit"); st_view(TM=.,.,"time"); st_view(SH=.,.,"share")
T0 = 18; res = J(30,4,.)
for (r=1; r<=200; r++) {
idx = selectindex(REP:==r)
un = UN[idx]; tm = TM[idx]; sh = SH[idx]
P = J(11,30,.)
for (i=1;i<=rows(idx);i++) P[un[i], tm[i]] = sh[i]
y = P[1,]'; Dm = P[2..11,]'
X = (J(T0,1,1), Dm[1..T0,]); b = invsym(X'X)*X'y[1..T0]
cf_uc = (J(30,1,1), Dm)*b
if (colmin(cf_uc) < 0) {
cf_sc = Dm*scm_w(Dm, y, T0, 3000, 5e-3)
res = ((1::30), y, cf_sc, cf_uc)
break
}
}
st_matrix("M", res)
end
clear
svmat M
rename (M1 M2 M3 M4) (time observed convex uncon)
twoway (line observed time, lcolor(gs6)) ///
(line convex time, lcolor("31 78 121")) ///
(line uncon time, lcolor("178 34 34")), ///
yline(0) xline(18, lpattern(dot)) ///
legend(order(1 "Observed" 2 "Convex SCM" 3 "Unconstrained")) ///
xtitle("Period") ytitle("Proportion")
quietly graph export "sc_figs/fig-prop-stata.svg", replace width(1600)(encoding automatically selected: ISO-8859-2)
(4 vars, 66,000 obs)
------------------------------------------------- mata (type end to exit) ------------------------------------------------------------------
: real colvector scm_w(real matrix Dm, real colvector y, real scalar T0,
> real scalar iters, real scalar lr) {
> real scalar J, k, s
> real colvector w
> real matrix Dp
> Dp = Dm[1..T0,]; J = cols(Dm); w = J(J,1,1/J)
> for (k=1; k<=iters; k++) {
> w = w :- lr * (Dp' * (Dp*w :- y[1..T0]))
> w = w :* (w :> 0)
> s = colsum(w); if (s>0) w = w / s
> }
> return(w)
> }
: st_view(REP=.,.,"rep"); st_view(UN=.,.,"unit"); st_view(TM=.,.,"time"); st_view(SH=.,.,"share")
: T0 = 18; res = J(30,4,.)
: for (r=1; r<=200; r++) {
> idx = selectindex(REP:==r)
> un = UN[idx]; tm = TM[idx]; sh = SH[idx]
> P = J(11,30,.)
> for (i=1;i<=rows(idx);i++) P[un[i], tm[i]] = sh[i]
> y = P[1,]'; Dm = P[2..11,]'
> X = (J(T0,1,1), Dm[1..T0,]); b = invsym(X'X)*X'y[1..T0]
> cf_uc = (J(30,1,1), Dm)*b
> if (colmin(cf_uc) < 0) {
> cf_sc = Dm*scm_w(Dm, y, T0, 3000, 5e-3)
> res = ((1::30), y, cf_sc, cf_uc)
> break
> }
> }
: st_matrix("M", res)
: end
--------------------------------------------------------------------------------------------------------------------------------------------
number of observations will be reset to 30
Press any key to continue, or Break to abort
Number of observations (_N) was 0, now 30.
Choose the Scale Deliberately
| Approach | When it is useful | Main caution |
|---|---|---|
| Raw proportion | Effects are naturally percentage points and the fit is sensible | Denominators may imply very different precision across observations |
| Logit-transformed proportion | Effects are more naturally multiplicative, or values crowd near 0 and 1 | Zeros and ones need a documented adjustment; back-transformation is nonlinear |
| Count or rate framing | You observe numerators and denominators and sampling precision matters | The estimator and interpretation may need to change, not just the outcome column |
The raw scale is usually the right place to start because it is interpretable: a gap of \(0.04\) is four percentage points. Transform when it changes the fit, the modeling problem, or the substantive interpretation in a useful way, not because a textbook made proportions sound dangerous (Papke and Wooldridge 1996; Warton and Hui 2011).
Denominators Are Part of the Data-Generating Process
A turnout rate based on 200 voters is not measured as precisely as one based on two million. A deforestation share from a tiny polygon is not the same object as one from a large landscape. If denominators vary sharply across units or time, say so and consider whether a count- or rate-aware analysis is warranted.
Report the numerator, denominator, range of denominators, occurrence of zeros and ones, outcome scale, and how any transformed estimate translates back into substantive units.
A Practitioner’s Checklist
Work through this in order. Items one through four happen before you look at any post-treatment result.
Design, before estimation
- Define the treatment precisely. What changed, for whom, and when did behavior plausibly begin to respond? Anticipation moves the effective date earlier than the statute.
- Fix the treatment date and justify it. Enactment, implementation, and crediting dates are different dates.
- Build the donor pool from substantive knowledge. Exclude units exposed to the same or similar interventions, units subject to spillover or displacement, and units whose measurement changed. Write the exclusion rule down before estimating.
- Fix the predictor set. Decide what enters \(X\) and on what pre-treatment window, aware that the full outcome path can crowd out covariates (Kaul et al. 2022).
Fit and diagnose
- Check pre-treatment fit first. Plot it. Report pre-treatment RMSPE on the outcome’s own scale. If the paths do not align, stop and diagnose rather than proceeding.
- Inspect the donor weights. Name the units carrying the comparison and say whether that composition is substantively sensible.
- Report predictor balance against the raw donor average.
- Confirm the treated unit is inside donor support. If it is not, say so, and consider augmentation (Ben-Michael, Feller, and Rothstein 2021) rather than pretending the convex fit succeeded.
- Check the outcome scale. If the outcome is bounded, verify that fitted counterfactuals stay in range, especially with estimators that permit negative weights.
Inference and robustness
- Run placebo tests, state the test statistic, state the pre-fit screening rule, and state whether the comparison is one- or two-sided.
- Consider design-appropriate intervals rather than placebo p-values alone, such as conformal inference (Chernozhukov, Wüthrich, and Zhu 2021) or prediction intervals (Cattaneo, Feng, and Titiunik 2021).
- Run leave-one-out checks on every high-weight donor.
- Vary the donor pool across defensible exclusion rules and report whether conclusions move.
Reporting
- Show the path plot and the gap plot. These are the evidence, not illustrations of it.
- Report the estimate on an interpretable scale, in the outcome’s own units.
- Record software, package versions, optimization settings, and a random seed. Synthetic control weights are not always unique (Klößner et al. 2018).
- Name the most credible failure mode in the limitations. If a referee could break the design one way, say which way.
The goal is not to make the synthetic control look inevitable. It is to show that it is the most credible available counterfactual, while being candid about where that claim could break.
FAQ
When should I use SCM instead of difference-in-differences?
Use SCM when one or a few important units receive treatment and no single comparison unit is credible. Use DiD when many treated and untreated units support a defensible parallel-trends argument. When both are plausible, show both; agreement is reassuring and disagreement is information.
How many pre-treatment periods do I need?
Enough to assess persistent fit. There is no universal threshold, because it depends on outcome volatility, the number of latent factors, and donor-pool size. Two or three periods can fit perfectly by accident, which is overfitting rather than validation.
What if pre-treatment fit is poor?
Diagnose before estimating effects. Reconsider the donor pool, outcome scale, treatment date, and predictor set, and check whether the treated unit lies outside donor support. Augmented SCM or a factor-model estimator may help. If fit remains poor, reporting that the design is not feasible is a legitimate result, as the poor-support simulation above shows.
Can I include covariates?
Yes, using pre-treatment covariates that predict the untreated outcome. Do not include post-treatment variables that treatment could have affected, and do not expect covariates to rescue a synthetic control that misses the pre-treatment outcome path.
What if donors are affected by spillovers?
Exclude them, define a less exposed donor pool, or show sensitivity analyses without them. A partly treated donor cannot cleanly reveal the untreated potential outcome.
What if treatment timing varies across units?
Classical Synth is not the right first tool. Use a counterfactual-imputation estimator such as fect(..., method = "gsynth"), and see the staggered-adoption literature for estimand choices (Ben-Michael, Feller, and Rothstein 2022).
How do I interpret a placebo p-value?
As a design-based rank comparison, not a conventional test. It says how unusual the treated unit’s gap statistic is among eligible fake-treated donors, and it is meaningful only when placebo units have reasonable pre-treatment fit and the donor pool is large enough.
Cross-language translation
This appendix records the functional equivalents each tab relies on. Where an estimator’s optimizer differs across implementations, the estimand agrees but the exact weights or standard errors need not (see the Proposition 99 note); where a step is deterministic, all three agree to floating-point precision on the shared data.
| Task | R | Python | Stata |
|---|---|---|---|
| Data wrangling | dplyr/tidyr, read.csv |
pandas |
import delimited, native |
| Plain regression / linear algebra | lm, matrix ops |
numpy, statsmodels |
regress, Mata |
| Classical SCM | Synth::synth |
pysyncon.Synth |
synth (Abadie) |
| Augmented SCM | augsynth::augsynth |
pysyncon.AugSynth |
allsynth (Wiltshire) |
| Generalized SCM / factor imputation | fect::fect |
fect_python † |
fect (Liu–Wang–Xu) |
| Figures | ggplot2 |
matplotlib |
twoway |
| Simulated data | R generators (folded “Simulation design”) write CSVs | read the shared CSVs | read the shared CSVs |
† fect_python is an rpy2 wrapper around R’s fect; its import currently conflicts with this environment’s rpy2/ggplot2 versions, so the fect Python tabs are shown with eval: false (the code is the intended call). Because random-number streams differ across languages, every simulated panel is generated once in R, written to synthetic-control_files/data/, and read identically by all three languages, so results match downstream rather than depending on a shared RNG.
Reproducibility
All analyses were run in R with the packages attached in the setup chunk at the top of this document. Three points need a note.
Installing augsynth. The augmented synthetic control estimator is not on CRAN; it is distributed on GitHub (Ben-Michael, Feller, and Rothstein 2021). Install it once before rendering, using either of:
pacman::p_load_gh("ebenmichael/augsynth") # via pacman
remotes::install_github("ebenmichael/augsynth") # or via remotesThe setup chunk then only attaches it, with no network call, so the render stays deterministic. augsynth suggests MCPanel (also GitHub-only) for its matrix-completion option; this guide does not use that option, so a failed MCPanel build is harmless.
Stata and Python engines. Statamarkdown registers a stata knitr engine, which requires a local Stata installation and usually an explicit engine.path. reticulate provides the python engine; point it at an interpreter with reticulate::use_python() or use_condaenv(). Because a global RETICULATE_PYTHON in ~/.Rprofile would override the guide’s own environment, the setup chunk clears it and binds the local .venv.
Stata packages. The Stata tabs use several user-written commands. Install them once with:
ssc install synth // classical SCM (Abadie, Diamond, Hainmueller)
ssc install allsynth // augmented / bias-corrected SCM (Wiltshire)
ssc install distinct // allsynth dependency
ssc install elasticregress // allsynth ridge/lasso dependency
net install fect, from("https://raw.githubusercontent.com/xuyiqing/fect_stata/master/")
ssc install reghdfe // fect dependencies
ssc install ftools
ssc install _gwtmean
ssc install moremata
ssc install requireStatamarkdown does not auto-embed Stata graphs, so each Stata figure block ends with graph export to sc_figs/ followed by an image include. Intermediate .dta handoffs between Stata blocks are written to synthetic-control_files/stata/ (created by the setup chunk) rather than the project root.
Python packages. The .venv holds numpy, pandas, scipy, statsmodels, matplotlib, and pysyncon. The generalized-SCM (fect) Python tabs use a compact native NumPy implementation rather than fect_python, whose rpy2 bridge conflicts with this environment.
Data. The datasets are cached as CSVs under synthetic-control_files/data/, exported from tidysynth and fect, with full provenance in that folder’s PROVENANCE.md and in the footnotes where each dataset is introduced. Simulated panels are generated once in R (the folded “Simulation design” blocks) and written there so all three languages read identical inputs.
The exact package versions used to render this document:
Session info
if (requireNamespace("sessioninfo", quietly = TRUE)) {
sessioninfo::session_info()
} else {
sessionInfo()
}R version 4.6.1 (2026-06-24)
Platform: aarch64-apple-darwin23
Running under: macOS Sequoia 15.7.7
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US/C/en_US/C/en_US/en_US
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] augsynth_0.2.0 reticulate_1.46.0 Statamarkdown_0.9.6
[4] ggplot2_4.0.3 tidyr_1.3.2 dplyr_1.2.1
[7] tidysynth_0.2.1 fect_2.4.5 Synth_1.1-10
loaded via a namespace (and not attached):
[1] gtable_0.3.6 xfun_0.60 htmlwidgets_1.6.4
[4] GGally_2.4.0 lattice_0.22-9 numDeriv_2016.8-1.1
[7] vctrs_0.7.3 tools_4.6.1 generics_0.1.4
[10] parallel_4.6.1 sandwich_3.1-2 tibble_3.3.1
[13] pacman_0.5.1 pkgconfig_2.0.3 Matrix_1.7-5
[16] RColorBrewer_1.1-3 S7_0.2.2 rngtools_1.5.2
[19] stringmagic_1.2.0 optimx_2025-4.9 lifecycle_1.0.5
[22] stringr_1.6.0 compiler_4.6.1 farver_2.1.2
[25] codetools_0.2-20 htmltools_0.5.9 yaml_2.3.12
[28] Formula_1.2-5 pracma_2.4.6 pillar_1.11.1
[31] nloptr_2.2.1 MASS_7.3-65 doRNG_1.8.6.3
[34] iterators_1.0.14 rgenoud_5.9-0.11 abind_1.4-8
[37] foreach_1.5.2 nlme_3.1-169 parallelly_1.48.0
[40] ggstats_0.13.0 tidyselect_1.2.1 digest_0.6.39
[43] stringi_1.8.7 mvtnorm_1.4-2 future_1.75.0
[46] reshape2_1.4.5 kernlab_0.9-33 purrr_1.2.2
[49] listenv_1.0.0 labeling_0.4.3 fastmap_1.2.0
[52] grid_4.6.1 cli_3.6.6 magrittr_2.0.5
[55] future.apply_1.20.2 withr_3.0.3 dreamerr_1.5.0
[58] scales_1.4.0 estimability_2.0.0 rmarkdown_2.31
[61] emmeans_2.0.4 globals_0.19.1 otel_0.2.0
[64] gridExtra_2.3.1 png_0.1-9 zoo_1.8-15
[67] evaluate_1.0.5 knitr_1.51 fixest_0.14.2
[70] rlang_1.3.0 Rcpp_1.1.2 xtable_1.8-8
[73] glue_1.8.1 osqp_1.0.0 jsonlite_2.0.0
[76] plyr_1.8.9 R6_2.6.1 Code
import sys, numpy, pandas, scipy, matplotlib, pysyncon
print("Python", sys.version.split()[0])Python 3.12.2Code
for m in (numpy, pandas, scipy, matplotlib, pysyncon):
print(f"{m.__name__:12s} {m.__version__}")numpy 2.5.2
pandas 3.0.5
scipy 1.18.0
matplotlib 3.11.1
pysyncon 1.5.1aboutStataNow/MP 19.5 for Mac (Apple Silicon)
Revision 03 Jun 2026
Copyright 1985-2025 StataCorp LLC
Total physical memory: 32.01 GB
Stata license: Single-user 2-core , expiring 1 Jul 2027
Serial number: 501909374535
Licensed to: Ellsworth, Ted
Yale UniversityReferences
Abadie, Alberto. 2021. “Using Synthetic Controls: Feasibility, Data Requirements, and Methodological Aspects.” Journal of Economic Literature 59 (2): 391–425. https://doi.org/10.1257/jel.20191450.
Abadie, Alberto, Alexis Diamond, and Jens Hainmueller. 2010. “Synthetic Control Methods for Comparative Case Studies: Estimating the Effect of California’s Tobacco Control Program.” Journal of the American Statistical Association 105 (490): 493–505. https://doi.org/10.1198/jasa.2009.ap08746.
———. 2015. “Comparative Politics and the Synthetic Control Method.” American Journal of Political Science 59 (2): 495–510. https://doi.org/10.1111/ajps.12116.
Abadie, Alberto, and Javier Gardeazabal. 2003. “The Economic Costs of Conflict: A Case Study of the Basque Country.” American Economic Review 93 (1): 113–32. https://doi.org/10.1257/000282803321455188.
Ben-Michael, Eli, Avi Feller, and Jesse Rothstein. 2021. “The Augmented Synthetic Control Method.” Journal of the American Statistical Association 116 (536): 1789–1803. https://doi.org/10.1080/01621459.2021.1929245.
———. 2022. “Synthetic Controls with Staggered Adoption.” Journal of the Royal Statistical Society Series B: Statistical Methodology 84 (2): 351–81. https://doi.org/10.1111/rssb.12448.
Born, Benjamin, Gernot J. Müller, Moritz Schularick, and Petr Sedláček. 2019. “The Costs of Economic Nationalism: Evidence from the Brexit Experiment.” The Economic Journal 129 (623): 2722–44. https://doi.org/10.1093/ej/uez020.
Cattaneo, Matias D., Yingjie Feng, and Rocio Titiunik. 2021. “Prediction Intervals for Synthetic Control Methods.” Journal of the American Statistical Association 116 (536): 1865–80. https://doi.org/10.1080/01621459.2021.1979561.
Chernozhukov, Victor, Kaspar Wüthrich, and Yinchu Zhu. 2021. “An Exact and Robust Conformal Inference Method for Counterfactual and Synthetic Controls.” Journal of the American Statistical Association 116 (536): 1849–64. https://doi.org/10.1080/01621459.2021.1920957.
Dunford, Eric. 2025. tidysynth: A Tidy Implementation of the Synthetic Control Method (R package version 0.2.1). Comprehensive R Archive Network. https://doi.org/10.32614/CRAN.package.tidysynth.
Ferman, Bruno, Cristine Pinto, and Vitor Possebom. 2020. “Cherry Picking with Synthetic Controls.” Journal of Policy Analysis and Management 39 (2): 510–32. https://doi.org/10.1002/pam.22206.
Hainmueller, Jens, and Alexis Diamond. 2026. Synth: Synthetic Control Group Method for Comparative Case Studies (R package version 1.1-10). Comprehensive R Archive Network. https://doi.org/10.32614/CRAN.package.Synth.
Kaul, Ashok, Stefan Klößner, Gregor Pfeifer, and Manuel Schieler. 2022. “Standard Synthetic Control Methods: The Case of Using All Preintervention Outcomes Together with Covariates.” Journal of Business & Economic Statistics 40 (3): 1362–76. https://doi.org/10.1080/07350015.2021.1930012.
Klößner, Stefan, Ashok Kaul, Gregor Pfeifer, and Manuel Schieler. 2018. “Comparative Politics and the Synthetic Control Method Revisited: A Note on Abadie et al. (2015).” Swiss Journal of Economics and Statistics 154: 11. https://doi.org/10.1186/s41937-017-0004-9.
Papke, Leslie E., and Jeffrey M. Wooldridge. 1996. “Econometric Methods for Fractional Response Variables with an Application to 401(k) Plan Participation Rates.” Journal of Applied Econometrics 11 (6): 619–32. https://doi.org/10.1002/(SICI)1099-1255(199611)11:6<619::AID-JAE418>3.0.CO;2-1.
Warton, David I., and Francis K. C. Hui. 2011. “The Arcsine Is Asinine: The Analysis of Proportions in Ecology.” Ecology 92 (1): 3–10. https://doi.org/10.1890/10-0340.1.
West, Thales A. P., Sven Wunder, Erin O. Sills, Jan Börner, Sami W. Rifai, Alexandra N. Neidermeier, Gabriel P. Frey, and Andreas Kontoleon. 2023. “Action Needed to Make Carbon Offsets from Forest Conservation Work for Climate Change Mitigation.” Science 381 (6660): 873–77. https://doi.org/10.1126/science.ade3535.
Xu, Yiqing. 2017. “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models.” Political Analysis 25 (1): 57–76. https://doi.org/10.1017/pan.2016.2.
Xu, Yiqing, Licheng Liu, Ye Wang, Ziyi Liu, Shijian Liu, Tianzhu Qin, Jinwen Wu, and Rivka Lipkovitz. 2026. fect: Fixed Effects Counterfactual Estimators – User Manual (v2.4.5). https://yiqingxu.org/packages/fect/.
Footnotes
This guide reads a cached copy,
synthetic-control_files/data/smoking.csv, exported fromtidysynth::smoking(v0.2.1); seePROVENANCE.mdin that folder for full provenance.↩︎Read here as a cached copy,
synthetic-control_files/data/simgsynth.csv, exported fromfect::simgsynth(v2.4.5); seePROVENANCE.md.↩︎Read here as a cached copy,
synthetic-control_files/data/turnout.csv, exported fromfect::turnout(v2.4.5); seePROVENANCE.md.↩︎
Citation
BibTeX citation:
@online{elkobi2026,
author = {Elkobi, Jonathan},
publisher = {Yale Library StatLab},
title = {Synthetic {Control} {Methods}},
date = {2026-07-16},
langid = {en},
abstract = {Synthetic control methods estimate what would have
happened to a treated unit without treatment by constructing a
weighted comparison from untreated units. This guide develops the
method around a single running example, California’s Proposition 99,
and uses simulated panels with known answers to show what good and
bad designs actually look like. It covers the design logic,
diagnostics, a reproducible Synth workflow, generalized synthetic
control in fect, augmented synthetic controls, and bounded or
proportion outcomes.}
}
For attribution, please cite this work as:
Elkobi, Jonathan. 2026. “Synthetic Control Methods.” Yale
Library StatLab. July 16, 2026.