{
 "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-l3-0",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 10.3 · Fuzzy RDD, Manipulation Tests and Regression Kinks\n",
    "\n",
    "*Accelerated Foundation in Applied Causal Inference*\n",
    "\n",
    "**By the end of this lesson you can:**\n",
    "\n",
    "- Distinguish sharp from fuzzy designs and estimate a fuzzy RDD as IV at the cutoff\n",
    "- Test for manipulation of the running variable with a density test (McCrary / rddensity)\n",
    "- Interpret a density test carefully on real data with heaping and mass points\n",
    "- Recognise a regression kink design, where the slope of a policy changes at a threshold\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-l3-1",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fuzzy RDD: eligibility is not participation"
   ]
  },
  {
   "id": "m10-l3-2",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "!pip install -q numpy pandas rdrobust rddensity matplotlib causaldata\n",
    "import numpy as np, pandas as pd\n",
    "from rdrobust import rdrobust\n",
    "from rddensity import rddensity\n",
    "import matplotlib.pyplot as plt\n",
    "rng = np.random.default_rng(1)"
   ]
  },
  {
   "id": "m10-l3-3",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** What is the effect of a social programme on household consumption when eligibility at an income cutoff only partly determines take-up?"
   ]
  },
  {
   "id": "m10-l3-4",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "N = 6_000\n",
    "income = rng.uniform(-1, 1, N)                       # centred: eligible if income < 0\n",
    "eligible = (income < 0).astype(int)\n",
    "p_take = np.where(eligible == 1, 0.70, 0.10)         # take-up jumps by 60 points\n",
    "takeup = rng.binomial(1, p_take)\n",
    "consumption = 50 - 8 * income + 6 * takeup + rng.normal(0, 4, N)\n",
    "\n",
    "# Itemise the two jumps, then their ratio; rdrobust's fuzzy= option does it with correct inference\n",
    "itt = rdrobust(y=consumption, x=-income, c=0)       # flip sign so eligibility is on the right\n",
    "fs = rdrobust(y=takeup, x=-income, c=0)\n",
    "fz = rdrobust(y=consumption, x=-income, c=0, fuzzy=takeup)\n",
    "print(f\"Jump in consumption (ITT):  {itt.coef.iloc[0, 0]:.2f}\")\n",
    "print(f\"Jump in take-up (1st stage): {fs.coef.iloc[0, 0]:.2f}\")\n",
    "print(f\"Ratio:                       {itt.coef.iloc[0, 0] / fs.coef.iloc[0, 0]:.2f}\")\n",
    "print(f\"Fuzzy RD (rdrobust):         {fz.coef.iloc[0, 0]:.2f}, robust CI \"\n",
    "      f\"[{fz.ci.iloc[2, 0]:.2f}, {fz.ci.iloc[2, 1]:.2f}]  (truth 6)\")"
   ]
  },
  {
   "id": "m10-l3-5",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manipulation and the density test"
   ]
  },
  {
   "id": "m10-l3-6",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Question it answers:** Did firms manipulate their distress scores to qualify for the grant?"
   ]
  },
  {
   "id": "m10-l3-7",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "def simulate_scores(manipulate, n=10_000):\n",
    "    s = rng.normal(50, 12, n)\n",
    "    if manipulate:   # 30% of firms scoring 45-50 nudge themselves just above 50\n",
    "        movers = (s > 45) & (s < 50) & (rng.uniform(size=n) < 0.3)\n",
    "        s[movers] = 50 + rng.uniform(0, 1.5, movers.sum())\n",
    "    return s\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(13, 4), sharey=True)\n",
    "for ax, manip in zip(axes, [False, True]):\n",
    "    s = simulate_scores(manip)\n",
    "    test = rddensity(X=s, c=50)\n",
    "    ax.hist(s, bins=np.arange(20, 81, 1), color=\"#94a3b8\")\n",
    "    ax.axvline(50, color=\"red\", ls=\"--\")\n",
    "    ax.set(title=f\"{'Manipulated' if manip else 'No manipulation'}: density test p = {test.test['p_jk']:.4f}\")\n",
    "plt.tight_layout(); plt.show()"
   ]
  },
  {
   "id": "m10-l3-8",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real data: Uruguay's PANES transfer programme"
   ]
  },
  {
   "id": "m10-l3-9",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "from causaldata import gov_transfers, gov_transfers_density\n",
    "g = gov_transfers.load_pandas().data.dropna()\n",
    "r = rdrobust(y=g.Support, x=-g.Income_Centered, c=0)      # eligibility (income below 0) on the right\n",
    "print(f\"Jump in government support at eligibility: {r.coef.iloc[0, 0]:.3f}, \"\n",
    "      f\"robust CI [{r.ci.iloc[2, 0]:.3f}, {r.ci.iloc[2, 1]:.3f}]\")\n",
    "\n",
    "dens = gov_transfers_density.load_pandas().data.Income_Centered.to_numpy()\n",
    "test = rddensity(X=dens, c=0)\n",
    "print(f\"Density test on {len(dens):,} scores: p = {test.test['p_jk']:.2g}\")\n",
    "plt.hist(dens[np.abs(dens) < 0.05], bins=100, color=\"#94a3b8\"); plt.axvline(0, color=\"red\", ls=\"--\")\n",
    "plt.title(\"PANES income score near the cutoff\"); plt.show()"
   ]
  },
  {
   "id": "m10-l3-10",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression kink designs"
   ]
  },
  {
   "id": "m10-l3-11",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "n_k = 40_000\n",
    "earnings = rng.uniform(20, 80, n_k)                       # past earnings, $000\n",
    "benefit = np.where(earnings < 50, 0.5 * earnings, 25.0)  # replacement rate 50% up to a $50k cap\n",
    "duration = 10 + 0.4 * benefit + 0.05 * earnings + rng.normal(0, 2, n_k)   # weeks unemployed\n",
    "\n",
    "slope_y = rdrobust(y=duration, x=earnings, c=50, deriv=1)   # change in slope of the outcome\n",
    "slope_b = -0.5                                             # known change in the benefit slope at the kink\n",
    "print(f\"Kink in duration slope: {slope_y.coef.iloc[0, 0]:.3f} | implied effect of $1k more benefit: \"\n",
    "      f\"{slope_y.coef.iloc[0, 0] / slope_b:.2f} weeks (truth 0.40)\")"
   ]
  },
  {
   "id": "m10-l3-12",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ✍️ Exercise\n",
    "\n",
    "In the fuzzy simulation, shrink the take-up jump from 60 to 15 percentage points (for example, 0.25 vs\n",
    "0.10). What happens to the fuzzy RD confidence interval? Connect your answer to weak instruments in Module 9."
   ]
  },
  {
   "id": "m10-l3-13",
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Your answer here"
   ]
  },
  {
   "id": "m10-l3-14",
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checklist and recap"
   ]
  }
 ]
}