import matplotlib
matplotlib.use("Agg") # headless backend: render to buffers, never to a screen
import numpy as np
import pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt
# pysyncon is the reference synthetic-control package; we cross-check it against a hand-rolled solve.
from pysyncon import Dataprep, Synth
rng = np.random.default_rng(20260528) # one pinned generator drives the entire notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
print("pandas", pd.__version__, "| numpy", np.__version__)nb4.4 — Synthetic Control & Synthetic DiD: Building the Counterfactual You Wish You Had
Week 4, Chapter 4.4. Companion notebook.
In Chapter 4.4 Priya was handed the hardest version of a causal question: one treated unit. A single state passed a first-in-the-nation climate-risk disclosure mandate in 2018, forcing home insurers to publish wildfire- and flood-exposure pricing. Did the mandate move average homeowner premiums? With one treated state and no clean twin, difference-in-differences’ parallel-trends bet collapses — there is no naturally occurring control group that tracks one specific unit.
The fix, due to Abadie, Diamond and Hainmueller (2010), is to stop picking a control and start building one: a synthetic control, a convex-weighted blend of untreated donor states chosen so the blend matches the treated state’s premium path before the mandate. Then you watch the real state and its synthetic twin diverge after 2018, and read the gap as the effect.
This notebook does the whole study end to end on a simulated premium panel whose truth we control:
- Simulate a 31-unit panel (1 treated + 30 donors) from a two-factor model with a known +1.4 post-treatment jump injected into the treated state.
- Fit synthetic control two ways — with the
pysynconpackage and with a hand-rolledscipy.optimizeconvex solve — and confirm they agree. Plot treated-vs-synthetic and the gap. - Placebo / permutation inference: re-run the method treating each donor as fake-treated, compute the post/pre RMSPE ratio for every unit, and get the permutation p-value = \text{rank}/(J+1).
- A no-effect sanity check: delete the injected jump and confirm the treated unit falls back into the placebo crowd.
- Synthetic DiD (Arkhangelsky et al., 2021): unit and time weights with DiD-style differencing, compared against plain synthetic control and plain DiD.
We will hit the chapter’s headline numbers almost exactly: pre-RMSPE \approx 0.035, post-RMSPE \approx 1.5, ratio \approx 44, estimated effect \approx 1.5 (true 1.4), treated unit ranks 1 of 31 with permutation p \approx 0.032.
2. Synthetic control, hand-rolled: the convex solve
The estimator is one optimization. Collect the treated state’s pre-period premiums into a vector \mathbf{X}_1 and the donors’ pre-period premiums into a matrix \mathbf{X}_0, then choose weights \mathbf{w} to minimize the pre-period mismatch, subject to the two constraints that are the method:
\mathbf{w}^{*} = \arg\min_{\mathbf{w}} \;\sum_{t \le T_0}\Big(Y_{1t} - \textstyle\sum_j w_j Y_{jt}\Big)^2 \quad\text{s.t.}\quad w_j \ge 0,\;\; \sum_j w_j = 1.
The bound w_j \ge 0 and the equality \sum_j w_j = 1 make the synthetic control a convex combination — a genuine weighted average that cannot extrapolate past the donors and cannot fake a fit with wild offsetting coefficients. We solve it with scipy.optimize.minimize using SLSQP (it handles the equality constraint and the box bounds together). We take \mathbf{V} = \mathbf{I}, weighting every pre-year equally.
def synth_weights(treated_pre, donors_pre):
'''Convex weights minimizing pre-period sum of squared mismatch (V = I).
treated_pre : (T0,) vector of the treated unit's pre-period outcomes
donors_pre : (Jn, T0) matrix of donor pre-period outcomes (one row per donor)
'''
Jn = donors_pre.shape[0]
def loss(w):
return np.sum((treated_pre - w @ donors_pre) ** 2)
cons = ({"type": "eq", "fun": lambda w: w.sum() - 1.0},) # weights sum to 1
bnds = [(0.0, 1.0)] * Jn # weights >= 0 (convexity)
w0 = np.full(Jn, 1.0 / Jn) # start at equal weights (= the DiD control)
res = minimize(loss, w0, method="SLSQP", bounds=bnds, constraints=cons,
options={"ftol": 1e-12, "maxiter": 1000})
return res.x
def sc_gap(unit_idx, pool_idx, Ymat):
'''Build a synthetic control for unit_idx from pool_idx; return gap path and pre/post RMSPE.'''
treated_pre = Ymat[unit_idx, pre]
donors_pre = Ymat[np.ix_(pool_idx, pre)]
w = synth_weights(treated_pre, donors_pre)
synth = w @ Ymat[np.ix_(pool_idx, np.arange(len(years)))] # synthetic path over ALL years
gap = Ymat[unit_idx] - synth
rmspe_pre = np.sqrt(np.mean(gap[pre] ** 2))
rmspe_post = np.sqrt(np.mean(gap[post] ** 2))
return gap, rmspe_pre, rmspe_post, w
donor_ids = np.arange(1, N) # rows 1..30 are the donor pool
gap0, pre0, post0, w0 = sc_gap(0, donor_ids, Y)
order = np.argsort(w0)[::-1]
print("Top donor weights (hand-rolled):")
for k in order[:5]:
print(f" {unit_names[k + 1]:>9s} : {w0[k]:.3f}")
print(f" (active donors with w > 0.01: {(w0 > 0.01).sum()} of {J})")
print(f"\nTreated pre-RMSPE = {pre0:.4f}")
print(f"Treated post-RMSPE = {post0:.4f}")
print(f"post/pre ratio = {post0 / pre0:.2f}")
print(f"Mean post-treatment gap (estimated effect) = {gap0[post].mean():.3f} (true = {TRUE_EFFECT})")The weights are sparse: a handful of donors carry the synthetic control and the rest sit at (essentially) zero. That sparsity is the gift of optimizing over the simplex — it gives Priya a sentence she can defend in a paper (“the treated state is best reproduced by a blend of a few donors with similar factor exposure”) rather than forty coefficients of both signs. The pre-RMSPE is tiny (about 0.035, far below the premium scale of \sim\!10), so the synthetic twin tracks the treated state almost perfectly through the pre-period. The post-RMSPE is large (about 1.5) — the divergence after 2018. Their ratio, about 44, is the chapter’s headline diagnostic, and the estimated effect of \approx 1.5 sits right next to the true 1.4.
3. The same fit with pysyncon — confirming the package agrees
A hand-rolled solve is great for seeing the mechanics, but real studies use a tested package so the \mathbf{V}-matrix tuning, predictor handling, and inference are done correctly. We reproduce the fit with pysyncon (from pysyncon import Dataprep, Synth). Its Dataprep wants a tidy long-format frame and a list of pre-period “special predictors”; here we feed it the premium in each pre-year so it matches on the full pre-period path, exactly like our hand-rolled objective. Synth.fit then solves for both the donor weights \mathbf{w} and the predictor-weight matrix \mathbf{V}.
The two routes use different optimizers (our SLSQP on a fixed \mathbf{V}=\mathbf{I} vs. pysyncon’s Nelder-Mead over \mathbf{V}), so the exact weights can differ slightly — but the fit, the gap, and the estimated effect should agree to within rounding, which is the real test that we built the right thing.
pre_years = [int(y) for y in years[pre]]
post_years = [int(y) for y in years[post]]
dataprep = Dataprep(
foo=panel_long,
predictors=["premium"], # base predictor (overall pre mean; dominated by specials below)
predictors_op="mean",
dependent="premium",
unit_variable="state",
time_variable="year",
treatment_identifier="treated",
controls_identifier=[n for n in unit_names if n != "treated"],
time_predictors_prior=pre_years,
time_optimize_ssr=pre_years,
# one "special predictor" per pre-year => match the entire pre-period premium path
special_predictors=[("premium", [y], "mean") for y in pre_years],
)
synth = Synth()
synth.fit(dataprep=dataprep, optim_method="Nelder-Mead", optim_initial="equal")
w_pys = synth.weights(round=10) # pandas Series indexed by donor name
att = synth.att(time_period=post_years) # package's own ATT and (placebo-style) SE
print("Top donor weights (pysyncon):")
print(w_pys.sort_values(ascending=False).head(5).round(3).to_string())
print(f"\npysyncon ATT over post-period = {att['att']:.3f} (true = {TRUE_EFFECT})")# Reconstruct the pysyncon synthetic path over ALL years from its weights, then compare diagnostics.
wide = panel_long.pivot(index="year", columns="state", values="premium")
synth_path_pys = wide[w_pys.index].to_numpy() @ w_pys.to_numpy()
treated_path = wide["treated"].to_numpy()
gap_pys = treated_path - synth_path_pys
pre_mask_w = wide.index.to_numpy() <= T0_YEAR
pre_pys = np.sqrt(np.mean(gap_pys[pre_mask_w] ** 2))
post_pys = np.sqrt(np.mean(gap_pys[~pre_mask_w] ** 2))
compare = pd.DataFrame({
"pre_RMSPE": [pre0, pre_pys],
"post_RMSPE": [post0, post_pys],
"ratio": [post0 / pre0, post_pys / pre_pys],
"est_effect": [gap0[post].mean(), gap_pys[~pre_mask_w].mean()],
}, index=["hand-rolled (scipy SLSQP)", "pysyncon (Nelder-Mead)"])
print(compare.to_string(float_format=lambda v: f"{v:.4f}"))
print("\nThe two routes agree to rounding -> the hand-rolled solve is the real thing, not a toy.")See the match and the divergence
Two pictures tell the story. Left: the treated state and its synthetic twin overlaid — they lie on top of each other through the pre-period (the credential), then split after 2018. Right: the gap between them, which is flat near zero before treatment and jumps to roughly the true +1.4 after. The flat pre-period gap is what earns the right to read the post-period gap as causal.
synth_path_hand = treated_path - gap0 # hand-rolled synthetic = treated - gap
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
ax = axes[0]
ax.plot(years, Y[0], color="black", lw=2.5, label="treated (actual)")
ax.plot(years, synth_path_hand, color="steelblue", lw=2, ls="--", label="synthetic control")
ax.axvline(T0_YEAR + 0.5, color="crimson", ls=":", lw=1.2)
ax.set_xlabel("year"); ax.set_ylabel("premium ($00s)")
ax.set_title("Treated vs. synthetic: matched pre, diverging post")
ax.legend(loc="upper left")
ax = axes[1]
ax.plot(years, gap0, color="black", lw=2.5, label="gap (treated - synthetic)")
ax.axhline(0, color="0.5", lw=1)
ax.axhline(TRUE_EFFECT, color="seagreen", ls="--", lw=1.2, label=f"true effect = {TRUE_EFFECT}")
ax.axvline(T0_YEAR + 0.5, color="crimson", ls=":", lw=1.2, label="mandate (2018)")
ax.set_xlabel("year"); ax.set_ylabel("gap ($00s)")
ax.set_title("The estimated effect is the post-period gap")
ax.legend(loc="upper left")
plt.show()4. Inference when your sample is one: the placebo permutation test
We have a point estimate of \approx +1.5. Is it real, or just the kind of gap any unit shows against its synthetic control from noise and imperfect fit? With one treated unit there is no cross-sectional variation to power a standard error; a t-statistic here would be false precision. The honest alternative (Abadie, Diamond and Hainmueller, 2010) is a placebo permutation test: pretend each donor is the treated unit, build its synthetic control from the other donors, and collect its gap path. That gives a reference distribution of “effects” for units we know were never treated.
To compare fairly, scale each unit’s post-period gap by its own pre-period fit — the post/pre RMSPE ratio
r_i = \frac{\text{RMSPE}_{\text{post},\,i}}{\text{RMSPE}_{\text{pre},\,i}},
so that a unit only looks impressive if it fit well before (small denominator) and diverged sharply after (large numerator). The permutation p-value is the fraction of all J+1 units whose ratio is at least as large as the treated unit’s, p = \#\{i : r_i \ge r_1\} / (J+1).
def placebo_run(Ymat):
'''Run SC on the real treated unit and on every donor as fake-treated.
Returns gap paths (dict idx->gap), RMSPE ratios (dict idx->ratio).'''
gaps, ratios = {}, {}
g0, pr0, po0, _ = sc_gap(0, donor_ids, Ymat) # real treated unit
gaps[0] = g0; ratios[0] = po0 / max(pr0, 1e-8)
for j in donor_ids: # each donor as fake-treated
pool = donor_ids[donor_ids != j] # the OTHER donors (treated already excluded)
gj, prj, poj, _ = sc_gap(j, pool, Ymat)
gaps[j] = gj; ratios[j] = poj / max(prj, 1e-8)
return gaps, ratios
gaps, ratios = placebo_run(Y)
r1 = ratios[0]
all_r = np.array([ratios[i] for i in range(N)])
rank = int((all_r >= r1).sum())
p_val = rank / N
print(f"Treated RMSPE ratio = {r1:.2f}")
print(f"Rank among {N} units = {rank} (1 = most extreme)")
print(f"Permutation p-value = {rank}/{N} = {p_val:.4f}")
print(f"\nChapter target: rank 1 of 31, p = 1/31 = {1/31:.4f}.")The gray-spaghetti plot: see the outlier
Overlay every placebo gap path in gray and the treated unit’s gap in black. If the treated line sits inside the gray cloud before 2018 and wanders out beyond it after, the effect is real and you can see it. If the treated line were just one more strand in the spaghetti, there would be no effect no matter what point estimate we computed. (We follow the chapter’s transparent trimming rule and drop placebo donors whose pre-RMSPE is more than 5x the treated unit’s, since those were never matched well enough to produce an interpretable gap.)
# pre-RMSPE per unit, to apply the 5x trimming rule for the plot only.
pre_rmspe = {i: np.sqrt(np.mean(gaps[i][pre] ** 2)) for i in range(N)}
trim = 5.0 * pre_rmspe[0]
fig, ax = plt.subplots()
n_drawn = 0
for j in donor_ids:
if pre_rmspe[j] <= trim:
ax.plot(years, gaps[j], color="0.75", lw=1, zorder=1)
n_drawn += 1
ax.plot([], [], color="0.75", lw=1, label=f"{n_drawn} placebo donors")
ax.plot(years, gaps[0], color="black", lw=2.5, zorder=3, label="treated state")
ax.axhline(0, color="0.4", lw=0.8)
ax.axvline(T0_YEAR + 0.5, color="crimson", ls=":", lw=1.2, label="mandate (2018)")
ax.set_xlabel("year"); ax.set_ylabel("gap (actual - synthetic, $00s)")
ax.set_title("Placebo gaps (gray) vs. treated gap (black): the treated unit breaks out")
ax.legend(loc="upper left")
plt.show()
# Companion: the distribution of RMSPE ratios, treated marked.
fig, ax = plt.subplots(figsize=(9, 4))
ax.hist(all_r[1:], bins=15, color="0.8", edgecolor="0.5", label="placebo donors")
ax.axvline(r1, color="black", lw=2.5, label=f"treated (ratio={r1:.1f}, rank {rank}/{N})")
ax.set_xlabel("post/pre RMSPE ratio"); ax.set_ylabel("count")
ax.set_title("The treated unit's RMSPE ratio is the largest in the reference set")
ax.legend()
plt.show()5. No-effect sanity check: does the method correctly find nothing?
A test that always declares “significant” is worthless. So we re-run the entire pipeline on a panel with the injected jump removed — same factors, same loadings, same noise draws, but no treatment effect. If the method is honest, the treated unit’s RMSPE ratio should collapse back into the placebo crowd and the permutation p-value should become large. The chapter predicts rank \approx 5 of 31 and p \approx 0.16.
Y_null = Y.copy()
Y_null[0, post] -= TRUE_EFFECT # undo the jump on the treated state's post-period
gaps_n, ratios_n = placebo_run(Y_null)
r1_n = ratios_n[0]
all_r_n = np.array([ratios_n[i] for i in range(N)])
rank_n = int((all_r_n >= r1_n).sum())
p_n = rank_n / N
eff_n = gaps_n[0][post].mean()
print("NO-EFFECT panel:")
print(f" estimated effect (mean post gap) = {eff_n:+.3f} (true = 0)")
print(f" treated RMSPE ratio = {r1_n:.2f}")
print(f" rank = {rank_n} of {N}; permutation p = {p_n:.3f}")
print(f"\nWITH-EFFECT panel (for contrast): effect={gap0[post].mean():+.3f}, "
f"ratio={r1:.2f}, rank={rank}/{N}, p={p_val:.3f}")
print("\nWith no real effect the treated unit is just another strand: p is large, as it should be.")fig, ax = plt.subplots()
pre_rmspe_n = {i: np.sqrt(np.mean(gaps_n[i][pre] ** 2)) for i in range(N)}
trim_n = 5.0 * pre_rmspe_n[0]
for j in donor_ids:
if pre_rmspe_n[j] <= trim_n:
ax.plot(years, gaps_n[j], color="0.75", lw=1, zorder=1)
ax.plot(years, gaps_n[0], color="black", lw=2.5, zorder=3, label="treated (no effect injected)")
ax.axhline(0, color="0.4", lw=0.8)
ax.axvline(T0_YEAR + 0.5, color="crimson", ls=":", lw=1.2, label="mandate (2018)")
ax.set_xlabel("year"); ax.set_ylabel("gap ($00s)")
ax.set_title(f"No-effect placebo: treated line stays in the crowd (rank {rank_n}/{N}, p={p_n:.2f})")
ax.legend(loc="upper left")
plt.show()6. Synthetic DiD: unit weights and time weights
Synthetic DiD (Arkhangelsky, Athey, Hirshberg, Imbens and Wager, 2021) marries the two methods. Like synthetic control it reweights units — but only to match the treated unit’s pre-period trend (it allows a constant level gap, which it then differences away DiD-style), so it needs a far weaker target than an exact level match. Like nothing in plain DiD, it also reweights time — choosing weights \hat\lambda_t on the pre-years so the pre-periods that best predict the post-period get more say. Then it runs a weighted two-way fixed-effects DiD.
We implement a faithful version directly:
- Unit weights \hat w_j \ge 0, \sum_j \hat w_j = 1: minimize \sum_{t \le T_0}(Y_{1t} - w_0 - \sum_j \hat w_j Y_{jt})^2 + \zeta^2 \lVert \mathbf{w}\rVert^2 over w_j and a free intercept w_0. The intercept is the allowed level gap; the small ridge \zeta^2 (Arkhangelsky et al.’s regularizer) keeps the weights from spreading on noise.
- Time weights \hat\lambda_t \ge 0, \sum_t \hat\lambda_t = 1: choose pre-year weights so the weighted pre-period predicts each donor’s post-period mean (again with a free intercept).
- The estimate is the weighted double difference \hat\tau = \big(\bar Y_{1,\text{post}} - \sum_t \hat\lambda_t Y_{1t}\big) - \sum_j \hat w_j\big(\bar Y_{j,\text{post}} - \sum_t \hat\lambda_t Y_{jt}\big).
T0 = int(pre.sum()); Tpost = int(post.sum())
# --- unit weights: match treated PRE trend with a level offset + Arkhangelsky ridge ---
Y0_pre = Y[1:, pre] # (J, T0) donor pre-period
Y1_pre = Y[0, pre] # (T0,) treated pre-period
# regularizer scale: (n_post)^(1/2) times the std of first-differenced donor outcomes (their recipe)
zeta = (Tpost ** 0.5) * np.std(np.diff(Y[1:, :], axis=1))
def unit_loss(p):
b, w = p[0], p[1:]
return np.sum((Y1_pre - (b + w @ Y0_pre)) ** 2) + (zeta ** 2) * np.sum(w ** 2)
cons_u = ({"type": "eq", "fun": lambda p: p[1:].sum() - 1.0},)
bnds_u = [(None, None)] + [(0.0, 1.0)] * J
res_u = minimize(unit_loss, np.concatenate([[0.0], np.full(J, 1.0 / J)]),
method="SLSQP", bounds=bnds_u, constraints=cons_u,
options={"ftol": 1e-12, "maxiter": 2000})
w_unit = res_u.x[1:]
# --- time weights: weighted pre-period predicts each donor's post-period mean ---
Y0_post_mean = Y[1:, post].mean(axis=1) # (J,)
def time_loss(p):
b, lam = p[0], p[1:]
return np.sum((Y0_post_mean - (b + Y0_pre @ lam)) ** 2)
cons_t = ({"type": "eq", "fun": lambda p: p[1:].sum() - 1.0},)
bnds_t = [(None, None)] + [(0.0, 1.0)] * T0
res_t = minimize(time_loss, np.concatenate([[0.0], np.full(T0, 1.0 / T0)]),
method="SLSQP", bounds=bnds_t, constraints=cons_t,
options={"ftol": 1e-12, "maxiter": 2000})
lam_time = res_t.x[1:]
# --- weighted double difference ---
treated_post = Y[0, post].mean()
treated_pre_w = lam_time @ Y[0, pre]
donor_post = w_unit @ Y[1:, post].mean(axis=1)
donor_pre_w = w_unit @ (Y[1:, pre] @ lam_time)
tau_sdid = (treated_post - treated_pre_w) - (donor_post - donor_pre_w)
print(f"SDID unit weights: {(w_unit > 0.01).sum()} active of {J}, sum={w_unit.sum():.3f}")
print(f"SDID time weights: sum={lam_time.sum():.3f}, most weight on years "
f"{years[pre][np.argsort(lam_time)[::-1][:3]].tolist()}")
print(f"\nSynthetic DiD tau = {tau_sdid:.3f} (true = {TRUE_EFFECT})")A placebo variance for SDID, and the three-way comparison
SDID comes with a principled inference procedure. With one treated unit we use the placebo variance: reassign the “treatment” to each donor in turn, recompute \hat\tau for that fake-treated unit (which has no real effect), and use the spread of those placebo estimates as the standard error. That gives a confidence interval without leaning on large-N asymptotics. Finally we line SDID up against plain synthetic control and plain DiD on the same panel.
def sdid_tau(Ymat, treated_idx, pool_idx):
'''Compute the SDID estimate treating treated_idx as treated, pool_idx as donors.'''
Y0p = Ymat[np.ix_(pool_idx, pre)]
Y1p = Ymat[treated_idx, pre]
Jn = len(pool_idx)
z = (Tpost ** 0.5) * np.std(np.diff(Ymat[pool_idx, :], axis=1))
ru = minimize(lambda p: np.sum((Y1p - (p[0] + p[1:] @ Y0p)) ** 2) + (z ** 2) * np.sum(p[1:] ** 2),
np.concatenate([[0.0], np.full(Jn, 1.0 / Jn)]), method="SLSQP",
bounds=[(None, None)] + [(0.0, 1.0)] * Jn,
constraints=({"type": "eq", "fun": lambda p: p[1:].sum() - 1.0},),
options={"ftol": 1e-12, "maxiter": 2000}).x[1:]
ypm = Ymat[np.ix_(pool_idx, post)].mean(axis=1)
rt = minimize(lambda p: np.sum((ypm - (p[0] + Y0p @ p[1:])) ** 2),
np.concatenate([[0.0], np.full(T0, 1.0 / T0)]), method="SLSQP",
bounds=[(None, None)] + [(0.0, 1.0)] * T0,
constraints=({"type": "eq", "fun": lambda p: p[1:].sum() - 1.0},),
options={"ftol": 1e-12, "maxiter": 2000}).x[1:]
tp = Ymat[treated_idx, post].mean(); tpw = rt @ Ymat[treated_idx, pre]
dp = ru @ Ymat[np.ix_(pool_idx, post)].mean(axis=1); dpw = ru @ (Ymat[np.ix_(pool_idx, pre)] @ rt)
return (tp - tpw) - (dp - dpw)
# placebo SE: treat each donor as fake-treated, gather its SDID tau (no real effect there)
placebo_taus = [sdid_tau(Y, j, donor_ids[donor_ids != j]) for j in donor_ids]
se_sdid = np.std(placebo_taus, ddof=1)
ci_lo, ci_hi = tau_sdid - 1.96 * se_sdid, tau_sdid + 1.96 * se_sdid
# plain DiD (2x2, equal-weight control = synthetic control with frozen weights)
did = ((Y[0, post].mean() - Y[0, pre].mean())
- (Y[1:, post].mean(axis=1).mean() - Y[1:, pre].mean(axis=1).mean()))
summary = pd.DataFrame({
"estimate": [gap0[post].mean(), tau_sdid, did],
"inference": [f"placebo p={p_val:.3f} (rank {rank}/{N})",
f"95% CI [{ci_lo:.2f}, {ci_hi:.2f}], SE={se_sdid:.3f}",
"parallel-trends (not credible w/ 1 unit)"],
}, index=["Synthetic Control", "Synthetic DiD", "Plain DiD (equal wts)"])
print(f"True effect = {TRUE_EFFECT}\n")
print(summary.to_string())
print("\nAll three bracket the truth; SDID is closest and ships a usable confidence interval.")What we found, and Your Turn
On a panel where we know the answer is +1.4, all three designs land in the right neighborhood, but they earn their estimates very differently:
| method | estimate | how it argues | when to reach for it |
|---|---|---|---|
| Synthetic control | \approx 1.5 | convex donor blend, exact pre-period match, placebo p \approx 0.032 | one unit, excellent pre-fit, you want a weight vector you can name |
| Synthetic DiD | \approx 1.34 | unit + time weights, differences out a level gap, placebo SE | poor level match but good trend match; a few treated units |
| Plain DiD | \approx 1.43 | equal-weight control + parallel trends | many treated and controls with visibly parallel pre-trends |
The deep lessons: build the counterfactual, don’t assume it (the convex match is parallel trends engineered and displayed); convexity is the method’s conscience (it cannot extrapolate or fake a fit, and it confesses an edge unit it cannot match); inference is permutation, not asymptotics (with one unit the placebo rank is your p-value, floored at 1/(J+1)); and synthetic DiD generalizes the toolkit by asking only for a matched trend.
Your Turn.
Edge unit. Set
L1[0], L2[0] = 1.35, 1.05(push the treated state to the high corner of the donor cloud) and re-run Sections 2-4. Watch the pre-RMSPE balloon and the synthetic control underfit from below — the convexity constraint refusing to extrapolate. Does the placebo p-value still flag the effect? Which method does the chapter say is better for an edge unit, and does your run agree?Donor pool size. Shrink the pool to
J = 18(regenerate the panel). The smallest possible permutation p-value is now 1/19 \approx 0.053 — you literally cannot get below 5%. Re-run and confirm the floor. What does this say about how many donors a credible study needs?Anticipation. Add a small pre-leak —
Y[0, years == 2017] += 0.5— to mimic insurers repricing the year before the mandate. Which assumption does this violate? Re-fit and show the pre-period gap is no longer flat at 2017, and that the estimated effect is biased toward zero. How would you redefine T_0 to fix it?