{ "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 29: OLS Asymptotics & Robust (Sandwich) Variance\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "Large-sample theory for OLS without normality or homoskedasticity, and the heteroskedasticity- and cluster-robust standard errors applied work relies on.\n", "\n", "**Learning outcomes**\n", "\n", "- Prove consistency of OLS using the law of large numbers and the continuous mapping theorem.\n", "- Derive the asymptotic sandwich variance and show when it collapses to σ²Q⁻¹.\n", "- Compute HC0/HC1 and cluster-robust variances and explain the Moulton problem.\n", "- Run Wald tests of joint hypotheses with a robust covariance matrix." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Assumptions for large samples\n", "\n", "Drop normality and homoskedasticity. Assume iid observations $(y_i, x_i)$ with\n", "\n", "$$\n", "y_i = x_i^{\\top}\\beta + \\varepsilon_i, \\qquad E[x_i\\varepsilon_i] = 0, \\qquad Q = E[x_ix_i^{\\top}] \\text{ invertible}.\n", "$$\n", "\n", "### Consistency\n", "\n", "$$\n", "\\hat\\beta - \\beta = \\Big(\\tfrac1n\\sum_ix_ix_i^{\\top}\\Big)^{-1}\\Big(\\tfrac1n\\sum_ix_i\\varepsilon_i\\Big) \\xrightarrow{p} Q^{-1}\\cdot 0 = 0,\n", "$$\n", "\n", "by the LLN applied to each average and the **continuous mapping theorem** (continuous functions preserve convergence in probability).\n", "\n", "### Asymptotic normality: the sandwich\n", "\n", "By the CLT, $\\frac{1}{\\sqrt n}\\sum x_i\\varepsilon_i \\xrightarrow{d} N(0, \\Omega)$ with $\\Omega = E[\\varepsilon_i^2x_ix_i^{\\top}]$. By **Slutsky's theorem**,\n", "\n", "$$\n", "\\sqrt n(\\hat\\beta - \\beta) \\xrightarrow{d} N\\big(0,\\; Q^{-1}\\Omega Q^{-1}\\big).\n", "$$\n", "\n", "\"Bread–meat–bread\". Its sample analogue is the **heteroskedasticity-robust (HC0, White)** estimator\n", "\n", "$$\n", "\\widehat{\\operatorname{Var}}(\\hat\\beta) = (X^{\\top}X)^{-1}\\Big(\\sum_ie_i^2x_ix_i^{\\top}\\Big)(X^{\\top}X)^{-1},\n", "$$\n", "\n", "with HC1 multiplying by $n/(n - k)$ (Stata's `robust`). Under homoskedasticity, $E[\\varepsilon^2 \\mid x] = \\sigma^2$, so $\\Omega = \\sigma^2Q$ and the sandwich collapses to $\\sigma^2Q^{-1}$, the classical formula.\n", "\n", "### Clustered data\n", "\n", "If errors are correlated within clusters $g = 1, \\dots, G$ but independent across them,\n", "\n", "$$\n", "\\widehat{\\operatorname{Var}}_{CR}(\\hat\\beta) = (X^{\\top}X)^{-1}\\Big(\\sum_gX_g^{\\top}e_ge_g^{\\top}X_g\\Big)(X^{\\top}X)^{-1}.\n", "$$\n", "\n", "Its reliability depends on the number of **clusters** $G$, not the number of observations. With few clusters (say under 30–50), use small-sample corrections or the wild cluster bootstrap.\n", "\n", "> **Causal link.** Treatment is often assigned at the cluster level (schools, villages, states). Inference must then be clustered at the level of assignment, or standard errors can be badly understated.\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · A scalar sandwich\n", "\n", "In a no-intercept regression $y = \\beta x + \\varepsilon$, write the classical and the HC0 variance of $\\hat\\beta$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Classical: $s^2/\\sum x_i^2$.\n", "\n", "**Step 2.** HC0: $\\sum x_i^2e_i^2/(\\sum x_i^2)^2$.\n", "\n", "**Step 3.** The robust version weights each squared residual by its $x_i^2$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · Numbers\n", "\n", "Take $x = (1, 2, 3)$ and residuals $e = (0.5, -1, 0.5)$ (no intercept, $k = 1$). Compute both variances.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\sum x^2 = 1 + 4 + 9 = 14$.\n", "\n", "**Step 2.** Classical: $s^2 = (0.25 + 1 + 0.25)/(3 - 1) = 0.75$, so $\\operatorname{Var} = 0.75/14 \\approx 0.0536$.\n", "\n", "**Step 3.** HC0: $\\sum x^2e^2 = 1(0.25) + 4(1) + 9(0.25) = 6.5$, so $\\operatorname{Var} = 6.5/14^2 = 6.5/196 \\approx 0.0332$.\n", "\n", "**Answer:** 0.0536 (classical) versus 0.0332 (robust). (These residuals satisfy $\\sum xe = 0.5 - 2 + 1.5 = 0$, as OLS residuals must.)\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · HC1\n", "\n", "Apply the HC1 correction to Exercise 2.\n", "\n", "
Show solution\n", "\n", "$n/(n - k) = 3/2$: $0.0332 \\times 1.5 \\approx 0.0497$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · When do they agree?\n", "\n", "Under what condition do classical and robust SEs agree in large samples, and what does a large gap suggest?\n", "\n", "
Show solution\n", "\n", "Under homoskedasticity (more precisely, when $\\varepsilon^2$ is uncorrelated with $x_ix_i^{\\top}$). A large gap signals heteroskedasticity, or sometimes model misspecification.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · Counting clusters\n", "\n", "A study has 10,000 pupils in 20 schools, with treatment assigned by school. How much information is there for inference on the treatment effect?\n", "\n", "
Show solution\n", "\n", "Effectively 20 clusters (10 treated, 10 control). Cluster-robust inference treats each school as one draw; with so few, use a $t_{G-1}$ reference distribution or the wild cluster bootstrap.\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · Consistency\n", "\n", "Prove $\\hat\\beta \\xrightarrow{p} \\beta$ under the assumptions above.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Substitute $y_i = x_i^{\\top}\\beta + \\varepsilon_i$ into $\\hat\\beta$:\n", "\n", "$$\n", "\\hat\\beta - \\beta = \\Big(\\tfrac1n\\textstyle\\sum_ix_ix_i^{\\top}\\Big)^{-1}\\Big(\\tfrac1n\\sum_ix_i\\varepsilon_i\\Big).\n", "$$\n", "\n", "**Step 2.** LLN, applied entry by entry: $\\frac1n\\sum x_ix_i^{\\top} \\xrightarrow{p} Q$.\n", "\n", "**Step 3.** LLN again: $\\frac1n\\sum x_i\\varepsilon_i \\xrightarrow{p} E[x\\varepsilon] = 0$.\n", "\n", "**Step 4.** Matrix inversion is continuous at the invertible $Q$, so by the continuous mapping theorem the product converges to $Q^{-1}\\cdot 0 = 0$.\n", "\n", "**Step 5.** Note what was *not* needed: normality, homoskedasticity, or even $E[\\varepsilon \\mid x] = 0$; only $E[x\\varepsilon] = 0$. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · The sandwich\n", "\n", "Derive $\\sqrt n(\\hat\\beta - \\beta) \\xrightarrow{d} N(0, Q^{-1}\\Omega Q^{-1})$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Scale Step 1 of Exercise 6 by $\\sqrt n$: $\\sqrt n(\\hat\\beta - \\beta) = \\hat Q^{-1}\\cdot\\frac{1}{\\sqrt n}\\sum_ix_i\\varepsilon_i$, with $\\hat Q = \\frac1n\\sum x_ix_i^{\\top}$.\n", "\n", "**Step 2.** The vectors $x_i\\varepsilon_i$ are iid with mean 0 and covariance $\\Omega = E[\\varepsilon^2xx^{\\top}]$. The multivariate CLT gives $\\frac{1}{\\sqrt n}\\sum x_i\\varepsilon_i \\xrightarrow{d} N(0, \\Omega)$.\n", "\n", "**Step 3.** $\\hat Q^{-1} \\xrightarrow{p} Q^{-1}$ (Exercise 6).\n", "\n", "**Step 4.** Slutsky: $\\sqrt n(\\hat\\beta - \\beta) \\xrightarrow{d} Q^{-1}Z$ with $Z \\sim N(0, \\Omega)$.\n", "\n", "**Step 5.** By Day 22, $\\operatorname{Var}(Q^{-1}Z) = Q^{-1}\\Omega Q^{-1}$ ($Q^{-1}$ is symmetric). So the limit is $N(0, Q^{-1}\\Omega Q^{-1})$. $\\blacksquare$\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · Homoskedasticity collapses the sandwich\n", "\n", "Show that if $E[\\varepsilon^2 \\mid x] = \\sigma^2$ then $Q^{-1}\\Omega Q^{-1} = \\sigma^2Q^{-1}$.\n", "\n", "
Show solution\n", "\n", "By the LIE, $\\Omega = E[E[\\varepsilon^2 \\mid x]xx^{\\top}] = \\sigma^2E[xx^{\\top}] = \\sigma^2Q$. So $Q^{-1}\\sigma^2QQ^{-1} = \\sigma^2Q^{-1}$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · The linear probability model is heteroskedastic\n", "\n", "For a binary outcome with $P(y = 1 \\mid x) = x^{\\top}\\beta$, show $\\operatorname{Var}(\\varepsilon \\mid x) = x^{\\top}\\beta(1 - x^{\\top}\\beta)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\varepsilon = y - x^{\\top}\\beta$ and $y \\mid x$ is Bernoulli($p(x)$), so $\\operatorname{Var}(\\varepsilon \\mid x) = \\operatorname{Var}(y \\mid x) = p(x)(1 - p(x))$, which varies with $x$.\n", "\n", "**Step 2.** Robust SEs are mandatory for linear probability models.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Robust SE of a difference in means\n", "\n", "In $y = \\beta_0 + \\beta_1T + \\varepsilon$, show that the HC0 variance of $\\hat\\beta_1$ equals $\\hat\\sigma_1^2/n_1 + \\hat\\sigma_0^2/n_0$ with plug-in (divide-by-$n_g$) group variances.\n", "\n", "
Show solution\n", "\n", "**Step 1.** By FWL (Day 28), the slope uses the demeaned treatment $\\tilde T_i = T_i - p$, where $p = n_1/n$: $\\hat\\beta_1 = \\sum\\tilde T_iy_i/\\sum\\tilde T_i^2$. Its HC0 variance is $\\sum\\tilde T_i^2e_i^2\\big/\\big(\\sum\\tilde T_i^2\\big)^2$.\n", "\n", "**Step 2.** Ingredients: $\\tilde T_i = 1 - p$ for treated and $-p$ for controls; $e_i = y_i - \\bar y_{g(i)}$ (deviation from own group mean); $\\sum\\tilde T_i^2 = n_1(1 - p)^2 + n_0p^2 = np(1 - p)$.\n", "\n", "**Step 3.** Numerator, writing $\\hat\\sigma_g^2 = \\frac{1}{n_g}\\sum_{i \\in g}e_i^2$: $(1 - p)^2n_1\\hat\\sigma_1^2 + p^2n_0\\hat\\sigma_0^2$.\n", "\n", "**Step 4.** Divide by $n^2p^2(1 - p)^2$ and use $n_1 = np$, $n_0 = n(1 - p)$:\n", "\n", "$$\n", "\\frac{(1 - p)^2\\,np\\,\\hat\\sigma_1^2}{n^2p^2(1 - p)^2} + \\frac{p^2\\,n(1 - p)\\,\\hat\\sigma_0^2}{n^2p^2(1 - p)^2} = \\frac{\\hat\\sigma_1^2}{np} + \\frac{\\hat\\sigma_0^2}{n(1 - p)} = \\frac{\\hat\\sigma_1^2}{n_1} + \\frac{\\hat\\sigma_0^2}{n_0}.\n", "$$\n", "\n", "**Answer:** the robust SE of a regression on a treatment dummy is Neyman's unequal-variance SE from Day 9.\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · Consumption, income and heteroskedasticity\n", "\n", "In household data, consumption varies much more among high-income than low-income households. Model this as $c_i = 1 + y_i + \\varepsilon_i$ with $\\operatorname{sd}(\\varepsilon_i \\mid y_i) = |y_i|$ (income $y$ standardised). Predict whether classical SEs for the marginal propensity to consume over- or under-state its uncertainty, and check by simulation (code below).\n", "\n", "
Show solution\n", "\n", "**Step 1.** The slope is driven by households far from mean income, where $(y_i - \\bar y)^2$ is large.\n", "\n", "**Step 2.** Here the error variance is largest exactly for those households.\n", "\n", "**Step 3.** The classical formula uses one average error variance for everyone, so it ignores that the influential observations are also the noisiest.\n", "\n", "**Answer:** classical SEs **understate** the uncertainty. In the simulation, classical 95% intervals cover only about 76% of the time, while HC1 robust intervals cover about 94.5%.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · The Moulton factor\n", "\n", "The regressor is constant within clusters of size $m = 50$ and the error's intra-cluster correlation is 0.1. By how much are classical SEs too small?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Moulton factor for the variance: $1 + (m - 1)\\rho_x\\rho_\\varepsilon = 1 + 49(1)(0.1) = 5.9$.\n", "\n", "**Step 2.** SE factor: $\\sqrt{5.9} \\approx 2.43$.\n", "\n", "**Answer:** classical SEs are about 2.4 times too small: a reported $t = 4$ might really be $t \\approx 1.6$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · A joint Wald test\n", "\n", "Two coefficients are estimated at $(0.3, 0.2)$ with robust covariance $\\operatorname{diag}(0.01, 0.01)$. Test that both are zero.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $W = \\hat\\beta^{\\top}V^{-1}\\hat\\beta$ with $V = \\operatorname{diag}(0.01, 0.01)$, so $V^{-1} = \\operatorname{diag}(100, 100)$.\n", "\n", "**Step 2.** $W = 100(0.3^2) + 100(0.2^2) = 9 + 4 = 13$.\n", "\n", "**Step 3.** Reference: $\\chi^2_2$, 5% critical value 5.99.\n", "\n", "**Answer:** $13 > 5.99$: reject that both coefficients are zero. In general $W = (R\\hat\\beta - r)^{\\top}(RVR^{\\top})^{-1}(R\\hat\\beta - r)$ with $q$ degrees of freedom.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · Few clusters\n", "\n", "A state-level policy is studied with 6 treated and 44 control states. Why can cluster-robust SEs mislead, and what do practitioners do?\n", "\n", "
Show solution\n", "\n", "**Step 1.** The cluster-robust \"meat\" is a sum over clusters; it is reliable only with many clusters, especially many *treated* ones.\n", "\n", "**Step 2.** With 6 treated states, those few clusters' residuals dominate, so the variance estimate is noisy and typically biased downward.\n", "\n", "**Step 3.** Remedies: the wild cluster bootstrap (with restricted residuals), randomisation inference across states, or synthetic-control methods (Phase 6).\n", "\n", "**Answer:** standard cluster-robust SEs are too small with few treated clusters; use one of the remedies.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · Exogeneity versus orthogonality\n", "\n", "The asymptotic results need only $E[x\\varepsilon] = 0$. Why does causal interpretation usually require the stronger $E[\\varepsilon \\mid x] = 0$?\n", "\n", "
Show solution\n", "\n", "**Step 1.** For the *best linear predictor* of $y$ given $x$, the error is defined so that $E[x\\varepsilon] = 0$; OLS always estimates those projection coefficients consistently.\n", "\n", "**Step 2.** A causal reading needs a *structural* error: all other causes of $y$.\n", "\n", "**Step 3.** $E[\\varepsilon \\mid x] = 0$ for that error says the regressors are unrelated to every omitted cause: the identifying assumption.\n", "\n", "**Answer:** asymptotic theory delivers precision for whatever $\\hat\\beta$ estimates; whether that is a causal effect depends on identification.\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", "rng = np.random.default_rng(29)\n", "n, reps, beta = 200, 4000, 1.0\n", "cover_cl = cover_hc = 0\n", "for _ in range(reps):\n", " x = rng.normal(size=n); y = 1 + beta * x + np.abs(x) * rng.normal(size=n)\n", " X = np.column_stack([np.ones(n), x]); XtXi = np.linalg.inv(X.T @ X)\n", " b = XtXi @ X.T @ y; e = y - X @ b\n", " v_cl = (e @ e / (n - 2)) * XtXi[1, 1]\n", " v_hc = (XtXi @ (X.T * e**2) @ X @ XtXi)[1, 1] * n / (n - 2) # HC1\n", " cover_cl += abs(b[1] - beta) <= 1.96 * np.sqrt(v_cl)\n", " cover_hc += abs(b[1] - beta) <= 1.96 * np.sqrt(v_hc)\n", "print(f\"95% CI coverage: classical {cover_cl / reps:.3f}, robust HC1 {cover_hc / reps:.3f}\")\n", "\n", "x = np.array([1, 2, 3.]); e = np.array([0.5, -1, 0.5])\n", "print(\"classical:\", round((e @ e / 2) / (x @ x), 4), \" HC0:\", round((x**2 * e**2).sum() / (x @ x)**2, 4))\n", "print(\"Wald:\", 0.3**2 / 0.01 + 0.2**2 / 0.01, \" Moulton SE factor:\", round(np.sqrt(1 + 49 * 0.1), 2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- OLS is consistent and asymptotically normal under weak conditions; its variance is a sandwich.\n", "- Use heteroskedasticity-robust SEs by default, and cluster at the level of treatment assignment.\n", "- Asymptotic theory tells you how precise $\\hat\\beta$ is, not whether it is causal." ] } ] }