{
 "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 25: Unconstrained Optimisation & Convexity\n",
    "\n",
    "*Beginner to Ninja: Math for Causal Inference*\n",
    "\n",
    "First- and second-order conditions, why convexity makes optimisation easy, and the two algorithms behind every estimator.\n",
    "\n",
    "**Learning outcomes**\n",
    "\n",
    "- Find and classify critical points with first- and second-order conditions in one and several variables.\n",
    "- Prove that for a convex function every stationary point is a global minimum.\n",
    "- Run gradient descent and Newton's method by hand and state their convergence conditions.\n",
    "- Solve profit- and monopoly-pricing problems and recognise when an MLE does not exist."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "### Conditions for an optimum\n",
    "\n",
    "For a smooth $f: \\mathbb{R}^k \\to \\mathbb{R}$:\n",
    "\n",
    "- **First-order condition (FOC):** at an interior optimum, $\\nabla f(x^*) = 0$.\n",
    "- **Second-order condition (SOC):** if also $H(x^*)$ is positive definite, $x^*$ is a strict local minimum; negative definite, a strict local maximum; indefinite, a **saddle point**.\n",
    "\n",
    "### Convexity\n",
    "\n",
    "A function is **convex** if the chord lies above the graph:\n",
    "\n",
    "$$\n",
    "f(\\lambda x + (1 - \\lambda)y) \\le \\lambda f(x) + (1 - \\lambda)f(y), \\qquad \\lambda \\in [0, 1].\n",
    "$$\n",
    "\n",
    "For smooth functions, convex $\\iff$ the tangent plane lies below: $f(y) \\ge f(x) + \\nabla f(x)^{\\top}(y - x)$ $\\iff$ $H(x)$ is PSD everywhere. The payoff: **for a convex function, any point with zero gradient is a global minimum**, and strict convexity makes it unique. Least squares, ridge, lasso and the (negative) logit log-likelihood are all convex.\n",
    "\n",
    "### Algorithms\n",
    "\n",
    "When the FOC cannot be solved by hand:\n",
    "\n",
    "- **Gradient descent:** $x_{t+1} = x_t - \\eta\\nabla f(x_t)$. Cheap per step; needs a good step size $\\eta$. On a quadratic with Hessian eigenvalues $\\lambda_j$ it converges iff $\\eta < 2/\\lambda_{\\max}$, and slowly when the condition number is large.\n",
    "- **Newton's method:** $x_{t+1} = x_t - H(x_t)^{-1}\\nabla f(x_t)$. It minimises the local second-order Taylor approximation (Day 23) at each step; it solves a quadratic in one step and converges very fast near the optimum, but each step needs the Hessian.\n",
    "\n",
    "For maximum likelihood, Newton's method applied to $\\ell$ is the classic **Newton–Raphson** (for the logit it is called iteratively reweighted least squares).\n",
    "\n",
    "## Part 1 · Warm-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1 · A parabola\n",
    "\n",
    "Minimise $f(x) = x^2 - 6x + 11$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$f'(x) = 2x - 6 = 0 \\Rightarrow x = 3$; $f''= 2 > 0$, a minimum with $f(3) = 2$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 · A cubic\n",
    "\n",
    "Find and classify the critical points of $f(x) = x^3 - 3x$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $3x^2 - 3 = 0 \\Rightarrow x = \\pm1$.\n",
    "\n",
    "**Step 2.** $f'' = 6x$: local minimum at $x = 1$, local maximum at $x = -1$.\n",
    "\n",
    "**Step 3.** Neither is global: $f$ is unbounded in both directions.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3 · Two variables\n",
    "\n",
    "Minimise $f(x, y) = x^2 + xy + y^2 - 3x$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** FOC: $f_x = 2x + y - 3 = 0$ and $f_y = x + 2y = 0$.\n",
    "\n",
    "**Step 2.** From the second, $x = -2y$; substitute: $-4y + y = 3 \\Rightarrow y = -1$, $x = 2$.\n",
    "\n",
    "**Step 3.** Hessian $\\begin{pmatrix}2 & 1\\\\1 & 2\\end{pmatrix}$ has eigenvalues 3 and 1 (Day 21): PD, so $f$ is strictly convex and $(2, -1)$ is the global minimum.\n",
    "\n",
    "**Step 4.** $f(2, -1) = 4 - 2 + 1 - 6 = -3$.\n",
    "\n",
    "**Answer:** minimum −3 at $(2, -1)$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4 · A saddle\n",
    "\n",
    "Classify the critical point of $f(x, y) = x^2 - y^2$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$\\nabla f = (2x, -2y) = 0$ at the origin; $H = \\operatorname{diag}(2, -2)$ is indefinite: a saddle (minimum along $x$, maximum along $y$).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5 · Step sizes\n",
    "\n",
    "Run two steps of gradient descent on $f(x) = x^2$ from $x_0 = 4$ with $\\eta = 0.25$, then with $\\eta = 1$ and $\\eta = 1.1$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $f'(x) = 2x$, so $x_{t+1} = x_t - \\eta(2x_t) = (1 - 2\\eta)x_t$.\n",
    "\n",
    "**Step 2.** $\\eta = 0.25$: factor 0.5, so $4 \\to 2 \\to 1$ (converges).\n",
    "\n",
    "**Step 3.** $\\eta = 1$: factor −1, so $4 \\to -4 \\to 4$ (oscillates forever).\n",
    "\n",
    "**Step 4.** $\\eta = 1.1$: factor −1.2, so $4 \\to -4.8 \\to 5.76$ (diverges).\n",
    "\n",
    "**Answer:** here the Hessian is $\\lambda = 2$, and convergence needs $|1 - 2\\eta| < 1$, i.e. $\\eta < 2/\\lambda = 1$ (Exercise 10).\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 2 · Proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6 · Stationary means optimal under convexity\n",
    "\n",
    "Prove that if $f$ is convex and differentiable and $\\nabla f(x^*) = 0$, then $x^*$ is a global minimum.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "By the tangent-plane characterisation, for every $y$: $f(y) \\ge f(x^*) + \\nabla f(x^*)^{\\top}(y - x^*) = f(x^*)$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7 · Building convex functions\n",
    "\n",
    "Prove that the sum of two convex functions is convex, and that $g(x) = f(Ax + b)$ is convex when $f$ is.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Sum: add the two chord inequalities.\n",
    "\n",
    "**Step 2.** Composition: $g(\\lambda x + (1 - \\lambda)y) = f(\\lambda(Ax + b) + (1 - \\lambda)(Ay + b)) \\le \\lambda f(Ax + b) + (1 - \\lambda)f(Ay + b)$.\n",
    "\n",
    "**Step 3.** Hence $\\|y - X\\beta\\|^2 + \\lambda\\|\\beta\\|^2$ is convex in $\\beta$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8 · Uniqueness\n",
    "\n",
    "Prove that a strictly convex function has at most one minimiser.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "Suppose $x \\ne y$ both attain the minimum $m$. Strict convexity gives $f(\\tfrac12x + \\tfrac12y) < \\tfrac12m + \\tfrac12m = m$, contradicting minimality.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 9 · Newton solves quadratics in one step\n",
    "\n",
    "Show that Newton's method reaches the minimiser of $f(x) = \\frac12x^{\\top}Ax - b^{\\top}x$ ($A$ PD) in one step from any start.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\nabla f = Ax - b$, $H = A$.\n",
    "\n",
    "**Step 2.** Then $x_1 = x_0 - A^{-1}(Ax_0 - b) = A^{-1}b$, which solves $\\nabla f = 0$.\n",
    "\n",
    "**Step 3.** So Newton applied to least squares returns $(X^{\\top}X)^{-1}X^{\\top}y$ immediately.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 10 · When gradient descent converges\n",
    "\n",
    "For $f(x) = \\frac12x^{\\top}Ax$ with $A$ PD, show gradient descent converges to 0 iff $\\eta < 2/\\lambda_{\\max}$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $\\nabla f = Ax$, so $x_{t+1} = x_t - \\eta Ax_t = (I - \\eta A)x_t$, and by induction $x_t = (I - \\eta A)^tx_0$.\n",
    "\n",
    "**Step 2.** If $Aq_j = \\lambda_jq_j$, then $(I - \\eta A)q_j = (1 - \\eta\\lambda_j)q_j$: same eigenvectors, eigenvalues $1 - \\eta\\lambda_j$.\n",
    "\n",
    "**Step 3.** Write $x_0 = \\sum_jc_jq_j$. Then $x_t = \\sum_jc_j(1 - \\eta\\lambda_j)^tq_j$.\n",
    "\n",
    "**Step 4.** $x_t \\to 0$ for every start iff $|1 - \\eta\\lambda_j| < 1$ for all $j$, i.e. $0 < \\eta < 2/\\lambda_j$ for all $j$, i.e. $\\eta < 2/\\lambda_{\\max}$.\n",
    "\n",
    "**Step 5.** Speed: the slowest direction shrinks by $|1 - \\eta\\lambda_{\\min}|$ per step. With $\\eta$ near $2/\\lambda_{\\max}$ and $\\lambda_{\\min} \\ll \\lambda_{\\max}$ (ill-conditioning), this is close to 1, so progress is slow.\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 3 · Applied"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 11 · A competitive firm\n",
    "\n",
    "A firm sells at $p = 10$, produces $q = 2L^{0.5}$ and pays wage $w = 5$. Find the profit-maximising labour, output and profit.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Profit: $\\pi(L) = p\\,q - wL = 10(2L^{0.5}) - 5L = 20L^{0.5} - 5L$.\n",
    "\n",
    "**Step 2.** FOC: $\\pi'(L) = 10L^{-0.5} - 5 = 0 \\Rightarrow L^{-0.5} = 0.5 \\Rightarrow L = 4$.\n",
    "\n",
    "**Step 3.** SOC: $\\pi''(L) = -5L^{-1.5} < 0$: a maximum.\n",
    "\n",
    "**Step 4.** Output $q = 2\\sqrt4 = 4$; profit $= 40 - 20 = 20$.\n",
    "\n",
    "**Answer:** hire 4 units of labour, produce 4, earn profit 20 (the value of the marginal product equals the wage).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 12 · Monopoly pricing\n",
    "\n",
    "Inverse demand $P = 100 - 2Q$, cost $C = 20Q$. Find the monopoly quantity, price and profit.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Profit: $\\pi(Q) = (100 - 2Q)Q - 20Q = 80Q - 2Q^2$.\n",
    "\n",
    "**Step 2.** FOC (marginal revenue = marginal cost): $80 - 4Q = 0 \\Rightarrow Q = 20$. SOC: $\\pi'' = -4 < 0$.\n",
    "\n",
    "**Step 3.** Price: $P = 100 - 2(20) = 60$.\n",
    "\n",
    "**Step 4.** Profit: $(60 - 20) \\times 20 = 800$.\n",
    "\n",
    "**Answer:** $Q = 20$, $P = 60$, profit 800.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 13 · Newton–Raphson for a Poisson MLE\n",
    "\n",
    "Ten counts have $\\sum x_i = 30$. Starting from $\\lambda_0 = 1$, perform two Newton–Raphson steps on $\\ell(\\lambda) = \\sum x_i\\ln\\lambda - n\\lambda$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Derivatives: $\\ell'(\\lambda) = 30/\\lambda - 10$ and $\\ell''(\\lambda) = -30/\\lambda^2$. Newton update: $\\lambda \\leftarrow \\lambda - \\ell'(\\lambda)/\\ell''(\\lambda)$.\n",
    "\n",
    "**Step 2.** From $\\lambda_0 = 1$: $\\ell' = 20$, $\\ell'' = -30$, so $\\lambda_1 = 1 + 20/30 \\approx 1.667$.\n",
    "\n",
    "**Step 3.** From $\\lambda_1 = 5/3$: $\\ell' = 18 - 10 = 8$, $\\ell'' = -30/(25/9) = -10.8$, so $\\lambda_2 = 1.667 + 8/10.8 \\approx 2.407$.\n",
    "\n",
    "**Step 4.** Further steps give 2.883, 2.995, 3.000 (code below).\n",
    "\n",
    "**Answer:** the iterates converge to $\\hat\\lambda = \\bar x = 3$, slowly at first and then very fast near the optimum.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 14 · Estimating a CES production function\n",
    "\n",
    "An economist fits a CES production function $\\ln Y_i = \\ln A + \\frac{\\nu}{\\rho}\\ln\\big(\\delta K_i^{\\rho} + (1 - \\delta)L_i^{\\rho}\\big) + \\varepsilon_i$ by non-linear least squares. The objective can have several local minima and flat regions (e.g. near $\\rho = 0$, the Cobb–Douglas limit). What should she do?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Recognise the risk: a gradient method stops at *a* point with zero gradient, which need not be the global minimum.\n",
    "\n",
    "**Step 2.** Start the optimiser from several points (or a grid over $\\rho$ and $\\delta$) and compare the objective values reached.\n",
    "\n",
    "**Step 3.** Reparametrise where possible (e.g. profile out $\\ln A$, which enters linearly) and check the curvature (Hessian) at the solution.\n",
    "\n",
    "**Answer:** use multiple starts, report sensitivity, and prefer formulations that are convex. This is one reason linear and logit models, whose objectives are convex, dominate applied work.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 15 · When the MLE does not exist\n",
    "\n",
    "In a logit of treatment on one covariate, every unit with $x > 5$ is treated and every unit with $x \\le 5$ is not. What happens to Newton–Raphson, and what does it mean for propensity-score methods?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Perfect separation: every unit with $x > 5$ is treated and every unit with $x \\le 5$ is not.\n",
    "\n",
    "**Step 2.** A steeper slope always fits better: the log-likelihood keeps increasing as $\\beta_1 \\to \\infty$, so there is no finite maximum.\n",
    "\n",
    "**Step 3.** Newton–Raphson therefore diverges, with fitted probabilities heading to exactly 0 and 1 (the code shows the slope growing every iteration).\n",
    "\n",
    "**Answer:** no MLE exists. For causal inference this is a failure of **overlap**: some covariate values have no comparison units, so effects there cannot be estimated without extrapolation (Day 42).\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",
    "# Gradient descent on a quadratic with condition number 100 vs Newton\n",
    "A = np.diag([1.0, 100.0]); grad = lambda x: A @ x\n",
    "x = np.array([1.0, 1.0]); eta = 1.9 / 100\n",
    "for t in range(200):\n",
    "    x = x - eta * grad(x)\n",
    "print(\"GD after 200 steps:\", x.round(4), \"  (slow along the small-eigenvalue direction)\")\n",
    "x0 = np.array([1.0, 1.0]); print(\"Newton after 1 step:\", x0 - np.linalg.solve(A, grad(x0)))\n",
    "\n",
    "# Newton-Raphson for the Poisson MLE (Exercise 13)\n",
    "lam, path = 1.0, []\n",
    "for _ in range(6):\n",
    "    lam = lam - (30 / lam - 10) / (-30 / lam**2); path.append(round(lam, 4))\n",
    "print(\"Newton path:\", path)\n",
    "\n",
    "# Perfect separation: the logit slope runs away\n",
    "x = np.arange(1, 11.); T = (x > 5).astype(float); X = np.column_stack([np.ones(10), x]); b = np.zeros(2)\n",
    "for it in range(8):\n",
    "    p = 1 / (1 + np.exp(-X @ b)); W = p * (1 - p)\n",
    "    b = b + np.linalg.solve(X.T @ (X * W[:, None]) + 1e-12 * np.eye(2), X.T @ (T - p))\n",
    "    print(f\"iteration {it + 1}: slope = {b[1]:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key takeaways\n",
    "\n",
    "- FOC finds candidates; the Hessian classifies them; convexity turns local into global.\n",
    "- Gradient descent is cheap but slowed by ill-conditioning; Newton uses curvature and converges fast near the optimum.\n",
    "- Non-existence of an MLE (separation) is a data problem, often an overlap problem, not a software problem."
   ]
  }
 ]
}