{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "id": "m10-l1-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 10.1 · Sharp Regression Discontinuity: Grants at a Distress Threshold\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Explain why an administrative cutoff creates a local experiment\n",
    "- State the continuity assumption and what would violate it\n",
    "- Estimate a sharp RDD by local linear regression with different slopes on each side\n",
    "- Draw a binned RD plot that lets the reader see the jump\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": "m10-l1-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholds everywhere"
   ]
  },
  {
   "id": "m10-l1-2",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "!pip install -q numpy pandas statsmodels matplotlib seaborn\n",
    "import numpy as np, pandas as pd\n",
    "import statsmodels.formula.api as smf\n",
    "import matplotlib.pyplot as plt, seaborn as sns\n",
    "rng = np.random.default_rng(2026)\n",
    "sns.set_theme(style=\"whitegrid\")"
   ]
  },
  {
   "id": "m10-l1-3",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A municipal grant at a distress score of 50"
   ]
  },
  {
   "id": "m10-l1-4",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** Does a small-business grant raise firm revenue, judged by firms just either side of the eligibility cutoff?"
   ]
  },
  {
   "id": "m10-l1-5",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "N, CUTOFF = 3_000, 50\n",
    "score = rng.uniform(20, 80, N)\n",
    "grant = (score >= CUTOFF).astype(int)\n",
    "revenue = 100_000 - 500 * score + 4 * (score - 50) ** 2 + 12_000 * grant + rng.normal(0, 6_000, N)\n",
    "df = pd.DataFrame({\"score\": score, \"grant\": grant, \"revenue\": revenue})\n",
    "df[\"x\"] = df.score - CUTOFF                         # centre the running variable at the cutoff\n",
    "\n",
    "naive = df.revenue[df.grant == 1].mean() - df.revenue[df.grant == 0].mean()\n",
    "print(f\"Naive difference: {naive:,.0f} (true effect at the cutoff: 12,000)\")"
   ]
  },
  {
   "id": "m10-l1-6",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local linear estimation"
   ]
  },
  {
   "id": "m10-l1-7",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** How large is the jump in revenue at a distress score of exactly 50?"
   ]
  },
  {
   "id": "m10-l1-8",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def local_linear(data, h, kernel=\"triangular\"):\n",
    "    sub = data[data.x.abs() <= h].copy()\n",
    "    w = (1 - sub.x.abs() / h) if kernel == \"triangular\" else np.ones(len(sub))\n",
    "    fit = smf.wls(\"revenue ~ grant * x\", data=sub, weights=w).fit(cov_type=\"HC1\")\n",
    "    return fit.params[\"grant\"], fit.bse[\"grant\"], len(sub)\n",
    "\n",
    "for h in [5, 10, 15, 30]:\n",
    "    est, se, n = local_linear(df, h)\n",
    "    print(f\"h = {h:>2}: effect = {est:>8,.0f} (SE {se:,.0f}), n = {n:,}\")"
   ]
  },
  {
   "id": "m10-l1-9",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The RD plot"
   ]
  },
  {
   "id": "m10-l1-10",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "h = 15\n",
    "est, se, _ = local_linear(df, h)\n",
    "bins = pd.cut(df.score, np.arange(20, 81, 2))\n",
    "binned = df.groupby(bins, observed=True).agg(score=(\"score\", \"mean\"), revenue=(\"revenue\", \"mean\"))\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.scatter(binned.score, binned.revenue, color=\"#334155\", zorder=3, label=\"Bin means\")\n",
    "for side in (0, 1):\n",
    "    sub = df[(df.grant == side) & (df.x.abs() <= h)]\n",
    "    line = smf.ols(\"revenue ~ x\", data=sub).fit()\n",
    "    xs = np.linspace(sub.x.min(), sub.x.max(), 50)\n",
    "    plt.plot(xs + CUTOFF, line.params[\"Intercept\"] + line.params[\"x\"] * xs, lw=3,\n",
    "             color=\"#f97316\" if side else \"#2563eb\", label=\"Grant\" if side else \"No grant\")\n",
    "plt.axvline(CUTOFF, ls=\"--\", color=\"red\")\n",
    "plt.title(f\"Sharp RDD: estimated jump = ${est:,.0f} (SE {se:,.0f})\")\n",
    "plt.xlabel(\"Distress score\"); plt.ylabel(\"Revenue one year later ($)\"); plt.legend(); plt.show()"
   ]
  },
  {
   "id": "m10-l1-11",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: the danger of global polynomials"
   ]
  },
  {
   "id": "m10-l1-12",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "A common older practice fits a high-order polynomial to the whole sample instead of local lines. Estimate\n",
    "the jump with global polynomials of order 1 to 5 (interacted with `grant`) using all observations. How\n",
    "stable are the estimates? (Gelman and Imbens, 2019, argue against orders above 2.)"
   ]
  },
  {
   "id": "m10-l1-13",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m10-l1-14",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}