{
 "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 64: 2SLS in Matrix Form\n",
    "\n",
    "*Beginner to Ninja: Math for Causal Inference*\n",
    "\n",
    "Two-stage least squares generalises IV to many instruments and controls. Project the regressors on the instruments, regress the outcome on the projections, and get the variance right. Includes the GMM interpretation and the classic mistakes (manual second-stage standard errors, the forbidden regression).\n",
    "\n",
    "**Learning outcomes**\n",
    "\n",
    "- Derive the 2SLS estimator with the projection matrix and show it equals OLS on first-stage fitted values.\n",
    "- State the order and rank conditions for identification.\n",
    "- Compute 2SLS by hand in small examples and see it as a weighted average of single-instrument estimates.\n",
    "- Explain why manual second-stage standard errors are wrong and why nonlinear first stages must not be plugged in."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "### Set-up\n",
    "\n",
    "$\\mathbf y = \\mathbf X\\boldsymbol\\beta + \\boldsymbol\\varepsilon$, where $\\mathbf X$ ($n \\times k$) contains the endogenous regressors and the exogenous controls (including the constant). $\\mathbf Z$ ($n \\times l$) contains the excluded instruments **and the same exogenous controls**. Assumptions: $E[\\mathbf z_i\\varepsilon_i] = \\mathbf 0$ (instrument exogeneity) and $\\operatorname{rank}E[\\mathbf z_i\\mathbf x_i^\\top] = k$ (rank condition).\n",
    "\n",
    "- **Order condition:** $l \\ge k$: at least as many instruments as endogenous regressors. $l = k$: just identified; $l > k$: over-identified.\n",
    "\n",
    "### The estimator\n",
    "\n",
    "With the projection matrix $\\mathbf P_Z = \\mathbf Z(\\mathbf Z^\\top\\mathbf Z)^{-1}\\mathbf Z^\\top$,\n",
    "\n",
    "$$\n",
    "\\hat{\\boldsymbol\\beta}_{2SLS} = (\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\mathbf y.\n",
    "$$\n",
    "\n",
    "**Two stages.** (1) Regress each column of $\\mathbf X$ on $\\mathbf Z$: $\\hat{\\mathbf X} = \\mathbf P_Z\\mathbf X$ (exogenous columns are reproduced exactly). (2) Regress $\\mathbf y$ on $\\hat{\\mathbf X}$. Because $\\mathbf P_Z$ is symmetric and idempotent, $\\hat{\\mathbf X}^\\top\\hat{\\mathbf X} = \\mathbf X^\\top\\mathbf P_Z\\mathbf X$, so the second-stage OLS coefficients are exactly 2SLS.\n",
    "\n",
    "**Just identified:** $\\hat{\\boldsymbol\\beta} = (\\mathbf Z^\\top\\mathbf X)^{-1}\\mathbf Z^\\top\\mathbf y$, the IV estimator.\n",
    "\n",
    "### Variance\n",
    "\n",
    "Under homoskedasticity, $\\widehat{\\operatorname{Var}}(\\hat{\\boldsymbol\\beta}) = \\hat\\sigma^2(\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}$ with $\\hat\\sigma^2 = \\frac1n\\sum(y_i - \\mathbf x_i^\\top\\hat{\\boldsymbol\\beta})^2$, residuals computed with the **actual** $\\mathbf X$. The robust (sandwich) version is\n",
    "\n",
    "$$\n",
    "(\\hat{\\mathbf X}^\\top\\hat{\\mathbf X})^{-1}\\Big(\\sum_i\\hat e_i^2\\,\\hat{\\mathbf x}_i\\hat{\\mathbf x}_i^\\top\\Big)(\\hat{\\mathbf X}^\\top\\hat{\\mathbf X})^{-1}.\n",
    "$$\n",
    "\n",
    "Running the second stage by hand in software uses residuals $y_i - \\hat{\\mathbf x}_i^\\top\\hat{\\boldsymbol\\beta}$, which gives the wrong standard errors.\n",
    "\n",
    "### GMM view\n",
    "\n",
    "2SLS minimises $(\\mathbf Z^\\top(\\mathbf y - \\mathbf X\\mathbf b))^\\top(\\mathbf Z^\\top\\mathbf Z)^{-1}(\\mathbf Z^\\top(\\mathbf y - \\mathbf X\\mathbf b))$: GMM (Day 30) with weight matrix $(\\mathbf Z^\\top\\mathbf Z)^{-1}$, which is efficient under homoskedasticity. With heteroskedasticity, two-step GMM can be more efficient.\n",
    "\n",
    "### Many instruments: a weighted average\n",
    "\n",
    "With one endogenous regressor and orthogonal instruments, 2SLS is a weighted average of the single-instrument IV estimates, with weights proportional to the squared first-stage strength of each instrument (Exercise 11).\n",
    "\n",
    "### Two mistakes to avoid\n",
    "\n",
    "- **Controls in one stage only.** Every exogenous control must appear in both stages; otherwise the first-stage residual is correlated with the omitted control.\n",
    "- **The forbidden regression.** With a binary treatment, plugging fitted probabilities from a probit into the second stage is not 2SLS and is inconsistent unless the probit is exactly right. Use the fitted probabilities as an *instrument* instead, or just use the linear first stage.\n",
    "\n",
    "## Part 1 · Warm-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1 · Just-identified by hand\n",
    "\n",
    "Demeaned data (no constant): $z = (1, -1, 1, -1)$, $x = (2, 0, 3, -1)$, $y = (5, 1, 7, -1)$. Compute $\\hat\\beta_{IV}$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\mathbf z^\\top\\mathbf y = 5 - 1 + 7 + 1 = 12$.\n",
    "\n",
    "**Step 2.** $\\mathbf z^\\top\\mathbf x = 2 - 0 + 3 + 1 = 6$.\n",
    "\n",
    "**Step 3.** $\\hat\\beta = 12/6 = 2$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 · Properties of the projection\n",
    "\n",
    "Show that $\\mathbf P_Z$ is symmetric and idempotent.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\mathbf P_Z^\\top = \\mathbf Z(\\mathbf Z^\\top\\mathbf Z)^{-\\top}\\mathbf Z^\\top = \\mathbf P_Z$, since $\\mathbf Z^\\top\\mathbf Z$ is symmetric.\n",
    "\n",
    "**Step 2.** $\\mathbf P_Z\\mathbf P_Z = \\mathbf Z(\\mathbf Z^\\top\\mathbf Z)^{-1}\\mathbf Z^\\top\\mathbf Z(\\mathbf Z^\\top\\mathbf Z)^{-1}\\mathbf Z^\\top = \\mathbf P_Z$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3 · First-stage fitted values\n",
    "\n",
    "In Exercise 1, compute the first-stage coefficient and fitted values $\\hat x$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\hat\\pi = \\mathbf z^\\top\\mathbf x/\\mathbf z^\\top\\mathbf z = 6/4 = 1.5$.\n",
    "\n",
    "**Step 2.** $\\hat x = 1.5z = (1.5, -1.5, 1.5, -1.5)$.\n",
    "\n",
    "**Step 3.** Second stage: $\\hat x^\\top y/\\hat x^\\top\\hat x = 1.5(12)/(2.25 \\times 4) = 18/9 = 2$, matching Exercise 1.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4 · Order condition\n",
    "\n",
    "A wage equation has two endogenous regressors (schooling and experience) and one instrument (quarter of birth). Is it identified? What if you add age as an instrument for experience?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $k_{endog} = 2$, excluded instruments 1: order condition fails; not identified.\n",
    "\n",
    "**Step 2.** With age: 2 instruments for 2 endogenous regressors: just identified (if the rank condition holds: age must predict experience separately from quarter of birth).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5 · Controls in both stages\n",
    "\n",
    "Why must region dummies that appear in the wage equation also appear in the first stage?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** If region is omitted from the first stage, the first-stage residual $X - \\hat X$ contains the part of $X$ explained by region.\n",
    "\n",
    "**Step 2.** Region is in the second stage, so that residual is correlated with a second-stage regressor; the estimator no longer equals $(\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\mathbf y$ with a valid $\\mathbf Z$ and becomes inconsistent in general.\n",
    "\n",
    "**Answer:** exogenous controls are their own instruments; include them in $\\mathbf Z$.\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 2 · Proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6 · Two stages give 2SLS\n",
    "\n",
    "Prove that OLS of $\\mathbf y$ on $\\hat{\\mathbf X} = \\mathbf P_Z\\mathbf X$ equals $(\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\mathbf y$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** OLS: $(\\hat{\\mathbf X}^\\top\\hat{\\mathbf X})^{-1}\\hat{\\mathbf X}^\\top\\mathbf y$.\n",
    "\n",
    "**Step 2.** $\\hat{\\mathbf X}^\\top\\hat{\\mathbf X} = \\mathbf X^\\top\\mathbf P_Z^\\top\\mathbf P_Z\\mathbf X = \\mathbf X^\\top\\mathbf P_Z\\mathbf X$ (Exercise 2).\n",
    "\n",
    "**Step 3.** $\\hat{\\mathbf X}^\\top\\mathbf y = \\mathbf X^\\top\\mathbf P_Z\\mathbf y$. $\\blacksquare$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7 · Just identification\n",
    "\n",
    "If $l = k$ and $\\mathbf Z^\\top\\mathbf X$ is invertible, prove $\\hat{\\boldsymbol\\beta}_{2SLS} = (\\mathbf Z^\\top\\mathbf X)^{-1}\\mathbf Z^\\top\\mathbf y$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\mathbf X^\\top\\mathbf P_Z\\mathbf X = (\\mathbf X^\\top\\mathbf Z)(\\mathbf Z^\\top\\mathbf Z)^{-1}(\\mathbf Z^\\top\\mathbf X)$, a product of three square invertible matrices.\n",
    "\n",
    "**Step 2.** Its inverse is $(\\mathbf Z^\\top\\mathbf X)^{-1}(\\mathbf Z^\\top\\mathbf Z)(\\mathbf X^\\top\\mathbf Z)^{-1}$.\n",
    "\n",
    "**Step 3.** Multiply by $\\mathbf X^\\top\\mathbf P_Z\\mathbf y = (\\mathbf X^\\top\\mathbf Z)(\\mathbf Z^\\top\\mathbf Z)^{-1}\\mathbf Z^\\top\\mathbf y$: everything cancels except $(\\mathbf Z^\\top\\mathbf X)^{-1}\\mathbf Z^\\top\\mathbf y$. $\\blacksquare$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8 · Consistency\n",
    "\n",
    "Assume $\\frac1n\\mathbf Z^\\top\\mathbf Z \\to \\mathbf Q_{zz}$ (invertible), $\\frac1n\\mathbf Z^\\top\\mathbf X \\to \\mathbf Q_{zx}$ (rank $k$) and $\\frac1n\\mathbf Z^\\top\\boldsymbol\\varepsilon \\to \\mathbf 0$. Prove $\\hat{\\boldsymbol\\beta}_{2SLS} \\to \\boldsymbol\\beta$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Substitute $\\mathbf y = \\mathbf X\\boldsymbol\\beta + \\boldsymbol\\varepsilon$: $\\hat{\\boldsymbol\\beta} = \\boldsymbol\\beta + (\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\boldsymbol\\varepsilon$.\n",
    "\n",
    "**Step 2.** Scale by $1/n$: $\\frac1n\\mathbf X^\\top\\mathbf P_Z\\mathbf X \\to \\mathbf Q_{zx}^\\top\\mathbf Q_{zz}^{-1}\\mathbf Q_{zx}$, invertible by the rank condition.\n",
    "\n",
    "**Step 3.** $\\frac1n\\mathbf X^\\top\\mathbf P_Z\\boldsymbol\\varepsilon = (\\frac1n\\mathbf X^\\top\\mathbf Z)(\\frac1n\\mathbf Z^\\top\\mathbf Z)^{-1}(\\frac1n\\mathbf Z^\\top\\boldsymbol\\varepsilon) \\to \\mathbf Q_{zx}^\\top\\mathbf Q_{zz}^{-1}\\mathbf 0 = \\mathbf 0$.\n",
    "\n",
    "**Step 4.** By Slutsky, the second term vanishes. $\\blacksquare$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 9 · Why manual second-stage SEs are wrong\n",
    "\n",
    "Show that the manual residual $y_i - \\hat x_i\\hat\\beta$ equals $\\varepsilon_i + (x_i - \\hat x_i)\\beta$ asymptotically (scalar case), so $\\hat\\sigma^2_{manual} \\to \\sigma^2 + \\beta^2\\operatorname{Var}(v) + 2\\beta\\operatorname{Cov}(v, \\varepsilon)$, where $v = x - \\hat x$ is the first-stage error.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $y_i = x_i\\beta + \\varepsilon_i = \\hat x_i\\beta + (x_i - \\hat x_i)\\beta + \\varepsilon_i$.\n",
    "\n",
    "**Step 2.** With $\\hat\\beta \\to \\beta$, the manual residual is $\\varepsilon_i + v_i\\beta$.\n",
    "\n",
    "**Step 3.** Its mean square converges to $\\sigma^2 + \\beta^2\\operatorname{Var}(v) + 2\\beta\\operatorname{Cov}(v, \\varepsilon)$.\n",
    "\n",
    "**Answer:** unless $\\beta = 0$, this differs from $\\sigma^2$; the standard errors are wrong (usually too large, but too small when $\\beta\\operatorname{Cov}(v, \\varepsilon) < 0$ is large). $\\blacksquare$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 10 · 2SLS is GMM\n",
    "\n",
    "Show that minimising $J(\\mathbf b) = (\\mathbf y - \\mathbf X\\mathbf b)^\\top\\mathbf Z(\\mathbf Z^\\top\\mathbf Z)^{-1}\\mathbf Z^\\top(\\mathbf y - \\mathbf X\\mathbf b)$ gives the 2SLS estimator.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $J(\\mathbf b) = (\\mathbf y - \\mathbf X\\mathbf b)^\\top\\mathbf P_Z(\\mathbf y - \\mathbf X\\mathbf b)$.\n",
    "\n",
    "**Step 2.** Gradient: $-2\\mathbf X^\\top\\mathbf P_Z(\\mathbf y - \\mathbf X\\mathbf b) = \\mathbf 0$.\n",
    "\n",
    "**Step 3.** Solve: $\\mathbf b = (\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\mathbf y$. $\\blacksquare$\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 3 · Applied"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 11 · Two instruments for schooling\n",
    "\n",
    "With demeaned data and two orthogonal instruments (a quarter-of-birth indicator and a compulsory-schooling-law indicator), $\\mathbf Z^\\top\\mathbf Z = \\operatorname{diag}(100, 100)$, $\\mathbf Z^\\top\\mathbf x = (20, 10)^\\top$ and $\\mathbf Z^\\top\\mathbf y = (3, 1)^\\top$ (log wages). Compute the two single-instrument IV estimates and the 2SLS estimate, and show 2SLS is their weighted average.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Single instruments: $3/20 = 0.15$ and $1/10 = 0.10$.\n",
    "\n",
    "**Step 2.** $\\mathbf x^\\top\\mathbf P_Z\\mathbf x = (20^2 + 10^2)/100 = 5$.\n",
    "\n",
    "**Step 3.** $\\mathbf x^\\top\\mathbf P_Z\\mathbf y = (20 \\times 3 + 10 \\times 1)/100 = 0.7$.\n",
    "\n",
    "**Step 4.** $\\hat\\beta_{2SLS} = 0.7/5 = 0.14$.\n",
    "\n",
    "**Step 5.** Weights proportional to $(\\mathbf z_j^\\top\\mathbf x)^2$: $400/500 = 0.8$ and $100/500 = 0.2$; $0.8(0.15) + 0.2(0.10) = 0.14$.\n",
    "\n",
    "**Answer:** 14% per year; the stronger instrument gets four times the weight.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 12 · A first stage with controls\n",
    "\n",
    "A study of the effect of class size on test scores instruments class size with a rule that splits classes above 40 pupils (Maimonides' rule), controlling for enrolment. Write the two stages and say what identifies the effect.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** First stage: class size $= \\pi_0 + \\pi_1 f(\\text{enrolment}) + \\pi_2\\,\\text{enrolment} + v$, where $f$ is the rule's predicted class size.\n",
    "\n",
    "**Step 2.** Second stage: score $= \\beta_0 + \\beta_1\\widehat{\\text{class size}} + \\beta_2\\,\\text{enrolment} + \\varepsilon$.\n",
    "\n",
    "**Step 3.** Enrolment is in both stages; identification comes from the rule's discontinuous jumps in predicted class size, not from enrolment's smooth effect.\n",
    "\n",
    "**Answer:** a 2SLS with an RD flavour (Day 76).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 13 · Manual standard errors in a simulation\n",
    "\n",
    "In the simulation below (returns to schooling, true effect 0.08), compare the correct 2SLS standard error with the one from running the second stage by OLS.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** The correct SE uses residuals $y - X\\hat\\beta$.\n",
    "\n",
    "**Step 2.** The manual SE uses $y - \\hat X\\hat\\beta$, which adds $\\beta v$ to every residual (Exercise 9).\n",
    "\n",
    "**Step 3.** Here the first-stage error has a large variance relative to $\\sigma^2$, so the manual SE is wrong by a noticeable factor (the code prints both, and the Monte Carlo standard deviation for reference).\n",
    "\n",
    "**Answer:** always use the IV routine's standard errors (or compute residuals with the original regressors).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 14 · The forbidden regression\n",
    "\n",
    "To estimate the effect of hospital admission $D$ (binary) on costs with distance as instrument, an analyst fits a probit of $D$ on distance and covariates, then regresses costs on the fitted probabilities. Why is this risky, and what should be done instead?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** 2SLS's consistency relies on the first-stage residual being orthogonal to all second-stage regressors by construction (linear projection).\n",
    "\n",
    "**Step 2.** A probit's residual has no such property; if the probit is misspecified, the second stage inherits a correlation between the fitted value and the error.\n",
    "\n",
    "**Step 3.** Safe alternatives: linear first stage (2SLS), or use the probit fitted probability as an *instrument* for $D$ in a linear IV.\n",
    "\n",
    "**Answer:** the plug-in version is consistent only if the probit is exactly right; using the probit as an instrument stays consistent regardless.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 15 · Heteroskedasticity-robust standard errors\n",
    "\n",
    "In an IV regression of firm productivity on management training (instrumented by a subsidy lottery), residual variance is larger for small firms. Write the robust variance estimator and explain why it is preferred.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\widehat{\\operatorname{Var}} = (\\hat{\\mathbf X}^\\top\\hat{\\mathbf X})^{-1}\\big(\\sum_i\\hat e_i^2\\hat{\\mathbf x}_i\\hat{\\mathbf x}_i^\\top\\big)(\\hat{\\mathbf X}^\\top\\hat{\\mathbf X})^{-1}$ with $\\hat e_i = y_i - \\mathbf x_i^\\top\\hat{\\boldsymbol\\beta}$.\n",
    "\n",
    "**Step 2.** The homoskedastic formula assumes $E[\\varepsilon_i^2 \\mid \\mathbf z_i] = \\sigma^2$, false here.\n",
    "\n",
    "**Answer:** the sandwich is consistent either way (Day 29); with clustered assignment (lotteries by region), cluster the middle term.\n",
    "\n",
    "</details>\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",
    "# Exercises 1 and 3\n",
    "z = np.array([1, -1, 1, -1.]); x = np.array([2, 0, 3, -1.]); y = np.array([5, 1, 7, -1.])\n",
    "print(\"IV:\", z @ y / (z @ x), \" first stage:\", z @ x / (z @ z))\n",
    "\n",
    "# Exercise 11 from cross-products\n",
    "ZZ = np.diag([100., 100.]); Zx = np.array([20., 10.]); Zy = np.array([3., 1.])\n",
    "print(\"single IVs:\", Zy / Zx, \" 2SLS:\", round((Zx @ np.linalg.solve(ZZ, Zy)) / (Zx @ np.linalg.solve(ZZ, Zx)), 4))\n",
    "\n",
    "# 2SLS by matrices, correct versus manual standard errors, with a Monte Carlo check (Exercise 13)\n",
    "rng = np.random.default_rng(64)\n",
    "def draw(n=2000):\n",
    "    A = rng.normal(size=n); Z = rng.normal(size=(n, 2))\n",
    "    S = 12 + Z @ np.array([.5, .3]) + 1.5 * A + rng.normal(0, 2, n)\n",
    "    lw = 1 + .08 * S + .1 * A + rng.normal(0, .3, n)\n",
    "    return lw, np.column_stack([np.ones(n), S]), np.column_stack([np.ones(n), Z])\n",
    "def tsls(y, X, Z):\n",
    "    Xh = Z @ np.linalg.lstsq(Z, X, rcond=None)[0]\n",
    "    b = np.linalg.solve(Xh.T @ Xh, Xh.T @ y)\n",
    "    inv = np.linalg.inv(Xh.T @ Xh)\n",
    "    se_ok = np.sqrt((y - X @ b) @ (y - X @ b) / len(y) * inv[1, 1])\n",
    "    se_manual = np.sqrt((y - Xh @ b) @ (y - Xh @ b) / len(y) * inv[1, 1])\n",
    "    return b[1], se_ok, se_manual\n",
    "y, X, Z = draw()\n",
    "b, s1, s2 = tsls(y, X, Z)\n",
    "print(f\"2SLS {b:.4f}  correct SE {s1:.4f}  manual SE {s2:.4f}\")\n",
    "reps = [tsls(*draw())[0] for _ in range(500)]\n",
    "print(f\"Monte Carlo SD of 2SLS: {np.std(reps):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key takeaways\n",
    "\n",
    "- $\\hat{\\boldsymbol\\beta}_{2SLS} = (\\mathbf X^\\top\\mathbf P_Z\\mathbf X)^{-1}\\mathbf X^\\top\\mathbf P_Z\\mathbf y$: OLS on first-stage fitted values, with exogenous controls in both stages.\n",
    "- Identification needs as many instruments as endogenous regressors (order) and a full-rank first stage (rank).\n",
    "- Compute residuals with the actual regressors; never plug nonlinear first-stage predictions into the second stage."
   ]
  }
 ]
}