One treated province, many possible controls
When one province adopts a carbon cap-and-trade programme, DiD needs a control province with parallel trends, and usually no single province qualifies. The synthetic control method (Abadie, Diamond and Hainmueller) builds one: a weighted average of untreated units (the "donor pool") chosen so that it tracks the treated unit closely before the policy. After the policy, the synthetic unit shows what would have happened without it.
Two constraints make the method transparent:
- weights are non-negative, and
- weights sum to one,
so the synthetic unit is an interpolation of real units, never an extrapolation. You can name exactly which provinces make up "synthetic Province A" and in what proportion.
!pip install -q numpy pandas scipy matplotlib seaborn
import numpy as np, pandas as pd
from scipy.optimize import minimize
import matplotlib.pyplot as plt, seaborn as sns
rng = np.random.default_rng(42)
sns.set_theme(style="whitegrid")
A cap-and-trade programme in Province A
Province A introduces cap-and-trade in 2018. Ten other provinces/states form the donor pool. Each has its own emissions trajectory driven by a few common factors (industrial cycles, energy prices). The true effect is −15 units per year from 2018.
years = np.arange(2005, 2025)
T0 = 2018
pre = years < T0
J = 10 # donors
factors = np.column_stack([np.ones(len(years)), (years - 2005) / 10, np.sin((years - 2005) / 3)])
loadings = rng.uniform([80, 0, -5], [140, 30, 10], size=(J, 3))
donors = loadings @ factors.T + rng.normal(0, 1.5, (J, len(years))) # J x T
true_w = np.array([0.4, 0.5, 0.1] + [0] * (J - 3)) # A resembles donors 0, 1, 2
effect = np.where(years >= T0, -15.0, 0.0)
treated = true_w @ donors + effect + rng.normal(0, 1.0, len(years))
panel = pd.DataFrame(donors.T, index=years, columns=[f"Donor {j}" for j in range(J)])
panel["Province A"] = treated
panel.iloc[:, [0, 1, 2, -1]].plot(figsize=(10, 4), title="Emissions: Province A and three donors"); plt.show()
Solving for the weights
Find w ≥ 0 with Σw = 1 that minimises the pre-treatment squared gap between Province A and the weighted donors. This is a small quadratic programme; scipy's SLSQP solver handles it.
def synth_weights(y_treated_pre, Y_donors_pre):
J = Y_donors_pre.shape[0]
loss = lambda w: np.sum((y_treated_pre - w @ Y_donors_pre) ** 2)
res = minimize(loss, np.full(J, 1 / J), method="SLSQP",
bounds=[(0, 1)] * J, constraints={"type": "eq", "fun": lambda w: w.sum() - 1})
return res.x
w = synth_weights(treated[pre], donors[:, pre])
print(pd.Series(w, index=panel.columns[:J]).round(3).sort_values(ascending=False).head(5))
Most donors get (almost) zero weight: synthetic controls are typically sparse, which makes them easy to explain.
Fit, gap and interpretation
synthetic = w @ donors
gap = treated - synthetic
rmspe_pre = np.sqrt(np.mean(gap[pre] ** 2))
print(f"Pre-treatment RMSPE: {rmspe_pre:.2f}")
print(f"Average post-treatment gap: {gap[~pre].mean():.2f} (truth -15)")
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
axes[0].plot(years, treated, color="red", lw=2, label="Province A (actual)")
axes[0].plot(years, synthetic, color="blue", ls="--", lw=2, label="Synthetic Province A")
axes[0].axvline(T0, color="black", ls=":"); axes[0].legend(); axes[0].set_title("Actual vs synthetic")
axes[1].bar(years, gap, color=np.where(pre, "#94a3b8", "#ef4444"))
axes[1].axvline(T0 - 0.5, color="black", ls=":"); axes[1].set_title("Gap = estimated effect each year")
plt.tight_layout(); plt.show()
A good pre-period fit over many years is what makes the counterfactual credible: if a weighted mix of donors tracked Province A through several cycles, it is plausible it would have kept tracking it. A poor fit (large pre-RMSPE) means no convex combination of donors resembles the treated unit, and the method should not be used, rather than "fixed" by allowing negative weights.
Exercise
Compare the synthetic control with two alternatives: (1) the simple average of all donors as the control (a DiD-style comparison), and (2) unconstrained OLS weights (allowing negative weights and no sum-to-one). Compare pre-period fit and the estimated effect. Why are OLS weights dangerous with few pre-periods?
Show a solution
avg = donors.mean(axis=0)
did = (treated[~pre] - avg[~pre]).mean() - (treated[pre] - avg[pre]).mean()
ols_w, *_ = np.linalg.lstsq(donors[:, pre].T, treated[pre], rcond=None)
ols_gap = treated - ols_w @ donors
print(f"DiD vs donor average: {did:.2f}")
print(f"OLS weights: effect {ols_gap[~pre].mean():.2f}, pre-RMSPE {np.sqrt(np.mean(ols_gap[pre]**2)):.2f}, "
f"weights range [{ols_w.min():.2f}, {ols_w.max():.2f}]")
# OLS fits the pre-period almost perfectly by overfitting (13 periods, 10 weights) and extrapolates with
# negative weights; its post-period counterfactual is unreliable.
Checklist and recap
- I can explain the synthetic control idea and its two weight constraints.
- I can solve for weights with scipy and report which donors matter.
- I judge credibility by pre-treatment fit before reading the post-treatment gap.