{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "id": "m12-l1-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 12.1 · Synthetic Control by Convex Optimisation\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Explain why a single treated unit needs a weighted, data-driven control group\n",
    "- Solve for synthetic control weights with scipy (non-negative, summing to one)\n",
    "- Evaluate pre-treatment fit and read the post-treatment gap as the effect\n",
    "- Explain why the weights rule out extrapolation and why that matters\n",
    "\n",
    "> Read the full explanations on the course website, then run each cell here in order (**Runtime → Run all** works on a fresh runtime)."
   ]
  },
  {
   "id": "m12-l1-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One treated province, many possible controls"
   ]
  },
  {
   "id": "m12-l1-2",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "!pip install -q numpy pandas scipy matplotlib seaborn\n",
    "import numpy as np, pandas as pd\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt, seaborn as sns\n",
    "rng = np.random.default_rng(42)\n",
    "sns.set_theme(style=\"whitegrid\")"
   ]
  },
  {
   "id": "m12-l1-3",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A cap-and-trade programme in Province A"
   ]
  },
  {
   "id": "m12-l1-4",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** How much did cap-and-trade reduce Province A's emissions compared with what would have happened otherwise?"
   ]
  },
  {
   "id": "m12-l1-5",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "years = np.arange(2005, 2025)\n",
    "T0 = 2018\n",
    "pre = years < T0\n",
    "J = 10                                                  # donors\n",
    "\n",
    "factors = np.column_stack([np.ones(len(years)), (years - 2005) / 10, np.sin((years - 2005) / 3)])\n",
    "loadings = rng.uniform([80, 0, -5], [140, 30, 10], size=(J, 3))\n",
    "donors = loadings @ factors.T + rng.normal(0, 1.5, (J, len(years)))         # J x T\n",
    "\n",
    "true_w = np.array([0.4, 0.5, 0.1] + [0] * (J - 3))                          # A resembles donors 0, 1, 2\n",
    "effect = np.where(years >= T0, -15.0, 0.0)\n",
    "treated = true_w @ donors + effect + rng.normal(0, 1.0, len(years))\n",
    "\n",
    "panel = pd.DataFrame(donors.T, index=years, columns=[f\"Donor {j}\" for j in range(J)])\n",
    "panel[\"Province A\"] = treated\n",
    "panel.iloc[:, [0, 1, 2, -1]].plot(figsize=(10, 4), title=\"Emissions: Province A and three donors\"); plt.show()"
   ]
  },
  {
   "id": "m12-l1-6",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving for the weights"
   ]
  },
  {
   "id": "m12-l1-7",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** Which combination of untreated provinces best reproduces Province A's emissions before 2018?"
   ]
  },
  {
   "id": "m12-l1-8",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def synth_weights(y_treated_pre, Y_donors_pre):\n",
    "    J = Y_donors_pre.shape[0]\n",
    "    loss = lambda w: np.sum((y_treated_pre - w @ Y_donors_pre) ** 2)\n",
    "    res = minimize(loss, np.full(J, 1 / J), method=\"SLSQP\",\n",
    "                   bounds=[(0, 1)] * J, constraints={\"type\": \"eq\", \"fun\": lambda w: w.sum() - 1})\n",
    "    return res.x\n",
    "\n",
    "w = synth_weights(treated[pre], donors[:, pre])\n",
    "print(pd.Series(w, index=panel.columns[:J]).round(3).sort_values(ascending=False).head(5))"
   ]
  },
  {
   "id": "m12-l1-9",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit, gap and interpretation"
   ]
  },
  {
   "id": "m12-l1-10",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "synthetic = w @ donors\n",
    "gap = treated - synthetic\n",
    "rmspe_pre = np.sqrt(np.mean(gap[pre] ** 2))\n",
    "print(f\"Pre-treatment RMSPE: {rmspe_pre:.2f}\")\n",
    "print(f\"Average post-treatment gap: {gap[~pre].mean():.2f} (truth -15)\")\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n",
    "axes[0].plot(years, treated, color=\"red\", lw=2, label=\"Province A (actual)\")\n",
    "axes[0].plot(years, synthetic, color=\"blue\", ls=\"--\", lw=2, label=\"Synthetic Province A\")\n",
    "axes[0].axvline(T0, color=\"black\", ls=\":\"); axes[0].legend(); axes[0].set_title(\"Actual vs synthetic\")\n",
    "axes[1].bar(years, gap, color=np.where(pre, \"#94a3b8\", \"#ef4444\"))\n",
    "axes[1].axvline(T0 - 0.5, color=\"black\", ls=\":\"); axes[1].set_title(\"Gap = estimated effect each year\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "id": "m12-l1-11",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "id": "m12-l1-12",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "Compare the synthetic control with two alternatives: (1) the simple **average of all donors** as the\n",
    "control (a DiD-style comparison), and (2) **unconstrained OLS** weights (allowing negative weights and no\n",
    "sum-to-one). Compare pre-period fit and the estimated effect. Why are OLS weights dangerous with few\n",
    "pre-periods?"
   ]
  },
  {
   "id": "m12-l1-13",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m12-l1-14",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}