{ "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 27: OLS in Matrix Form\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "The linear model, the OLS estimator, its mean and variance, s², standard errors and the Gauss–Markov theorem, all in a few lines of matrix algebra.\n", "\n", "**Learning outcomes**\n", "\n", "- Compute β̂ = (X′X)⁻¹X′y, residuals, s² and standard errors by hand for a small data set.\n", "- Prove that OLS is unbiased with Var(β̂ | X) = σ²(X′X)⁻¹ and that E[e′e] = σ²(n − k).\n", "- Sketch the Gauss–Markov theorem and derive the omitted-variable-bias formula in matrix form.\n", "- Interpret dummy, interaction and log coefficients and see how covariates sharpen experiments." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### The model and the estimator\n", "\n", "$$\n", "y = X\\beta + \\varepsilon, \\qquad \\hat\\beta = (X^{\\top}X)^{-1}X^{\\top}y.\n", "$$\n", "\n", "Assumptions (the classical set):\n", "\n", "1. **Full rank:** $X$ has rank $k$ (Day 18).\n", "2. **Exogeneity:** $E[\\varepsilon \\mid X] = 0$.\n", "3. **Spherical errors:** $\\operatorname{Var}(\\varepsilon \\mid X) = \\sigma^2I$ (homoskedastic, uncorrelated).\n", "\n", "### Properties\n", "\n", "Substituting the model into the estimator gives the key identity\n", "\n", "$$\n", "\\hat\\beta = \\beta + (X^{\\top}X)^{-1}X^{\\top}\\varepsilon.\n", "$$\n", "\n", "- Under 1–2, $E[\\hat\\beta \\mid X] = \\beta$ (**unbiased**).\n", "- Under 1–3, $\\operatorname{Var}(\\hat\\beta \\mid X) = \\sigma^2(X^{\\top}X)^{-1}$.\n", "- $s^2 = e^{\\top}e/(n - k)$ is unbiased for $\\sigma^2$; standard errors are the square roots of the diagonal of $s^2(X^{\\top}X)^{-1}$, and $t_j = \\hat\\beta_j/\\operatorname{SE}_j$.\n", "- **Gauss–Markov:** under 1–3, OLS has the smallest variance among all linear unbiased estimators (BLUE).\n", "\n", "Assumption 3 is often doubtful; Day 29 replaces it with robust standard errors. Assumption 2 is the causal one: it fails whenever an omitted variable is correlated with $X$.\n", "\n", "### Omitted-variable bias, in matrices\n", "\n", "If the truth is $y = X_1\\beta_1 + X_2\\beta_2 + \\varepsilon$ but you regress on $X_1$ only,\n", "\n", "$$\n", "E[\\hat\\beta_1 \\mid X] = \\beta_1 + \\underbrace{(X_1^{\\top}X_1)^{-1}X_1^{\\top}X_2}_{\\text{slopes of } X_2 \\text{ on } X_1}\\beta_2.\n", "$$\n", "\n", "### Reading coefficients\n", "\n", "- Dummy: difference in conditional means between the groups.\n", "- Interaction $T \\times F$: how the effect of $T$ differs between groups.\n", "- Log outcome: approximately a proportional effect (Day 23 for the exact conversion).\n", "\n", "## Part 1 · Warm-up\n", "\n", "Data: five workers with $x = (0, 1, 2, 3, 4)$ years of experience earn $y = (1, 3, 2, 5, 4)$ (in \\$10/hour). Regress wages on an intercept and experience." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · The estimator\n", "\n", "Compute $X^{\\top}X$, $X^{\\top}y$ and $\\hat\\beta$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Sums: $n = 5$, $\\sum x = 10$, $\\sum x^2 = 30$, $\\sum y = 15$, $\\sum xy = 0 + 3 + 4 + 15 + 16 = 38$.\n", "\n", "**Step 2.** Cross-products:\n", "\n", "$$\n", "X^{\\top}X = \\begin{pmatrix}5 & 10\\\\10 & 30\\end{pmatrix}, \\qquad X^{\\top}y = \\begin{pmatrix}15\\\\38\\end{pmatrix}.\n", "$$\n", "\n", "**Step 3.** $\\det X^{\\top}X = 150 - 100 = 50$, so $(X^{\\top}X)^{-1} = \\frac{1}{50}\\begin{pmatrix}30 & -10\\\\-10 & 5\\end{pmatrix}$.\n", "\n", "**Step 4.** $\\hat\\beta = \\frac{1}{50}\\begin{pmatrix}30(15) - 10(38)\\\\-10(15) + 5(38)\\end{pmatrix} = \\frac{1}{50}\\begin{pmatrix}70\\\\40\\end{pmatrix} = \\begin{pmatrix}1.4\\\\0.8\\end{pmatrix}$.\n", "\n", "**Answer:** intercept 1.4 (\\$14/hour with no experience) and slope 0.8 (\\$8/hour per year of experience).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · Fitted values and residuals\n", "\n", "Compute $\\hat y$ and $e$, and check $X^{\\top}e = 0$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Fitted values $1.4 + 0.8x$: $1.4, 2.2, 3.0, 3.8, 4.6$.\n", "\n", "**Step 2.** Residuals $y - \\hat y$: $-0.4, 0.8, -1.0, 1.2, -0.6$.\n", "\n", "**Step 3.** $\\sum e = 0$ and $\\sum xe = 0 + 0.8 - 2.0 + 3.6 - 2.4 = 0$. ✓\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · The error variance\n", "\n", "Compute $s^2$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $e^{\\top}e = 0.16 + 0.64 + 1.00 + 1.44 + 0.36 = 3.6$.\n", "\n", "**Step 2.** Degrees of freedom $n - k = 5 - 2 = 3$.\n", "\n", "**Answer:** $s^2 = 3.6/3 = 1.2$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · Standard errors and t\n", "\n", "Compute $\\widehat{\\operatorname{Var}}(\\hat\\beta)$, the SE of the slope and its t-statistic.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\widehat{\\operatorname{Var}}(\\hat\\beta) = s^2(X^{\\top}X)^{-1} = 1.2 \\times \\frac{1}{50}\\begin{pmatrix}30 & -10\\\\-10 & 5\\end{pmatrix} = \\begin{pmatrix}0.72 & -0.24\\\\-0.24 & 0.12\\end{pmatrix}$.\n", "\n", "**Step 2.** SE of the slope $= \\sqrt{0.12} \\approx 0.346$.\n", "\n", "**Step 3.** $t = 0.8/0.346 \\approx 2.31$; with 3 degrees of freedom, $p \\approx 0.10$.\n", "\n", "**Answer:** SE ≈ 0.35 and $t \\approx 2.3$: suggestive, but five workers are too few for significance at 5%.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · A dummy regressor\n", "\n", "Show that in $y = \\beta_0 + \\beta_1T + \\varepsilon$, $\\hat\\beta_1 = \\bar y_1 - \\bar y_0$.\n", "\n", "
Show solution\n", "\n", "Fitted values are group means (Day 20, Exercise 14): $\\hat\\beta_0 = \\bar y_0$ and $\\hat\\beta_0 + \\hat\\beta_1 = \\bar y_1$. So $\\hat\\beta_1 = \\bar y_1 - \\bar y_0$.\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · Unbiasedness\n", "\n", "Prove $E[\\hat\\beta \\mid X] = \\beta$ under assumptions 1–2.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\hat\\beta = (X^{\\top}X)^{-1}X^{\\top}(X\\beta + \\varepsilon) = \\beta + (X^{\\top}X)^{-1}X^{\\top}\\varepsilon$.\n", "\n", "**Step 2.** Conditional on $X$, the matrix $(X^{\\top}X)^{-1}X^{\\top}$ is fixed, so $E[\\hat\\beta \\mid X] = \\beta + (X^{\\top}X)^{-1}X^{\\top}E[\\varepsilon \\mid X] = \\beta$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · The variance\n", "\n", "Prove $\\operatorname{Var}(\\hat\\beta \\mid X) = \\sigma^2(X^{\\top}X)^{-1}$ under assumptions 1–3.\n", "\n", "
Show solution\n", "\n", "With $A = (X^{\\top}X)^{-1}X^{\\top}$, Day 22 gives $\\operatorname{Var}(A\\varepsilon \\mid X) = A(\\sigma^2I)A^{\\top} = \\sigma^2(X^{\\top}X)^{-1}X^{\\top}X(X^{\\top}X)^{-1} = \\sigma^2(X^{\\top}X)^{-1}$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · Why divide by n − k\n", "\n", "Prove $E[e^{\\top}e \\mid X] = \\sigma^2(n - k)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $e = My = M\\varepsilon$ (since $MX = 0$), so $e^{\\top}e = \\varepsilon^{\\top}M\\varepsilon$ ($M$ is symmetric and idempotent).\n", "\n", "**Step 2.** A scalar equals its trace: $E[\\varepsilon^{\\top}M\\varepsilon] = E[\\operatorname{tr}(M\\varepsilon\\varepsilon^{\\top})] = \\operatorname{tr}(M\\sigma^2I) = \\sigma^2(n - k)$ (Day 20).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · Gauss–Markov\n", "\n", "Let $\\tilde\\beta = Cy$ be any linear unbiased estimator and write $C = (X^{\\top}X)^{-1}X^{\\top} + D$. Show $DX = 0$ and $\\operatorname{Var}(\\tilde\\beta) - \\operatorname{Var}(\\hat\\beta) = \\sigma^2DD^{\\top}$ (PSD).\n", "\n", "
Show solution\n", "\n", "**Step 1.** $E[\\tilde\\beta \\mid X] = CX\\beta = \\beta + DX\\beta$ for all $\\beta$ requires $DX = 0$.\n", "\n", "**Step 2.** Then $\\operatorname{Var}(\\tilde\\beta) = \\sigma^2CC^{\\top} = \\sigma^2[(X^{\\top}X)^{-1} + DX(X^{\\top}X)^{-1} + (X^{\\top}X)^{-1}X^{\\top}D^{\\top} + DD^{\\top}] = \\sigma^2(X^{\\top}X)^{-1} + \\sigma^2DD^{\\top}$, and $DD^{\\top}$ is PSD.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Back to Cov/Var\n", "\n", "Show that the slope from the matrix formula, in a model with an intercept and one regressor, is $\\sum(x_i - \\bar x)(y_i - \\bar y)/\\sum(x_i - \\bar x)^2$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** For $X = (\\mathbf 1, x)$: $(X^{\\top}X)^{-1} = \\dfrac{1}{n\\sum x^2 - (\\sum x)^2}\\begin{pmatrix}\\sum x^2 & -\\sum x\\\\-\\sum x & n\\end{pmatrix}$.\n", "\n", "**Step 2.** The slope is the second row times $X^{\\top}y = (\\sum y, \\sum xy)^{\\top}$:\n", "\n", "$$\n", "\\hat\\beta_1 = \\frac{n\\sum xy - \\sum x\\sum y}{n\\sum x^2 - (\\sum x)^2}.\n", "$$\n", "\n", "**Step 3.** Divide top and bottom by $n$ and use $\\sum x = n\\bar x$, $\\sum y = n\\bar y$: $\\hat\\beta_1 = \\dfrac{\\sum xy - n\\bar x\\bar y}{\\sum x^2 - n\\bar x^2} = \\dfrac{\\sum(x_i - \\bar x)(y_i - \\bar y)}{\\sum(x_i - \\bar x)^2}$.\n", "\n", "**Step 4.** Check with Exercise 1: $(38 - 5 \\cdot 2 \\cdot 3)/(30 - 5 \\cdot 4) = 8/10 = 0.8$. ✓\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · Interactions\n", "\n", "A wage regression gives $\\hat y = 20 + 3T + 2F - 1.5\\,T\\!\\times\\!F$ ($T$ = trained, $F$ = female). What is the estimated training effect for men and for women?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Effect of training: $\\partial\\hat y/\\partial T = 3 - 1.5F$.\n", "\n", "**Step 2.** Men ($F = 0$): 3. Women ($F = 1$): $3 - 1.5 = 1.5$.\n", "\n", "**Answer:** \\$30/hour for men and \\$15/hour for women (wages in \\$10). The interaction coefficient is the *difference* in effects; its SE tests whether the effects differ.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Omitted ability in matrix form\n", "\n", "The true wage model is $y = X_1\\beta_1 + X_2\\beta_2 + \\varepsilon$, where $X_1$ holds schooling and experience and $X_2$ holds unmeasured ability, with $E[\\varepsilon \\mid X] = 0$. You regress $y$ on $X_1$ only. Derive $E[\\hat\\beta_1^{\\text{short}} \\mid X]$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Short-regression estimator: $\\hat\\beta_1^{\\text{short}} = (X_1^{\\top}X_1)^{-1}X_1^{\\top}y$.\n", "\n", "**Step 2.** Substitute the true model:\n", "\n", "$$\n", "\\hat\\beta_1^{\\text{short}} = \\beta_1 + (X_1^{\\top}X_1)^{-1}X_1^{\\top}X_2\\beta_2 + (X_1^{\\top}X_1)^{-1}X_1^{\\top}\\varepsilon.\n", "$$\n", "\n", "**Step 3.** Take expectations given $X$; the last term has mean zero.\n", "\n", "**Answer:** $E[\\hat\\beta_1^{\\text{short}} \\mid X] = \\beta_1 + (X_1^{\\top}X_1)^{-1}X_1^{\\top}X_2\\beta_2$: (slopes of ability on schooling and experience) × (effect of ability). The bias vanishes only if $\\beta_2 = 0$ or $X_1^{\\top}X_2 = 0$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · Designing a fertiliser trial for precision\n", "\n", "An agronomist has 20 plots to estimate the yield response to fertiliser, and the slope's variance is $\\sigma^2/\\sum(x_i - \\bar x)^2$. Compare two designs: (A) half the plots at 0 kg and half at 10 kg; (B) two plots at each dose 0, 1, …, 9 kg. What is the trade-off?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Design A: $\\bar x = 5$ and every plot is 5 kg from the mean, so $\\sum(x - \\bar x)^2 = 20 \\times 25 = 500$.\n", "\n", "**Step 2.** Design B: $\\bar x = 4.5$; for 0–9 once, $\\sum(x - 4.5)^2 = 82.5$, so twice gives 165.\n", "\n", "**Step 3.** Variance ratio: $500/165 \\approx 3$.\n", "\n", "**Answer:** design A estimates the slope with about one-third of the variance. But it assumes the response is linear; design B can reveal curvature (e.g. diminishing returns), which design A cannot.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · Regression adjustment in an experiment\n", "\n", "In an RCT, a baseline covariate explains 50% of the outcome variance. By how much does including it shrink the SE of the treatment effect, and what sample-size increase is that worth?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Adding the covariate cuts the residual variance from $\\sigma^2$ to $\\sigma^2(1 - R^2) = 0.5\\sigma^2$.\n", "\n", "**Step 2.** Randomisation makes $T$ nearly uncorrelated with the covariate, so the treatment column's partialled-out variance (Day 28) is essentially unchanged.\n", "\n", "**Step 3.** So $\\operatorname{Var}(\\hat\\tau)$ halves and the SE falls by $\\sqrt{0.5} \\approx 0.71$.\n", "\n", "**Step 4.** Halving the variance by sample size alone would need twice the sample.\n", "\n", "**Answer:** about 29% smaller SEs, worth as much as doubling the sample.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · Correlated errors break the formula\n", "\n", "Errors are correlated within villages. Is $\\hat\\beta$ still unbiased? Is $s^2(X^{\\top}X)^{-1}$ still right?\n", "\n", "
Show solution\n", "\n", "**Step 1.** Unbiasedness uses only full rank and $E[\\varepsilon \\mid X] = 0$ (Exercise 6), so $\\hat\\beta$ is still unbiased.\n", "\n", "**Step 2.** With $\\operatorname{Var}(\\varepsilon \\mid X) = \\Omega \\ne \\sigma^2I$, Day 22's rule gives $\\operatorname{Var}(\\hat\\beta \\mid X) = (X^{\\top}X)^{-1}X^{\\top}\\Omega X(X^{\\top}X)^{-1}$.\n", "\n", "**Step 3.** This differs from $\\sigma^2(X^{\\top}X)^{-1}$; with positive within-village correlation the classical SEs are usually too small.\n", "\n", "**Answer:** unbiased, but the usual SEs are wrong. Day 29 fixes this with cluster-robust variances.\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 stats\n", "\n", "x = np.array([0, 1, 2, 3, 4.]); y = np.array([1, 3, 2, 5, 4.])\n", "X = np.column_stack([np.ones(5), x]); n, k = X.shape\n", "XtX_inv = np.linalg.inv(X.T @ X); b = XtX_inv @ X.T @ y\n", "e = y - X @ b; s2 = e @ e / (n - k); V = s2 * XtX_inv\n", "se = np.sqrt(np.diag(V)); t = b / se\n", "print(\"beta:\", b, \" residuals:\", e.round(2), \" s2:\", round(s2, 3))\n", "print(\"SE:\", se.round(4), \" t:\", t.round(3), \" p:\", (2 * stats.t.sf(np.abs(t), n - k)).round(3))\n", "\n", "# Monte Carlo: unbiasedness and the variance formula\n", "rng = np.random.default_rng(27)\n", "draws = np.array([XtX_inv @ X.T @ (X @ [1.4, 0.8] + rng.normal(0, 1.1, n)) for _ in range(50_000)])\n", "print(\"mean of beta-hat:\", draws.mean(0).round(3), \" simulated Var(slope):\", draws[:, 1].var().round(4),\n", " \" formula:\", round(1.1**2 * XtX_inv[1, 1], 4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- $\\hat\\beta = \\beta + (X^{\\top}X)^{-1}X^{\\top}\\varepsilon$ is the identity behind every OLS property.\n", "- Exogeneity gives unbiasedness; spherical errors give the textbook variance and Gauss–Markov.\n", "- Spreading $x$ and adding predictive covariates both buy precision." ] } ] }