{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" }, "colab": { "provenance": [] } }, "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 74: Local Polynomial Regression & Kernels\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "RD estimation means estimating a regression function at the edge of its support. Kernel weights, local constant versus local linear fits, why local constant estimates are biased at a boundary while local linear ones are not, and how bias and variance scale with the bandwidth.\n", "\n", "**Learning outcomes**\n", "\n", "- Compute uniform, triangular and Epanechnikov kernel weights.\n", "- Estimate a boundary value by kernel-weighted least squares (local constant and local linear).\n", "- Derive the order of the boundary bias for local constant and local linear fits, and the variance's dependence on n and h.\n", "- Apply local polynomial methods to birth-weight and scholarship cutoffs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Kernels\n", "\n", "A **kernel** $K(u)$ weights an observation by its distance from the evaluation point $c$, scaled by a **bandwidth** $h$: weight $K((X_i - c)/h)$. Common choices on $|u| \\le 1$:\n", "\n", "| Kernel | $K(u)$ |\n", "|---|---|\n", "| Uniform | $1/2$ |\n", "| Triangular | $1 - \\lvert u\\rvert$ |\n", "| Epanechnikov | $\\tfrac34(1 - u^2)$ |\n", "\n", "Observations beyond $h$ get weight 0. The triangular kernel is the usual default for RD because it is optimal (in mean squared error) for estimating at a boundary point.\n", "\n", "### Local polynomial regression\n", "\n", "On one side of the cutoff, solve the weighted least-squares problem\n", "\n", "$$\n", "\\min_{b_0, \\dots, b_p}\\sum_i K\\Big(\\frac{X_i - c}{h}\\Big)\\Big(Y_i - b_0 - b_1(X_i - c) - \\dots - b_p(X_i - c)^p\\Big)^2,\n", "$$\n", "\n", "and use $\\hat b_0$ as the estimate of $\\mu(c^+)$ (or $\\mu(c^-)$). The RD estimate is the difference of the two intercepts. In matrix form, $\\hat{\\mathbf b} = (\\mathbf X^\\top\\mathbf W\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf W\\mathbf y$ with $\\mathbf W = \\operatorname{diag}(K_i)$.\n", "\n", "- $p = 0$: **local constant** (Nadaraya–Watson): a weighted mean.\n", "- $p = 1$: **local linear**: the RD workhorse.\n", "- $p = 2$: local quadratic: used for bias correction (Day 75).\n", "\n", "### Boundary bias\n", "\n", "At an interior point, local constant and local linear have bias of order $h^2$. At a **boundary** (data on one side only), the local constant estimator has bias of order $h$: if the regression function slopes, the weighted mean is pulled towards the values inside the window. Local linear removes the first-order term:\n", "\n", "$$\n", "\\text{bias}_{LL} \\approx \\frac{h^2}{2}\\mu''(c)\\,B_K, \\qquad \\operatorname{Var}_{LL} \\approx \\frac{\\sigma^2(c)}{n f(c) h}V_K,\n", "$$\n", "\n", "where $f(c)$ is the density of $X$ at $c$ and $B_K, V_K$ are constants depending on the kernel.\n", "\n", "### Bias–variance trade-off\n", "\n", "A wider bandwidth uses more observations (variance falls like $1/(nh)$) but lets curvature in (bias grows like $h^2$). Choosing $h$ is Day 75's topic. The kernel shape matters much less than the bandwidth.\n", "\n", "### Plots are not estimates\n", "\n", "RD plots show binned means of $Y$ across $X$ with a fitted curve. They are for visual inspection; the estimate comes from the local polynomial fit at the cutoff.\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · Triangular weights\n", "\n", "With $h = 2$, compute triangular weights for observations at distances 0.5, 1, 1.5 and 2.5 from the cutoff.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $K(u) = 1 - |u|$ with $u = d/h$.\n", "\n", "**Step 2.** $1 - 0.25 = 0.75$; $1 - 0.5 = 0.5$; $1 - 0.75 = 0.25$; distance 2.5 > h: 0.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · A local constant estimate\n", "\n", "Right of the cutoff, observations $(x - c, y)$: (0.5, 12), (1.0, 14), (1.5, 16), which lie exactly on $y = 10 + 4(x - c)$. With the weights from Exercise 1, compute the local constant estimate of $\\mu(c^+)$ and its error.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Weighted mean: $\\dfrac{0.75(12) + 0.5(14) + 0.25(16)}{0.75 + 0.5 + 0.25} = \\dfrac{9 + 7 + 4}{1.5} = 13.33$.\n", "\n", "**Step 2.** True value at the cutoff: 10.\n", "\n", "**Answer:** 13.33, an error of +3.33 from the slope (boundary bias). A local linear fit recovers 10 exactly (Exercise 7).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · Epanechnikov weights\n", "\n", "Compute Epanechnikov weights at $u = 0, 0.5, 0.9$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\tfrac34(1 - 0) = 0.75$; $\\tfrac34(1 - 0.25) = 0.5625$; $\\tfrac34(1 - 0.81) = 0.1425$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · Halving the bandwidth\n", "\n", "If the bandwidth is halved, how do the leading bias and variance of a local linear estimator change?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Bias $\\propto h^2$: divided by 4.\n", "\n", "**Step 2.** Variance $\\propto 1/(nh)$: doubled.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · Linear or quadratic?\n", "\n", "Why not always use a local quadratic fit, which has smaller bias?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Each additional polynomial term increases variance at the boundary, where few observations pin down the intercept.\n", "\n", "**Answer:** local linear is the usual compromise; local quadratic is used mainly to estimate and remove the linear fit's bias (Day 75).\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · Local constant bias is of order h\n", "\n", "Let $\\mu(x) = a + b(x - c)$ for $x \\ge c$, with $X - c$ uniform on $[0, h]$ and a uniform kernel. Show that the local constant estimator converges to $a + bh/2$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** The local constant estimator is the mean of $Y$ over the window, which converges to $E[\\mu(X)] = a + bE[X - c]$.\n", "\n", "**Step 2.** $E[X - c] = h/2$ for a uniform distribution on $[0, h]$.\n", "\n", "**Answer:** bias $bh/2$, of order $h$ whenever $b \\ne 0$. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · Local linear is exact for lines\n", "\n", "Show that if $\\mu(x)$ is exactly linear on the window and there is no noise, weighted least squares with any positive weights recovers the intercept exactly.\n", "\n", "
Show solution\n", "\n", "**Step 1.** The data satisfy $\\mathbf y = \\mathbf X\\mathbf b^*$ exactly.\n", "\n", "**Step 2.** $\\hat{\\mathbf b} = (\\mathbf X^\\top\\mathbf W\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf W\\mathbf X\\mathbf b^* = \\mathbf b^*$. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · Local linear bias with curvature\n", "\n", "Let $\\mu(x) = a + b(x - c) + q(x - c)^2$ with $X - c$ uniform on $[0, h]$ and a uniform kernel. Show the local linear intercept converges to $a - qh^2/6$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** The fitted intercept is $a$ plus $q$ times the intercept of the best linear approximation to $u^2$ on $[0, h]$.\n", "\n", "**Step 2.** Projection of $u^2$ on $(1, u)$ for $u \\sim U[0, h]$: slope $\\operatorname{Cov}(u, u^2)/\\operatorname{Var}(u) = (h^3/12)/(h^2/12) = h$; intercept $E[u^2] - hE[u] = h^2/3 - h^2/2 = -h^2/6$.\n", "\n", "**Step 3.** So the intercept converges to $a - qh^2/6$: bias of order $h^2$. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · Variance falls like 1/(nh)\n", "\n", "With a uniform kernel, $X$ having density $f$ near $c$ and constant variance $\\sigma^2$, argue that the variance of the local constant estimator is about $\\sigma^2/(nf(c)h)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** The expected number of observations in $[c, c + h]$ is about $nf(c)h$.\n", "\n", "**Step 2.** The estimator is their mean, with variance $\\sigma^2/(\\text{number of observations})$.\n", "\n", "**Answer:** $\\approx \\sigma^2/(nf(c)h)$; local linear has the same order with a larger kernel constant. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Weighted least squares\n", "\n", "Derive $\\hat{\\mathbf b} = (\\mathbf X^\\top\\mathbf W\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf W\\mathbf y$ from the weighted objective.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Objective: $(\\mathbf y - \\mathbf X\\mathbf b)^\\top\\mathbf W(\\mathbf y - \\mathbf X\\mathbf b)$.\n", "\n", "**Step 2.** Gradient: $-2\\mathbf X^\\top\\mathbf W(\\mathbf y - \\mathbf X\\mathbf b) = \\mathbf 0$.\n", "\n", "**Step 3.** Solve. $\\blacksquare$ Equivalently, OLS on $\\sqrt{K_i}$-scaled data.\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · A boundary estimate by hand\n", "\n", "Mortality data right of a birth-weight cutoff (very-low-birth-weight babies just above 1,500 g get less intensive care): distances in hundreds of grams and one-year mortality per 1,000: (0.2, 21), (0.5, 22), (1.0, 25), (1.5, 26). With a triangular kernel and $h = 2$, give the weights and the local constant and local linear estimates of $\\mu(c^+)$ (the code computes them).\n", "\n", "
Show solution\n", "\n", "**Step 1.** Weights: $1 - d/2$: 0.9, 0.75, 0.5, 0.25.\n", "\n", "**Step 2.** Local constant: $\\dfrac{0.9(21) + 0.75(22) + 0.5(25) + 0.25(26)}{2.4} = \\dfrac{54.4}{2.4} \\approx 22.67$.\n", "\n", "**Step 3.** Local linear (weighted regression on distance): intercept about 20.10, slope about 4.30 per 100 g.\n", "\n", "**Answer:** the local constant estimate is pulled up by 2.6 deaths per 1,000 because mortality rises with distance from the cutoff on this side; local linear corrects it.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Why the slope matters in health RD\n", "\n", "In the very-low-birth-weight design, mortality decreases steeply with birth weight on both sides of 1,500 g. Explain how a local constant RD would bias the effect of extra care.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Below the cutoff (lighter babies, more care), the window mean averages over lighter, higher-mortality babies: the left limit is overstated.\n", "\n", "**Step 2.** Above, the window mean averages over heavier, lower-mortality babies: the right limit is understated.\n", "\n", "**Answer:** the jump (right minus left) is biased downwards, which would exaggerate the benefit of extra care for those below the cutoff. Local linear fits remove the first-order bias.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · Kernel versus bandwidth\n", "\n", "In the simulation below, compare the RD estimate's mean and SD across kernels (uniform, triangular, Epanechnikov) and across bandwidths (0.1, 0.3, 0.6). Which choice matters more?\n", "\n", "
Show solution\n", "\n", "**Step 1.** For a fixed bandwidth, the three kernels give nearly the same bias and SD.\n", "\n", "**Step 2.** Across bandwidths, the SD falls and the bias rises noticeably as $h$ grows.\n", "\n", "**Answer:** bandwidth choice dominates kernel choice.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · Local quadratic for a scholarship\n", "\n", "Admission probability rises sharply and non-linearly with the test score near a scholarship cutoff. When would a local quadratic fit be preferred to local linear?\n", "\n", "
Show solution\n", "\n", "**Step 1.** When curvature near the cutoff is strong relative to the bandwidth that the sample size allows.\n", "\n", "**Step 2.** And when the sample is large enough that the extra variance is acceptable.\n", "\n", "**Answer:** report both as a robustness check; Day 75's bias-corrected inference uses the quadratic fit systematically.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · Reading an RD plot\n", "\n", "A plot shows binned means with a global fourth-order polynomial and a visible jump. A referee asks for \"the estimate in the plot\". What should be reported instead?\n", "\n", "
Show solution\n", "\n", "**Step 1.** The plot's polynomial is for visualisation; its jump is sensitive to observations far from the cutoff.\n", "\n", "**Answer:** report the local linear estimate with a data-driven bandwidth and robust confidence interval, and show the plot with the local fits overlaid.\n", "\n", "
\n", "\n", "## Verify in Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "\n", "K = {\"uniform\": lambda u: .5 * (np.abs(u) <= 1), \"triangular\": lambda u: np.clip(1 - np.abs(u), 0, None),\n", " \"epanechnikov\": lambda u: .75 * np.clip(1 - u ** 2, 0, None)}\n", "def boundary(x, y, h, kern=\"triangular\", p=1):\n", " w = K[kern](x / h); m = w > 0\n", " X = np.column_stack([x[m] ** k for k in range(p + 1)])\n", " W = w[m]\n", " return np.linalg.solve(X.T @ (W[:, None] * X), X.T @ (W * y[m]))[0]\n", "\n", "print(\"Ex 1 weights:\", K[\"triangular\"](np.array([.5, 1, 1.5, 2.5]) / 2))\n", "x = np.array([.5, 1, 1.5]); y = 10 + 4 * x\n", "print(\"Ex 2: local constant\", round(boundary(x, y, 2, p=0), 3), \" local linear\", round(boundary(x, y, 2, p=1), 3))\n", "x = np.array([.2, .5, 1., 1.5]); y = np.array([21, 22, 25, 26.])\n", "print(\"Ex 11: local constant\", round(boundary(x, y, 2, p=0), 3), \" local linear\", round(boundary(x, y, 2, p=1), 3))\n", "\n", "# Ex 8: local linear bias -q h^2/6 with uniform design and kernel\n", "u = np.linspace(0, 1, 200001); print(\"Ex 8 intercept for u^2 on [0,1]:\", round(boundary(u, u ** 2, 1, \"uniform\"), 4), \"(-1/6 = -0.1667)\")\n", "\n", "# Exercise 13: kernels and bandwidths in a simulated RD (truth 5)\n", "rng = np.random.default_rng(74)\n", "mu0 = lambda x: 50 + 8 * x + 6 * x ** 2 - 4 * x ** 3\n", "def rd(x, y, h, kern):\n", " r, l = x >= 0, x < 0\n", " return boundary(x[r], y[r], h, kern) - boundary(-x[l], y[l], h, kern)\n", "for h in (.1, .3, .6):\n", " for kern in K:\n", " est = []\n", " for _ in range(400):\n", " xx = rng.uniform(-1, 1, 1000); yy = mu0(xx) + 5 * (xx >= 0) + rng.normal(0, 3, 1000)\n", " est.append(rd(xx, yy, h, kern))\n", " print(f\"h={h:.1f} {kern:12s} bias {np.mean(est) - 5:6.2f} sd {np.std(est):5.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- Kernels weight observations by distance from the cutoff; the bandwidth sets how far the fit reaches.\n", "- At a boundary, local constant estimates have bias of order h; local linear removes it, leaving order h² bias from curvature.\n", "- Variance scales like 1/(nh); bandwidth matters far more than kernel shape." ] } ] }