{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "colab": {
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "cells": [
  {
   "id": "m03-l3-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 3.3 · Statistical Inference and Robust Standard Errors\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Read standard errors, t-statistics, p-values and confidence intervals from a regression\n",
    "- Simulate the false-positive rate of a test to see what a p-value promises\n",
    "- Diagnose heteroskedasticity and fix inference with robust (HC) standard errors\n",
    "- Recognise when observations are clustered and why that matters for policy evaluation\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": "m03-l3-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From estimates to uncertainty"
   ]
  },
  {
   "id": "m03-l3-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(42)\n",
    "sns.set_theme(style=\"whitegrid\")"
   ]
  },
  {
   "id": "m03-l3-3",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What a p-value promises: a placebo simulation"
   ]
  },
  {
   "id": "m03-l3-4",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** If a programme does nothing, how often will a regression still call it \"significant\"?"
   ]
  },
  {
   "id": "m03-l3-5",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def placebo_study(n=200):\n",
    "    d = rng.binomial(1, 0.5, n)                 # a 'programme' with zero effect\n",
    "    y = 50 + rng.normal(0, 10, n)\n",
    "    p = smf.ols(\"y ~ d\", data=pd.DataFrame({\"y\": y, \"d\": d})).fit().pvalues[\"d\"]\n",
    "    return p < 0.05\n",
    "\n",
    "false_positive_rate = np.mean([placebo_study() for _ in range(1_000)])\n",
    "print(f\"Share of placebo studies with p < 0.05: {false_positive_rate:.3f}\")"
   ]
  },
  {
   "id": "m03-l3-6",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heteroskedasticity and robust standard errors"
   ]
  },
  {
   "id": "m03-l3-7",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** Are my standard errors too small because wage dispersion grows with education?"
   ]
  },
  {
   "id": "m03-l3-8",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "N = 1_000\n",
    "education = rng.uniform(8, 20, N)\n",
    "errors = rng.normal(0, 150 * (education - 7) ** 2)  # noise SD grows fast with education\n",
    "wages = 25_000 + 3_000 * education + errors\n",
    "df = pd.DataFrame({\"education\": education, \"wages\": wages})\n",
    "\n",
    "classic = smf.ols(\"wages ~ education\", data=df).fit()\n",
    "robust = smf.ols(\"wages ~ education\", data=df).fit(cov_type=\"HC3\")\n",
    "print(f\"Coefficient (identical in both): {classic.params['education']:,.0f}\")\n",
    "print(f\"Classical SE: {classic.bse['education']:,.1f} | Robust HC3 SE: {robust.bse['education']:,.1f}\")\n",
    "print(f\"Classical SEs understate the uncertainty by {100 * (1 - classic.bse['education'] / robust.bse['education']):.0f}%\")\n",
    "\n",
    "plt.figure(figsize=(9, 4))\n",
    "plt.scatter(classic.fittedvalues, classic.resid, alpha=.4)\n",
    "plt.axhline(0, color=\"red\"); plt.title(\"Residuals fan out: heteroskedasticity\")\n",
    "plt.xlabel(\"Fitted wage\"); plt.ylabel(\"Residual\"); plt.show()"
   ]
  },
  {
   "id": "m03-l3-9",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "from statsmodels.stats.diagnostic import het_breuschpagan\n",
    "lm, lm_p, _, _ = het_breuschpagan(classic.resid, classic.model.exog)\n",
    "print(f\"Breusch-Pagan p-value: {lm_p:.2g}  (small -> heteroskedasticity)\")"
   ]
  },
  {
   "id": "m03-l3-10",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clustered data: when observations are not independent"
   ]
  },
  {
   "id": "m03-l3-11",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** If a programme is assigned to whole schools, how wrong are standard errors that ignore it?"
   ]
  },
  {
   "id": "m03-l3-12",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def clustered_study(n_schools=40, pupils=25, cluster_se=False):\n",
    "    school = np.repeat(np.arange(n_schools), pupils)\n",
    "    treated_school = rng.permutation(np.repeat([0, 1], n_schools // 2))[school]  # half the schools\n",
    "    school_shock = rng.normal(0, 5, n_schools)[school]            # shared within a school\n",
    "    y = 50 + school_shock + rng.normal(0, 5, school.size)         # NO true effect\n",
    "    d = pd.DataFrame({\"y\": y, \"d\": treated_school, \"school\": school})\n",
    "    if cluster_se:\n",
    "        fit = smf.ols(\"y ~ d\", data=d).fit(cov_type=\"cluster\", cov_kwds={\"groups\": d[\"school\"]}, use_t=True)\n",
    "    else:\n",
    "        fit = smf.ols(\"y ~ d\", data=d).fit(cov_type=\"HC1\")\n",
    "    return fit.pvalues[\"d\"] < 0.05\n",
    "\n",
    "for flag in (False, True):\n",
    "    rate = np.mean([clustered_study(cluster_se=flag) for _ in range(500)])\n",
    "    print(f\"Cluster-robust SE = {flag!s:5}: false-positive rate = {rate:.3f}\")"
   ]
  },
  {
   "id": "m03-l3-13",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "id": "m03-l3-14",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "With only **8 schools**, compare cluster-robust standard errors that use normal critical values\n",
    "(`use_t=False`) with those that use the small-sample t correction (`use_t=True`, as above). Which keeps\n",
    "the false-positive rate near 5%? What does that tell you about evaluating a policy rolled out in only a\n",
    "handful of provinces?"
   ]
  },
  {
   "id": "m03-l3-15",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m03-l3-16",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}