{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2a2e4d28",
   "metadata": {},
   "source": [
    "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/rpasquini/econometrics-book-notebooks/blob/main/en/ch4_inferencia_causal/6_synthetic_control_notebook.ipynb)\n",
    "\n",
    "# Code companion: The Synthetic Control Method\n",
    "\n",
    "This notebook accompanies the section [The Synthetic Control Method](6_synthetic_control.md). That section's interactive dashboards let you *feel* how the weights are built and how permutation inference works by moving controls; this notebook shows how the same problem is solved in a **real analysis pipeline**, combining Dashboard 1's mechanics (weight construction) and Dashboard 3's (permutation inference) into one applied workflow over two datasets:\n",
    "\n",
    "1. **The same Uber/Austin simulation** that generates Dashboard 1 — frozen with a fixed seed into a `.csv` so the analysis is reproducible.\n",
    "2. **The classic case study**: California's Proposition 99 (Abadie, Diamond & Hainmueller, 2010) — the paper that introduced the synthetic control method.\n",
    "\n",
    "For each dataset we first implement a **from-scratch** version (to see the exact mechanics behind the estimator) and then the same idea using `scpi_pkg`, a package specialized in synthetic control and its prediction intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4db12015",
   "metadata": {},
   "source": [
    "## Citations\n",
    "\n",
    "- Abadie, A., Diamond, A., & Hainmueller, J. (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.\n",
    "- The Proposition 99 data comes from the R package `synthdid` (Arkhangelsky, Athey, Hirshberg, Imbens & Wager), BSD-3-Clause licensed — permits redistribution and reuse with attribution. Repository: [synth-inference/synthdid](https://github.com/synth-inference/synthdid).\n",
    "- Cattaneo, M. D., Feng, Y., Palomba, F., & Titiunik, R. (2025). *scpi: Uncertainty Quantification for Synthetic Control Estimators.* Journal of Statistical Software, 113(1), 1–38.\n",
    "- Cattaneo, M. D., Feng, Y., & Titiunik, R. (2021). *Prediction Intervals for Synthetic Control Methods.* Journal of the American Statistical Association, 116(536), 1865–1880.\n",
    "- (For future reference) Abadie, A., Diamond, A., & Hainmueller, J. (2015). *Comparative Politics and the Synthetic Control Method.* American Journal of Political Science, 59(2), 495–510 — the West Germany reunification study, whose data is already downloaded under `data/` but not yet used in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc340fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# scpi_pkg is not preinstalled on Colab -- installs quickly, safe to skip if already present\n",
    "try:\n",
    "    import scpi_pkg  # noqa: F401\n",
    "except ImportError:\n",
    "    %pip install -q scpi_pkg\n",
    "\n",
    "from scpi_pkg.scdata import scdata\n",
    "from scpi_pkg.scest import scest\n",
    "from scpi_pkg.scpi import scpi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164600ed",
   "metadata": {},
   "source": [
    "## Part 1 — The Uber/Austin simulation\n",
    "\n",
    "`rides_thousands` is weekly completed rides (thousands), in Austin and an 18-city donor pool. The first 20 weeks are pre-intervention; the driver-incentive program starts in week 21. This is the same data generated by Dashboard 1, frozen with the same random seed (`comparison_mode=\"optimized\"`, `true_effect=10%`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a31e584",
   "metadata": {},
   "outputs": [],
   "source": [
    "uber = pd.read_csv('https://raw.githubusercontent.com/rpasquini/econometrics-book-notebooks/main/en/ch4_inferencia_causal/data/synthetic_control_uber_austin.csv')\n",
    "uber.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c8f25a7",
   "metadata": {},
   "source": [
    "### Estimator 1a — Synthetic control weights, from scratch\n",
    "\n",
    "We reproduce the Dashboard 1 mechanics exactly: 8 predictors per city (4 demographic variables + 4 lagged ride observations at weeks 5, 10, 15, and 20), and solve\n",
    "\n",
    "$$\\min_W (X_1 - X_0 W)'V(X_1-X_0W) \\quad \\text{s.t. } W\\ge 0,\\ \\mathbf{1}'W=1$$\n",
    "\n",
    "with $V$ the \"simple selector\" (inverse-variance weighting of each predictor)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "524a8d25",
   "metadata": {},
   "outputs": [],
   "source": [
    "T0 = 20\n",
    "HISTORY_WEEKS = [5, 10, 15, 20]\n",
    "DEMO_COLS = [\"population_m\", \"income_k\", \"drivers_k\", \"lyft_share_pct\"]\n",
    "\n",
    "def build_predictor_matrix(df, cities):\n",
    "    rows = []\n",
    "    for city in cities:\n",
    "        sub = df[df.city == city]\n",
    "        demo = sub[DEMO_COLS].iloc[0].values\n",
    "        hist = sub[sub.week.isin(HISTORY_WEEKS)].sort_values(\"week\")[\"rides_thousands\"].values\n",
    "        rows.append(np.concatenate([demo, hist]))\n",
    "    return np.array(rows).T  # (8, n_cities)\n",
    "\n",
    "def solve_weights(x1, x0):\n",
    "    \"\"\"Constrained least squares: minimize ||x1 - x0 @ w||_V s.t. w >= 0, sum(w) = 1.\"\"\"\n",
    "    J = x0.shape[1]\n",
    "    V = 1.0 / np.maximum(x0.var(axis=1), 1e-6)\n",
    "\n",
    "    def loss(w):\n",
    "        diff = x1 - x0 @ w\n",
    "        return float(diff @ (V * diff))\n",
    "\n",
    "    res = minimize(loss, x0=np.full(J, 1.0 / J), method=\"SLSQP\",\n",
    "                    bounds=[(0.0, 1.0)] * J,\n",
    "                    constraints=[{\"type\": \"eq\", \"fun\": lambda w: w.sum() - 1.0}],\n",
    "                    options={\"maxiter\": 300, \"ftol\": 1e-10})\n",
    "    w = np.clip(res.x, 0, None)\n",
    "    total = w.sum()\n",
    "    return w / total if total > 0 else np.full(J, 1.0 / J)\n",
    "\n",
    "donor_cities = sorted(uber.loc[uber.is_treated_unit == 0, \"city\"].unique())\n",
    "X = build_predictor_matrix(uber, [\"Austin\"] + donor_cities)\n",
    "w_uber = solve_weights(X[:, 0], X[:, 1:])\n",
    "\n",
    "weights_table = pd.Series(w_uber, index=donor_cities).sort_values(ascending=False)\n",
    "weights_table[weights_table > 0.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5af70c5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pivot_uber = uber.pivot(index=\"city\", columns=\"week\", values=\"rides_thousands\")\n",
    "y_austin = pivot_uber.loc[\"Austin\"].sort_index().values\n",
    "y_donors = pivot_uber.loc[donor_cities].sort_index(axis=1).values\n",
    "synthetic_uber = y_donors.T @ w_uber\n",
    "gap_uber = y_austin - synthetic_uber\n",
    "\n",
    "rmspe_pre_uber = np.sqrt(np.mean(gap_uber[:T0] ** 2))\n",
    "tau_hat_uber = gap_uber[T0:].mean()\n",
    "print(f\"Pre-treatment RMSPE: {rmspe_pre_uber:.2f}\")\n",
    "print(f\"tau_hat (mean post-treatment gap): {tau_hat_uber:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c13ca453",
   "metadata": {},
   "outputs": [],
   "source": [
    "weeks = np.arange(1, len(y_austin) + 1)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))\n",
    "\n",
    "axes[0].plot(weeks, y_austin, label=\"Austin\", color=\"#8e44ad\", linewidth=2)\n",
    "axes[0].plot(weeks, synthetic_uber, label=\"Synthetic control\", color=\"#e74c3c\", linestyle=\"--\", linewidth=2)\n",
    "axes[0].axvline(T0 + 0.5, color=\"gray\", linestyle=\":\")\n",
    "axes[0].set_xlabel(\"Week\"); axes[0].set_ylabel(\"Rides (thousands)\")\n",
    "axes[0].set_title(\"Austin vs. synthetic control\"); axes[0].legend()\n",
    "\n",
    "axes[1].plot(weeks, gap_uber, color=\"#2c3e50\")\n",
    "axes[1].axhline(0, color=\"lightgray\")\n",
    "axes[1].axvline(T0 + 0.5, color=\"gray\", linestyle=\":\")\n",
    "axes[1].set_xlabel(\"Week\"); axes[1].set_ylabel(\"Austin - Synthetic\")\n",
    "axes[1].set_title(\"Estimated gap\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa0bc7ac",
   "metadata": {},
   "source": [
    "**Interpretation.** The pre-treatment fit is nearly exact and the weights concentrate on a handful of cities (Kansas City, Boise, Columbus...) — the same sparsity seen in Dashboard 1 under `comparison_mode=\"optimized\"`. $\\hat\\tau \\approx 9.8$ is close to the true simulated effect (10%).\n",
    "\n",
    "### Estimator 1b — The same problem with `scpi_pkg`\n",
    "\n",
    "`scpi_pkg` is designed around outcome *features* and adjustment covariates (`cov_adj`), not a hand-built ADH-style predictor matrix. With its default configuration (no `features` specified), it uses **the entire pre-treatment outcome trajectory** as the predictor — here, all 20 weeks — rather than the specific set of 8 predictors (demographics + 4 lagged weeks) that Dashboard 1 builds.\n",
    "\n",
    "**Important note:** because of this, don't expect `scpi_pkg`'s weight vector to exactly reproduce the one from the previous cell — they solve different specifications of the matching problem, even though both are reasonable ways to implement synthetic control. This is a good example of \"the same idea, implemented with different tools\" not always giving an identical answer. `scpi_pkg`'s default $V$ is also `'separate'` (a partitioned least-squares variant), different from the \"simple selector\" (inverse variance) Dashboard 1 uses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e738130",
   "metadata": {},
   "outputs": [],
   "source": [
    "period_pre_uber = np.arange(1, T0 + 1)\n",
    "period_post_uber = np.arange(T0 + 1, uber.week.max() + 1)\n",
    "\n",
    "data_uber = scdata(df=uber, id_var=\"city\", time_var=\"week\", outcome_var=\"rides_thousands\",\n",
    "                    period_pre=period_pre_uber, period_post=period_post_uber,\n",
    "                    unit_tr=\"Austin\", unit_co=donor_cities, constant=False, verbose=False)\n",
    "est_uber = scest(data_uber, w_constr={\"name\": \"simplex\"})\n",
    "\n",
    "w_scpi_uber = est_uber.w.iloc[:, 0].sort_values(ascending=False)\n",
    "w_scpi_uber[w_scpi_uber > 0.01]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32fa8a90",
   "metadata": {},
   "source": [
    "## Part 2 — A classic case study: Proposition 99 in California\n",
    "\n",
    "We now apply the same pipeline to the real data from Abadie, Diamond & Hainmueller (2010): annual cigarette consumption (`PacksPerCapita`, packs per capita) for California and 38 donor states, 1970–2000. California passed Proposition 99 (a large tobacco tax and control campaign) in 1988, effective from 1989."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ac82227",
   "metadata": {},
   "outputs": [],
   "source": [
    "prop99 = pd.read_csv('https://raw.githubusercontent.com/rpasquini/econometrics-book-notebooks/main/en/ch4_inferencia_causal/data/california_prop99.csv', sep=';')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ec6f1c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "PRE_YEARS = list(range(1970, 1989))\n",
    "POST_YEARS = list(range(1989, 2001))\n",
    "ALL_YEARS = PRE_YEARS + POST_YEARS\n",
    "T0_PROP99 = len(PRE_YEARS)\n",
    "\n",
    "pivot_prop99 = prop99.pivot(index=\"State\", columns=\"Year\", values=\"PacksPerCapita\")\n",
    "donor_states = [s for s in pivot_prop99.index if s != \"California\"]\n",
    "prop99.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "839bf77f",
   "metadata": {},
   "source": [
    "### Estimator 1a (repeated) — From-scratch weights, on real data\n",
    "\n",
    "Here the only predictor is the outcome itself (packs per capita) in each of the 19 pre-treatment years — the classic ADH approach when no additional covariates are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e67dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_prop99 = pivot_prop99.loc[\"California\", PRE_YEARS].values\n",
    "x0_prop99 = pivot_prop99.loc[donor_states, PRE_YEARS].values.T\n",
    "w_prop99 = solve_weights(x1_prop99, x0_prop99)\n",
    "\n",
    "weights_prop99 = pd.Series(w_prop99, index=donor_states).sort_values(ascending=False)\n",
    "weights_prop99[weights_prop99 > 0.01]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fec471b",
   "metadata": {},
   "outputs": [],
   "source": [
    "synthetic_prop99 = pivot_prop99.loc[donor_states, ALL_YEARS].values.T @ w_prop99\n",
    "actual_prop99 = pivot_prop99.loc[\"California\", ALL_YEARS].values\n",
    "gap_prop99 = actual_prop99 - synthetic_prop99\n",
    "\n",
    "rmspe_pre_prop99 = np.sqrt(np.mean(gap_prop99[:T0_PROP99] ** 2))\n",
    "print(f\"Pre-treatment RMSPE: {rmspe_pre_prop99:.2f}\")\n",
    "print(f\"Gap by 2000: {gap_prop99[-1]:.1f} packs per capita\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1778631",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(ALL_YEARS, actual_prop99, label=\"California\", color=\"#8e44ad\", linewidth=2)\n",
    "ax.plot(ALL_YEARS, synthetic_prop99, label=\"Synthetic control\", color=\"#e74c3c\", linestyle=\"--\", linewidth=2)\n",
    "ax.axvline(1988.5, color=\"gray\", linestyle=\":\")\n",
    "ax.set_xlabel(\"Year\"); ax.set_ylabel(\"Cigarette packs per capita\")\n",
    "ax.set_title(\"Proposition 99: California vs. synthetic control (ADH 2010)\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "711cb62a",
   "metadata": {},
   "source": [
    "**Reproducing the classic result.** The weights concentrate on Utah, Montana, Nevada, Connecticut, and New Hampshire — the same donor set (with very similar weights) reported in Table 2 of ADH (2010). The gap by 2000 (~-27 packs per capita) also matches the published result: Proposition 99 substantially reduced cigarette consumption in California relative to its synthetic control.\n",
    "\n",
    "### Estimator 1b (repeated) — `scpi_pkg` on Prop 99\n",
    "\n",
    "With a single *feature* (the outcome) and no additional covariates, `scpi_pkg`'s default specification coincides with the classic ADH approach — so here we do expect weights very close to the previous cell's (unlike the Uber/Austin case, where demographics introduced a different specification)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c2150f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "period_pre_p99 = np.arange(1970, 1989)\n",
    "period_post_p99 = np.arange(1989, 2001)\n",
    "\n",
    "data_prop99 = scdata(df=prop99, id_var=\"State\", time_var=\"Year\", outcome_var=\"PacksPerCapita\",\n",
    "                      period_pre=period_pre_p99, period_post=period_post_p99,\n",
    "                      unit_tr=\"California\", unit_co=donor_states, constant=False, verbose=False)\n",
    "est_prop99 = scest(data_prop99, w_constr={\"name\": \"simplex\"})\n",
    "\n",
    "w_scpi_prop99 = est_prop99.w.iloc[:, 0].sort_values(ascending=False)\n",
    "w_scpi_prop99[w_scpi_prop99 > 0.01]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08f3caeb",
   "metadata": {},
   "source": [
    "## Estimator 2 — Permutation inference\n",
    "\n",
    "We reproduce Dashboard 3's logic: reassign \"treatment\" to each donor unit in turn (using the rest as its own donor pool), compute $r_j = \\text{RMSPE}_{post,j}/\\text{RMSPE}_{pre,j}$ for each, and get a p-value based on Austin/California's rank within that distribution.\n",
    "\n",
    "### On the Uber/Austin simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3136489a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def leave_one_out_placebo(pivot, all_units, predictor_fn, T0):\n",
    "    \"\"\"Run the leave-one-out placebo procedure for every unit; return r_j and gaps.\"\"\"\n",
    "    gaps, rmspe_pre, rmspe_post = {}, {}, {}\n",
    "    for unit in all_units:\n",
    "        donors = [u for u in all_units if u != unit]\n",
    "        x1, x0 = predictor_fn(unit, donors)\n",
    "        w = solve_weights(x1, x0)\n",
    "        synthetic = pivot.loc[donors].values.T @ w\n",
    "        actual = pivot.loc[unit].values\n",
    "        gap = actual - synthetic\n",
    "        gaps[unit] = gap\n",
    "        rmspe_pre[unit] = np.sqrt(np.mean(gap[:T0] ** 2))\n",
    "        rmspe_post[unit] = np.sqrt(np.mean(gap[T0:] ** 2))\n",
    "    r = {u: rmspe_post[u] / max(rmspe_pre[u], 1e-9) for u in all_units}\n",
    "    return gaps, rmspe_pre, r\n",
    "\n",
    "all_cities = [\"Austin\"] + donor_cities\n",
    "\n",
    "def uber_predictors(unit, donors):\n",
    "    X = build_predictor_matrix(uber, [unit] + donors)\n",
    "    return X[:, 0], X[:, 1:]\n",
    "\n",
    "gaps_uber, rmspe_pre_uber_all, r_uber = leave_one_out_placebo(\n",
    "    pivot_uber, all_cities, uber_predictors, T0)\n",
    "\n",
    "r_series_uber = pd.Series(r_uber).sort_values(ascending=False)\n",
    "p_value_uber = float((r_series_uber >= r_uber[\"Austin\"]).mean())\n",
    "print(f\"Austin's rank: {list(r_series_uber.index).index('Austin') + 1} of {len(r_series_uber)}\")\n",
    "print(f\"p-value: {p_value_uber:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86b536aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(9, 4.5))\n",
    "colors = [\"#8e44ad\" if u == \"Austin\" else \"#2980b9\" for u in r_series_uber.index]\n",
    "ax.scatter(r_series_uber.index, r_series_uber.values, c=colors)\n",
    "ax.set_xticks(range(len(r_series_uber)))\n",
    "ax.set_xticklabels(r_series_uber.index, rotation=45, ha=\"right\")\n",
    "ax.set_ylabel(r\"$r_j = RMSPE_{post}/RMSPE_{pre}$\")\n",
    "ax.set_title(f\"Distribution of $r_j$ (placebo) - Uber/Austin (p = {p_value_uber:.2f})\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cad32e06",
   "metadata": {},
   "source": [
    "Austin lands near the upper tail of the placebo distribution — consistent with the simulated 10% effect and with what Dashboard 3 shows at `true_effect=10%` under a reasonable fit cutoff.\n",
    "\n",
    "### On Proposition 99 — closing the loop with ADH's classic placebo distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310695b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def prop99_predictors(unit, donors):\n",
    "    x1 = pivot_prop99.loc[unit, PRE_YEARS].values\n",
    "    x0 = pivot_prop99.loc[donors, PRE_YEARS].values.T\n",
    "    return x1, x0\n",
    "\n",
    "all_states = [\"California\"] + donor_states\n",
    "gaps_prop99, rmspe_pre_prop99_all, r_prop99 = leave_one_out_placebo(\n",
    "    pivot_prop99[ALL_YEARS], all_states, prop99_predictors, T0_PROP99)\n",
    "\n",
    "r_series_prop99 = pd.Series(r_prop99).sort_values(ascending=False)\n",
    "p_value_prop99 = float((r_series_prop99 >= r_prop99[\"California\"]).mean())\n",
    "print(f\"California's rank: {list(r_series_prop99.index).index('California') + 1} of {len(r_series_prop99)}\")\n",
    "print(f\"p-value (all 39 states, no fit cutoff): {p_value_prop99:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb4e5835",
   "metadata": {},
   "source": [
    "We apply the same kind of fit cutoff Dashboard 3 uses (exclude donors whose pre-treatment RMSPE exceeds $k$ times California's own):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c14e1493",
   "metadata": {},
   "outputs": [],
   "source": [
    "rmspe_pre_ca = rmspe_pre_prop99_all[\"California\"]\n",
    "for k in [2, 3, 5, float(\"inf\")]:\n",
    "    kept = [s for s in all_states if rmspe_pre_prop99_all[s] <= k * rmspe_pre_ca]\n",
    "    r_kept = pd.Series({s: r_prop99[s] for s in kept}).sort_values(ascending=False)\n",
    "    p = float((r_kept >= r_prop99[\"California\"]).mean())\n",
    "    label = \"all\" if k == float(\"inf\") else f\"{k}x\"\n",
    "    print(f\"cutoff={label:>6}: donors kept={len(kept):>2}  p-value={p:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d60a0217",
   "metadata": {},
   "source": [
    "**Unlike Dashboard 3's clean simulation, the cutoff here doesn't necessarily sharpen the p-value.** Missouri and Virginia — the two states with a higher $r_j$ than California — actually fit *better* than California in the pre-treatment period (their pre-treatment RMSPE is 0.27x and 0.49x California's), so no \"exclude poor fit\" rule will remove them from the calculation, no matter how strict. They are states that fit well but had their own large, idiosyncratic move in the post-treatment period — a reminder that the RMSPE cutoff filters out placebos that never fit well to begin with, but doesn't guarantee excluding every large post-treatment gap if it comes from a well-fitting donor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10a8599",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(9, 5))\n",
    "for state in all_states:\n",
    "    color = \"#8e44ad\" if state == \"California\" else \"#bdc3c7\"\n",
    "    lw = 2.2 if state == \"California\" else 0.8\n",
    "    zorder = 5 if state == \"California\" else 1\n",
    "    ax.plot(ALL_YEARS, gaps_prop99[state], color=color, linewidth=lw, zorder=zorder)\n",
    "ax.axhline(0, color=\"black\", linewidth=0.8)\n",
    "ax.axvline(1988.5, color=\"gray\", linestyle=\":\")\n",
    "ax.set_xlabel(\"Year\"); ax.set_ylabel(\"Gap (actual - synthetic)\")\n",
    "ax.set_title(\"Placebo distribution: California vs. the 38 donor states\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87445d92",
   "metadata": {},
   "source": [
    "This is, in essence, ADH (2010)'s Figure 4: California (thick line) departs much more sharply than most placebos toward the end of the post-treatment period.\n",
    "\n",
    "### Contrast: `scpi_pkg`'s prediction intervals\n",
    "\n",
    "`scpi_pkg` implements its own inference via `scpi()` — but this is a **materially different** approach, not a reproduction of the permutation p-value above. Instead of a placebo-based ranking, `scpi()` builds **prediction intervals** (Cattaneo, Feng & Titiunik, 2021) using a combination of simulation and bias/variance adjustment of $\\hat Y_{1t}^N$. The two approaches answer related but non-interchangeable questions: one asks \"how unusual is this gap compared to placebos?\"; the other asks \"what is the plausible range for the unobserved outcome, given estimation uncertainty?\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a7d1ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_prop99 = scpi(data_prop99, w_constr={\"name\": \"simplex\"}, sims=200, cores=1, verbose=False)\n",
    "ci = pi_prop99.CI_all_gaussian\n",
    "ci.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7665ad2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(ALL_YEARS, actual_prop99, label=\"California\", color=\"#8e44ad\", linewidth=2)\n",
    "ax.plot(ALL_YEARS, synthetic_prop99, label=\"Synthetic control\", color=\"#e74c3c\", linestyle=\"--\", linewidth=2)\n",
    "ax.fill_between(POST_YEARS, ci[\"Lower\"].values, ci[\"Upper\"].values,\n",
    "                 color=\"#e74c3c\", alpha=0.15, label=\"Prediction interval (scpi)\")\n",
    "ax.axvline(1988.5, color=\"gray\", linestyle=\":\")\n",
    "ax.set_xlabel(\"Year\"); ax.set_ylabel(\"Cigarette packs per capita\")\n",
    "ax.set_title(\"Proposition 99 with scpi_pkg's prediction interval\")\n",
    "ax.legend(); plt.tight_layout(); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67b58597",
   "metadata": {},
   "source": [
    "## For further exploration\n",
    "\n",
    "- Change the fit cutoff (`k`) in Prop 99's permutation inference section and see how the p-value shifts as more poorly-fitting donors get excluded.\n",
    "- Try restricting the from-scratch Uber/Austin pass to demographics only (dropping the 4 lagged weeks) and compare the resulting weights to Dashboard 1's under `predictor_set=\"demographics\"`.\n",
    "- Adjust `w_constr` in `scest`/`scpi` (e.g. `{\"name\": \"lasso\"}` instead of `\"simplex\"`) and see how the weights' sparsity changes.\n",
    "- With the West Germany data already downloaded under `data/scpi_germany.csv` (Abadie, Diamond & Hainmueller, 2015), replicate the same pipeline as an extra exercise."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}