{ "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 30: Phase 2 Capstone: GMM & Numerical Optimisation\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "One framework for OLS, IV, 2SLS and MLE: choose moment conditions, minimise a quadratic form, and solve it numerically.\n", "\n", "**Learning outcomes**\n", "\n", "- Write OLS, IV and maximum likelihood as method-of-moments estimators.\n", "- Set up the GMM objective, derive its first-order condition and the optimal weighting matrix.\n", "- Show that 2SLS is GMM with a particular weighting matrix and interpret the J test of over-identifying restrictions.\n", "- Run Newton–Raphson by hand for a logit and use numerical derivatives safely." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Method of moments\n", "\n", "Many estimators come from **population moment conditions** $E[g(w_i, \\theta_0)] = 0$, where $w_i$ is the data and $\\theta$ has $K$ elements. The **method of moments** replaces the expectation by a sample average and solves\n", "\n", "$$\n", "\\bar g(\\theta) = \\frac1n\\sum_{i=1}^n g(w_i, \\theta) = 0.\n", "$$\n", "\n", "| Estimator | Moment condition |\n", "|---|---|\n", "| Mean | $E[y - \\mu] = 0$ |\n", "| OLS | $E[x(y - x^{\\top}\\beta)] = 0$ |\n", "| IV | $E[z(y - x^{\\top}\\beta)] = 0$ |\n", "| MLE | $E[\\partial\\ln f(w; \\theta)/\\partial\\theta] = 0$ (the score) |\n", "\n", "### Generalised method of moments (GMM)\n", "\n", "With $L > K$ moment conditions (over-identification), they cannot all be set to zero in the sample. GMM minimises a weighted distance from zero:\n", "\n", "$$\n", "\\hat\\theta = \\arg\\min_\\theta\\; Q(\\theta) = \\bar g(\\theta)^{\\top}W\\bar g(\\theta).\n", "$$\n", "\n", "With $G = E[\\partial g/\\partial\\theta^{\\top}]$ and $S = E[g\\,g^{\\top}]$:\n", "\n", "- FOC: $\\hat G^{\\top}W\\bar g(\\hat\\theta) = 0$.\n", "- Asymptotic variance: $(G^{\\top}WG)^{-1}G^{\\top}WSWG(G^{\\top}WG)^{-1}$ (another sandwich).\n", "- **Optimal weighting** $W = S^{-1}$ gives the smallest variance, $(G^{\\top}S^{-1}G)^{-1}$: trust precise moments more (compare Day 26's inverse-variance weights).\n", "- **J test:** $J = n\\,\\bar g(\\hat\\theta)^{\\top}\\hat S^{-1}\\bar g(\\hat\\theta) \\xrightarrow{d} \\chi^2_{L-K}$ under correct specification.\n", "\n", "**2SLS** is GMM with moments $z_i(y_i - x_i^{\\top}\\beta)$ and $W = (Z^{\\top}Z/n)^{-1}$: optimal under homoskedasticity.\n", "\n", "### Solving it numerically\n", "\n", "Most GMM and ML problems have no closed form. **Newton–Raphson** (Day 25) iterates $\\theta \\leftarrow \\theta - H^{-1}\\nabla Q$; **Gauss–Newton** approximates the Hessian with first derivatives of $\\bar g$. When derivatives are not available analytically, the **central difference** $[f(x + h) - f(x - h)]/2h$ has error of order $h^2$, but $h$ too small amplifies rounding error; $h \\approx 10^{-5}$ to $10^{-6}$ (relative to $x$) is a sensible default.\n", "\n", "### Phase 2 in one table\n", "\n", "| Tool | Where it returns |\n", "|---|---|\n", "| Dot products, norms (16) | covariances, matching distances |\n", "| $X^{\\top}X$, $X^{\\top}y$ (17) | every regression |\n", "| Rank (18) | identification, collinearity, fixed effects |\n", "| Inverses, conditioning (19) | IV, weak instruments |\n", "| Projections (20) | residuals, $R^2$, leverage |\n", "| Eigenvalues (21) | PCA, condition numbers, stability |\n", "| Covariance matrices (22) | SEs of combinations, simulation |\n", "| Derivatives, Taylor (23) | marginal effects, delta method |\n", "| Matrix calculus (24) | first-order conditions |\n", "| Optimisation (25–26) | MLE, weighting, design |\n", "| OLS, FWL, sandwich (27–29) | regression, fixed effects, robust inference |\n", "| GMM (30) | IV, 2SLS, efficient estimation |\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · Exponential by moments\n", "\n", "For $X \\sim \\text{Exponential}(\\lambda)$, $E[X] = 1/\\lambda$. Find the method-of-moments estimator.\n", "\n", "
Show solution\n", "\n", "Set $\\bar x - 1/\\lambda = 0$: $\\hat\\lambda = 1/\\bar x$, the same as the MLE (Day 11).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · Uniform by moments\n", "\n", "For $X \\sim U(0, \\theta)$, find the method-of-moments estimator. Is it always sensible?\n", "\n", "
Show solution\n", "\n", "**Step 1.** $E[X] = \\theta/2 \\Rightarrow \\hat\\theta = 2\\bar x$.\n", "\n", "**Step 2.** Not always: with data $(1, 1, 7)$, $\\hat\\theta = 6$ though an observation equals 7, impossible under the model.\n", "\n", "**Step 3.** (The MLE is $\\max x_i$.) Moment estimators are consistent but not always efficient or coherent.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · OLS as a moment estimator\n", "\n", "Show that solving the sample analogue of $E[x(y - x^{\\top}\\beta)] = 0$ gives OLS.\n", "\n", "
Show solution\n", "\n", "$\\frac1n\\sum x_i(y_i - x_i^{\\top}\\beta) = 0 \\iff X^{\\top}y = X^{\\top}X\\beta$, the normal equations (Day 24).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · IV as a moment estimator\n", "\n", "With as many instruments as regressors, solve $\\frac1nZ^{\\top}(y - X\\beta) = 0$.\n", "\n", "
Show solution\n", "\n", "$Z^{\\top}y = Z^{\\top}X\\beta \\Rightarrow \\hat\\beta_{IV} = (Z^{\\top}X)^{-1}Z^{\\top}y$ (Day 19, Exercise 15), which requires $Z^{\\top}X$ invertible: instruments must be relevant.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · Counting\n", "\n", "A model has $K = 2$ parameters. How many moment conditions are needed for identification, and what does the J test need?\n", "\n", "
Show solution\n", "\n", "At least 2 (order condition; the rank condition requires $G$ to have full column rank). The J test needs over-identification: with $L$ moments it has $L - 2$ degrees of freedom, so $L \\ge 3$.\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · The GMM first-order condition\n", "\n", "Derive $\\hat G^{\\top}W\\bar g(\\hat\\theta) = 0$ from $Q(\\theta) = \\bar g^{\\top}W\\bar g$ with symmetric $W$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** By the chain rule and Day 24's quadratic-form rule, $\\nabla Q = 2\\big(\\partial\\bar g/\\partial\\theta^{\\top}\\big)^{\\top}W\\bar g = 2\\hat G^{\\top}W\\bar g$.\n", "\n", "**Step 2.** Set it to zero.\n", "\n", "**Step 3.** With $L = K$ and $\\hat G$ invertible, this forces $\\bar g = 0$: just-identified GMM is method of moments, whatever $W$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · 2SLS is GMM\n", "\n", "With $\\bar g(\\beta) = \\frac1nZ^{\\top}(y - X\\beta)$ and $W = (Z^{\\top}Z/n)^{-1}$, show the GMM estimator is $\\hat\\beta_{2SLS} = (X^{\\top}P_ZX)^{-1}X^{\\top}P_Zy$ with $P_Z = Z(Z^{\\top}Z)^{-1}Z^{\\top}$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\hat G = -Z^{\\top}X/n$.\n", "\n", "**Step 2.** The FOC $X^{\\top}Z(Z^{\\top}Z)^{-1}Z^{\\top}(y - X\\beta) = 0$ gives $X^{\\top}P_ZX\\beta = X^{\\top}P_Zy$.\n", "\n", "**Step 3.** Since $P_ZX = \\hat X$ (first-stage fitted values) and $P_Z$ is idempotent, this is OLS of $y$ on $\\hat X$: the \"two stages\".\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · The optimal weighting matrix\n", "\n", "Show that with $W = S^{-1}$ the sandwich simplifies to $(G^{\\top}S^{-1}G)^{-1}$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Substitute $W = S^{-1}$: the middle becomes $G^{\\top}S^{-1}SS^{-1}G = G^{\\top}S^{-1}G$, so the variance is $(G^{\\top}S^{-1}G)^{-1}(G^{\\top}S^{-1}G)(G^{\\top}S^{-1}G)^{-1} = (G^{\\top}S^{-1}G)^{-1}$.\n", "\n", "**Step 2.** (That no other $W$ does better is a Gauss–Markov-style argument.)\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · Why J has L − K degrees of freedom\n", "\n", "Explain intuitively why the J statistic is $\\chi^2_{L-K}$ rather than $\\chi^2_L$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\sqrt n\\,\\hat S^{-1/2}\\bar g(\\theta_0)$ is approximately a vector of $L$ independent standard normals, so its squared length would be $\\chi^2_L$.\n", "\n", "**Step 2.** But we evaluate at $\\hat\\theta$, chosen to make $\\bar g$ small, which uses up $K$ directions (a projection, Day 20).\n", "\n", "**Step 3.** The remaining $L - K$ directions are free, and their squared length is $\\chi^2_{L-K}$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Two moments for doctor visits\n", "\n", "Annual doctor visits are modelled as Poisson($\\lambda$), which implies $E[X] = \\lambda$ and $E[X^2] = \\lambda + \\lambda^2$. A health survey finds $\\bar x = 3$ and $\\overline{x^2} = 13$ (so the sample variance is 4). Write the GMM objective $Q(\\lambda)$ with $W = I$ and find its minimiser numerically.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Sample moment conditions: $\\bar g(\\lambda) = (3 - \\lambda,\\; 13 - \\lambda - \\lambda^2)$.\n", "\n", "**Step 2.** With $W = I$: $Q(\\lambda) = (3 - \\lambda)^2 + (13 - \\lambda - \\lambda^2)^2$.\n", "\n", "**Step 3.** At $\\lambda = 3$ the first moment fits exactly but the second does not: $Q(3) = 0 + (13 - 12)^2 = 1$. The second moment wants a larger $\\lambda$.\n", "\n", "**Step 4.** Minimising numerically (code below): $\\hat\\lambda \\approx 3.14$.\n", "\n", "**Step 5.** With the optimal $W = \\hat S^{-1}$, the noisy second moment gets little weight and $\\hat\\lambda$ returns to $\\bar x = 3.00$.\n", "\n", "**Answer:** about 3.14 with $W = I$ and 3.00 with optimal weights. Variance (4) above the mean (3) is overdispersion (Day 4); in a large sample the J test would flag it.\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · Newton–Raphson for training take-up\n", "\n", "30% of eligible workers enrol in a training programme. Fit an intercept-only logit for enrolment. Starting from $\\beta_0 = 0$, do two Newton steps and compare with the exact MLE $\\ln(0.3/0.7)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Per worker, the score is $\\bar y - p$ and the Hessian is $-p(1 - p)$, where $p = \\Lambda(\\beta)$. Newton step: $\\beta \\leftarrow \\beta + (\\bar y - p)/(p(1 - p))$.\n", "\n", "**Step 2.** From $\\beta = 0$: $p = 0.5$, so $\\beta_1 = 0 + (0.3 - 0.5)/0.25 = -0.8$.\n", "\n", "**Step 3.** From $\\beta = -0.8$: $p = 1/(1 + e^{0.8}) \\approx 0.310$ and $p(1 - p) \\approx 0.2139$, so $\\beta_2 = -0.8 + (0.3 - 0.310)/0.2139 \\approx -0.847$.\n", "\n", "**Step 4.** Exact: $\\ln(0.3/0.7) = \\ln(0.4286) \\approx -0.8473$.\n", "\n", "**Answer:** two steps reach three decimals: Newton's quadratic convergence.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Two instruments disagree\n", "\n", "Two valid-looking instruments give Wald estimates 2.0 (SE 0.3) and 3.0 (SE 0.3). Would an over-identification test reject? What are the possible explanations?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Difference between the two estimates: $3.0 - 2.0 = 1.0$.\n", "\n", "**Step 2.** SE of the difference, ignoring correlation between them: $\\sqrt{0.3^2 + 0.3^2} = \\sqrt{0.18} \\approx 0.42$.\n", "\n", "**Step 3.** $z = 1.0/0.42 \\approx 2.4$, so an over-identification (J) test would likely reject.\n", "\n", "**Step 4.** Two explanations: one instrument is invalid; or both are valid but effects are heterogeneous and each identifies a LATE for different compliers (Day 63).\n", "\n", "**Answer:** likely rejection, but the J test cannot tell the two explanations apart.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · Pooling two evaluations is GMM\n", "\n", "Two independent RCTs of the same microcredit programme estimate its effect on business profits: $\\hat\\theta_1$ with variance $v_1$ and $\\hat\\theta_2$ with variance $v_2$. Treat $g = (\\hat\\theta_1 - \\theta,\\; \\hat\\theta_2 - \\theta)$ as moment conditions with $W = \\operatorname{diag}(1/v_1, 1/v_2)$. Show that GMM gives the inverse-variance average.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Objective: $Q(\\theta) = (\\hat\\theta_1 - \\theta)^2/v_1 + (\\hat\\theta_2 - \\theta)^2/v_2$.\n", "\n", "**Step 2.** FOC: $-2(\\hat\\theta_1 - \\theta)/v_1 - 2(\\hat\\theta_2 - \\theta)/v_2 = 0$.\n", "\n", "**Step 3.** Collect terms in $\\theta$: $\\theta(1/v_1 + 1/v_2) = \\hat\\theta_1/v_1 + \\hat\\theta_2/v_2$.\n", "\n", "**Answer:**\n", "\n", "$$\n", "\\hat\\theta = \\frac{\\hat\\theta_1/v_1 + \\hat\\theta_2/v_2}{1/v_1 + 1/v_2},\n", "$$\n", "\n", "Day 8's pooled estimator. Here $W = S^{-1}$, the optimal choice.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · A numerical demand elasticity\n", "\n", "A demand model is too complicated to differentiate by hand, so you compute derivatives numerically. Test the method on a case with a known answer: $f(p) = e^{p}$ at $p = 1$ (true derivative $e$). Use central differences $[f(p + h) - f(p - h)]/2h$ with $h = 10^{-1}$, $10^{-5}$ and $10^{-12}$. Which step is best?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Truncation error of the central difference is about $h^2f'''(p)/6$: roughly $0.01 \\times 2.718/6 \\approx 4.5 \\times 10^{-3}$ at $h = 0.1$, and about $5 \\times 10^{-11}$ at $h = 10^{-5}$.\n", "\n", "**Step 2.** Rounding error is about (machine precision ≈ $10^{-16}$) × $f/h$: negligible at $h = 10^{-5}$, but about $10^{-4}$ at $h = 10^{-12}$, because $f(p + h) - f(p - h)$ subtracts two almost equal numbers.\n", "\n", "**Step 3.** The code confirms errors of $4.5 \\times 10^{-3}$, $6 \\times 10^{-11}$ and $2 \\times 10^{-4}$.\n", "\n", "**Answer:** $h = 10^{-5}$. Too large a step is biased; too small a step is swamped by rounding.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · From Phase 2 to Phase 3\n", "\n", "Name the Phase 2 tool you would use for each task: (a) check whether a treatment dummy is collinear with fixed effects; (b) compute the SE of the sum of three event-study coefficients; (c) estimate an effect controlling for 50 covariates and explain what \"controlling\" does; (d) get valid SEs when treatment is assigned by school; (e) estimate a model with more instruments than endogenous regressors.\n", "\n", "
Show solution\n", "\n", "**Step 1.** (a) Is treatment in the span of the fixed effects? → **rank** (Day 18).\n", "\n", "**Step 2.** (b) Variance of a sum of coefficients → $\\mathbf 1^{\\top}V\\mathbf 1$ (Day 22).\n", "\n", "**Step 3.** (c) Many controls and what \"controlling\" does → **OLS and FWL** (Days 27–28).\n", "\n", "**Step 4.** (d) Treatment assigned by school → **cluster-robust sandwich** at the school level (Day 29).\n", "\n", "**Step 5.** (e) More instruments than endogenous regressors → **GMM/2SLS with a J test** (Day 30).\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", "from scipy import optimize\n", "\n", "Q = lambda l: (3 - l)**2 + (13 - l - l**2)**2\n", "lam_I = optimize.minimize_scalar(Q, bounds=(0.1, 10), method=\"bounded\").x\n", "print(f\"GMM with W = I: lambda = {lam_I:.4f}\")\n", "\n", "# Optimal W for Poisson moments at lambda = 3: S = Var(X, X^2)\n", "l = 3.0\n", "EX, EX2, EX3, EX4 = l, l + l**2, l**3 + 3 * l**2 + l, l**4 + 6 * l**3 + 7 * l**2 + l\n", "S = np.array([[EX2 - EX**2, EX3 - EX * EX2], [EX3 - EX * EX2, EX4 - EX2**2]])\n", "W = np.linalg.inv(S)\n", "Qopt = lambda l: (g := np.array([3 - l, 13 - l - l**2])) @ W @ g\n", "print(f\"GMM with optimal W: lambda = {optimize.minimize_scalar(Qopt, bounds=(0.1, 10), method='bounded').x:.4f}\")\n", "\n", "# Newton-Raphson for an intercept-only logit with ybar = 0.3\n", "b = 0.0\n", "for step in range(4):\n", " p = 1 / (1 + np.exp(-b)); b = b + (0.3 - p) / (p * (1 - p)); print(f\"step {step + 1}: {b:.6f}\")\n", "print(\"exact:\", round(np.log(0.3 / 0.7), 6))\n", "\n", "for h in (1e-1, 1e-5, 1e-12):\n", " print(f\"h={h:g}: error {abs((np.exp(1 + h) - np.exp(1 - h)) / (2 * h) - np.e):.2e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- OLS, IV, 2SLS and MLE are all method-of-moments estimators; GMM handles extra moments by weighting.\n", "- Optimal weighting trusts precise moments more; the J test checks whether over-identifying moments agree.\n", "- Newton's method converges fast near the optimum; numerical derivatives need a sensible step size." ] } ] }