{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ebb22dbd",
   "metadata": {},
   "source": [
    "# The Overlapping Generations Model\n",
    "\n",
    "In this lecture we study the famous overlapping generations (OLG) model, which\n",
    "is used by policy makers and researchers to examine\n",
    "\n",
    "- fiscal policy  \n",
    "- monetary policy  \n",
    "- long-run growth  \n",
    "\n",
    "\n",
    "and many other topics.\n",
    "\n",
    "The first rigorous version of the OLG model was developed by Paul Samuelson\n",
    "[[Samuelson, 1958](https://intro.quantecon.org/zreferences.html#id23)].\n",
    "\n",
    "Our aim is to gain a good understanding of a simple version of the OLG\n",
    "model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b65f2f2a",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "The dynamics of the OLG model are quite similar to those of the [Solow-Swan\n",
    "growth model](https://intro.quantecon.org/solow.html).\n",
    "\n",
    "At the same time, the OLG model adds an important new feature: the choice of\n",
    "how much to save is endogenous.\n",
    "\n",
    "To see why this is important, suppose, for example, that we are interested in\n",
    "predicting the effect of a new tax on long-run growth.\n",
    "\n",
    "We could add a tax to the Solow-Swan model and look at the change in the\n",
    "steady state.\n",
    "\n",
    "But this ignores the fact that households will change their savings and\n",
    "consumption behavior when they face the new tax rate.\n",
    "\n",
    "Such changes can substantially alter the predictions of the model.\n",
    "\n",
    "Hence, if we care about accurate predictions, we should model the decision\n",
    "problems of the agents.\n",
    "\n",
    "In particular, households in the model should decide how much to save and how\n",
    "much to consume, given the environment that they face (technology, taxes,\n",
    "prices, etc.)\n",
    "\n",
    "The OLG model takes up this challenge.\n",
    "\n",
    "We will present a simple version of the OLG model that clarifies the decision\n",
    "problem of households and studies the implications for long-run growth.\n",
    "\n",
    "Let’s start with some imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0ac5ba5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import optimize\n",
    "from collections import namedtuple\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43ca1eb1",
   "metadata": {},
   "source": [
    "## Environment\n",
    "\n",
    "We assume that time is discrete, so that $ t=0, 1, \\ldots $.\n",
    "\n",
    "An individual born at time $ t $ lives for two periods, $ t $ and $ t + 1 $.\n",
    "\n",
    "We call an agent\n",
    "\n",
    "- “young” during the first period of their lives and  \n",
    "- “old” during the second period of their lives.  \n",
    "\n",
    "\n",
    "Young agents work, supplying labor and earning labor income.\n",
    "\n",
    "They also decide how much to save.\n",
    "\n",
    "Old agents do not work, so all income is financial.\n",
    "\n",
    "Their financial income is from interest on their savings from wage income,\n",
    "which is then combined with the labor of the new young generation at $ t+1 $.\n",
    "\n",
    "The wage and interest rates are determined in equilibrium by supply and\n",
    "demand.\n",
    "\n",
    "To make the algebra slightly easier, we are going to assume a constant\n",
    "population size.\n",
    "\n",
    "We normalize the constant population size in each period to 1.\n",
    "\n",
    "We also suppose that each agent supplies one “unit” of labor hours, so total\n",
    "labor supply is 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b9f29e",
   "metadata": {},
   "source": [
    "## Supply of capital\n",
    "\n",
    "First let’s consider the household side."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a03f65eb",
   "metadata": {},
   "source": [
    "### Consumer’s problem\n",
    "\n",
    "Suppose that utility for individuals born at time $ t $ takes the form\n",
    "\n",
    "\n",
    "<a id='equation-eq-crra'></a>\n",
    "$$\n",
    "U_t = u(c_t) + \\beta u(c_{t+1}) \\tag{26.1}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ u: \\mathbb R_+ \\to \\mathbb R $ is called the “flow” utility function  \n",
    "- $ \\beta \\in (0, 1) $ is the discount factor  \n",
    "- $ c_t $ is time $ t $ consumption of the individual born at time $ t $  \n",
    "- $ c_{t+1} $ is time $ t+1 $ consumption of the same individual  \n",
    "\n",
    "\n",
    "We assume that $ u $ is strictly increasing.\n",
    "\n",
    "Savings behavior is determined by the optimization problem\n",
    "\n",
    "\n",
    "<a id='equation-max-sav-olg'></a>\n",
    "$$\n",
    "\\max_{c_t, c_{t+1}} \n",
    "    \\,  \\left \\{ u(c_t) + \\beta u(c_{t+1}) \\right \\} \\tag{26.2}\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "$$\n",
    "c_t + s_t \\le w_t \n",
    "     \\quad \\text{and} \\quad\n",
    "     c_{t+1}   \\le R_{t+1} s_t\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ s_t $ is savings by an individual born at time $ t $  \n",
    "- $ w_t $ is the wage rate at time $ t $  \n",
    "- $ R_{t+1} $ is the gross interest rate on savings invested at time $ t $, paid at time $ t+1 $  \n",
    "\n",
    "\n",
    "Since $ u $ is strictly increasing, both of these constraints will hold as equalities at the maximum.\n",
    "\n",
    "Using this fact and substituting $ s_t $ from the first constraint into the second we get\n",
    "$ c_{t+1} = R_{t+1}(w_t - c_t) $.\n",
    "\n",
    "The first-order condition for a maximum can be obtained\n",
    "by plugging $ c_{t+1} $ into the objective function, taking the derivative\n",
    "with respect to $ c_t $, and setting it to zero.\n",
    "\n",
    "This leads to the **Euler equation** of the OLG model, which describes the optimal intertemporal consumption dynamics:\n",
    "\n",
    "\n",
    "<a id='equation-euler-1-olg'></a>\n",
    "$$\n",
    "u'(c_t) = \\beta R_{t+1}  u'( R_{t+1} (w_t - c_t)) \\tag{26.3}\n",
    "$$\n",
    "\n",
    "From the first constraint we get $ c_t = w_t - s_t $, so the Euler equation\n",
    "can also be expressed as\n",
    "\n",
    "\n",
    "<a id='equation-euler-2-olg'></a>\n",
    "$$\n",
    "u'(w_t - s_t) = \\beta R_{t+1}  u'( R_{t+1} s_t) \\tag{26.4}\n",
    "$$\n",
    "\n",
    "Suppose that, for each $ w_t $ and $ R_{t+1} $, there is exactly one $ s_t $ that\n",
    "solves [(26.4)](#equation-euler-2-olg).\n",
    "\n",
    "Then savings can be written as a fixed function of $ w_t $ and $ R_{t+1} $.\n",
    "\n",
    "We write this as\n",
    "\n",
    "\n",
    "<a id='equation-saving-1-olg'></a>\n",
    "$$\n",
    "s_t = s(w_t, R_{t+1}) \\tag{26.5}\n",
    "$$\n",
    "\n",
    "The precise form of the function $ s $ will depend on the choice of flow utility\n",
    "function $ u $.\n",
    "\n",
    "Together, $ w_t $ and $ R_{t+1} $ represent the *prices* in the economy (price of\n",
    "labor and rental rate of capital).\n",
    "\n",
    "Thus, [(26.5)](#equation-saving-1-olg) states the quantity of savings given prices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2907d419",
   "metadata": {},
   "source": [
    "### Example: log preferences\n",
    "\n",
    "In the special case $ u(c) = \\log c $, the Euler equation simplifies to\n",
    "$ s_t= \\beta (w_t - s_t) $.\n",
    "\n",
    "Solving for saving, we get\n",
    "\n",
    "\n",
    "<a id='equation-saving-log-2-olg'></a>\n",
    "$$\n",
    "s_t = s(w_t, R_{t+1}) = \\frac{\\beta}{1+\\beta} w_t \\tag{26.6}\n",
    "$$\n",
    "\n",
    "In this special case, savings does not depend on the interest rate."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e84697d5",
   "metadata": {},
   "source": [
    "### Savings and investment\n",
    "\n",
    "Since the population size is normalized to 1, $ s_t $ is also total savings in\n",
    "the economy at time $ t $.\n",
    "\n",
    "In our closed economy, there is no foreign investment, so net savings equals\n",
    "total investment, which can be understood as supply of capital to firms.\n",
    "\n",
    "In the next section we investigate demand for capital.\n",
    "\n",
    "Equating supply and demand will allow us to determine equilibrium in the OLG\n",
    "economy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac99b83a",
   "metadata": {},
   "source": [
    "## Demand for capital\n",
    "\n",
    "First we describe the firm’s problem and then we write down an equation\n",
    "describing demand for capital given prices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf0a772",
   "metadata": {},
   "source": [
    "### Firm’s problem\n",
    "\n",
    "For each integer $ t \\geq 0 $, output $ y_t $ in period $ t $ is given by the\n",
    "**[Cobb-Douglas production function](https://en.wikipedia.org/wiki/Cobb%E2%80%93Douglas_production_function)**\n",
    "\n",
    "\n",
    "<a id='equation-cobb-douglas'></a>\n",
    "$$\n",
    "y_t = k_t^{\\alpha} \\ell_t^{1-\\alpha} \\tag{26.7}\n",
    "$$\n",
    "\n",
    "Here $ k_t $ is capital, $ \\ell_t $ is labor, and  $ \\alpha $ is a parameter\n",
    "(sometimes called the “output elasticity of capital”).\n",
    "\n",
    "The profit maximization problem of the firm is\n",
    "\n",
    "\n",
    "<a id='equation-opt-profit-olg'></a>\n",
    "$$\n",
    "\\max_{k_t, \\ell_t} \\{ k^{\\alpha}_t \\ell_t^{1-\\alpha} - R_t k_t -w_t \\ell_t  \\} \\tag{26.8}\n",
    "$$\n",
    "\n",
    "The first-order conditions are obtained by taking the derivative of the\n",
    "objective function with respect to capital and labor respectively and setting\n",
    "them to zero:\n",
    "\n",
    "$$\n",
    "(1-\\alpha)(k_t / \\ell_t)^{\\alpha} = w_t\n",
    "    \\quad \\text{and} \\quad\n",
    "    \\alpha (k_t / \\ell_t)^{\\alpha - 1} = R_t\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe96d18",
   "metadata": {},
   "source": [
    "### Demand\n",
    "\n",
    "Using our assumption $ \\ell_t = 1 $ allows us to write\n",
    "\n",
    "\n",
    "<a id='equation-wage-one'></a>\n",
    "$$\n",
    "w_t = (1-\\alpha)k_t^\\alpha \\tag{26.9}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-interest-rate-one'></a>\n",
    "$$\n",
    "R_t =\n",
    "    \\alpha k_t^{\\alpha - 1} \\tag{26.10}\n",
    "$$\n",
    "\n",
    "Rearranging [(26.10)](#equation-interest-rate-one) gives the aggregate demand for capital\n",
    "at time $ t+1 $\n",
    "\n",
    "\n",
    "<a id='equation-aggregate-demand-capital-olg'></a>\n",
    "$$\n",
    "k^d (R_{t+1}) \n",
    "    := \\left (\\frac{\\alpha}{R_{t+1}} \\right )^{1/(1-\\alpha)} \\tag{26.11}\n",
    "$$\n",
    "\n",
    "In Python code this is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e57d859a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def capital_demand(R, α):\n",
    "    return (α/R)**(1/(1-α)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4cbb4d5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def capital_supply(R, β, w):\n",
    "    R = np.ones_like(R)\n",
    "    return R * (β / (1 + β)) * w"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c01caa3e",
   "metadata": {},
   "source": [
    "The next figure plots the supply of capital, as in [(26.6)](#equation-saving-log-2-olg), as well as the demand for capital, as in [(26.11)](#equation-aggregate-demand-capital-olg), as functions of the interest rate $ R_{t+1} $.\n",
    "\n",
    "(For the special case of log utility, supply does not depend on the interest rate, so we have a constant function.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b86fd5b",
   "metadata": {},
   "source": [
    "## Equilibrium\n",
    "\n",
    "In this section we derive equilibrium conditions and investigate an example."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd65b475",
   "metadata": {},
   "source": [
    "### Equilibrium conditions\n",
    "\n",
    "In equilibrium, savings at time $ t $ equals investment at time $ t $, which\n",
    "equals capital supply at time $ t+1 $.\n",
    "\n",
    "Equilibrium is computed by equating these quantities, setting\n",
    "\n",
    "\n",
    "<a id='equation-equilibrium-1'></a>\n",
    "$$\n",
    "s(w_t, R_{t+1}) \n",
    "    = k^d(R_{t+1})\n",
    "    = \\left (\\frac{\\alpha}{R_{t+1}} \\right )^{1/(1-\\alpha)} \\tag{26.12}\n",
    "$$\n",
    "\n",
    "In principle, we can now solve for the equilibrium price $ R_{t+1} $ given $ w_t $.\n",
    "\n",
    "(In practice, we first need to specify the function $ u $ and hence $ s $.)\n",
    "\n",
    "When we solve this equation, which concerns time $ t+1 $ outcomes, time\n",
    "$ t $ quantities are already determined, so we can treat $ w_t $ as a constant.\n",
    "\n",
    "From equilibrium $ R_{t+1} $ and [(26.11)](#equation-aggregate-demand-capital-olg), we can obtain\n",
    "the equilibrium quantity $ k_{t+1} $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cb9bce6",
   "metadata": {},
   "source": [
    "### Example: log utility\n",
    "\n",
    "In the case of log utility, we can use [(26.12)](#equation-equilibrium-1) and [(26.6)](#equation-saving-log-2-olg) to obtain\n",
    "\n",
    "\n",
    "<a id='equation-equilibrium-2'></a>\n",
    "$$\n",
    "\\frac{\\beta}{1+\\beta} w_t\n",
    "    = \\left( \\frac{\\alpha}{R_{t+1}} \\right)^{1/(1-\\alpha)} \\tag{26.13}\n",
    "$$\n",
    "\n",
    "Solving for the equilibrium interest rate gives\n",
    "\n",
    "\n",
    "<a id='equation-equilibrium-price'></a>\n",
    "$$\n",
    "R_{t+1} = \n",
    "    \\alpha \n",
    "    \\left( \n",
    "        \\frac{\\beta}{1+\\beta} w_t\n",
    "    \\right)^{\\alpha-1} \\tag{26.14}\n",
    "$$\n",
    "\n",
    "In Python we can compute this via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81807f64",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def equilibrium_R_log_utility(α, β, w):\n",
    "    R = α * ( (β * w) / (1 + β))**(α - 1)\n",
    "    return R"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "448b1757",
   "metadata": {},
   "source": [
    "In the case of log utility, since capital supply does not depend on the interest rate, the equilibrium quantity is fixed by supply.\n",
    "\n",
    "That is,\n",
    "\n",
    "\n",
    "<a id='equation-equilibrium-quantity'></a>\n",
    "$$\n",
    "k_{t+1} = s(w_t, R_{t+1}) = \\frac{\\beta }{1+\\beta} w_t \\tag{26.15}\n",
    "$$\n",
    "\n",
    "Let’s redo our plot above but now inserting the equilibrium quantity and price."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7e1e8f6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "R_vals = np.linspace(0.3, 1)\n",
    "α, β = 0.5, 0.9\n",
    "w = 2.0\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(R_vals, capital_demand(R_vals, α), \n",
    "        label=\"aggregate demand\")\n",
    "ax.plot(R_vals, capital_supply(R_vals, β, w), \n",
    "        label=\"aggregate supply\")\n",
    "\n",
    "R_e = equilibrium_R_log_utility(α, β, w)\n",
    "k_e = (β / (1 + β)) * w\n",
    "\n",
    "ax.plot(R_e, k_e, 'o',label='equilibrium')\n",
    "\n",
    "ax.set_xlabel(\"$R_{t+1}$\")\n",
    "ax.set_ylabel(\"$k_{t+1}$\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c455b911",
   "metadata": {},
   "source": [
    "## Dynamics\n",
    "\n",
    "In this section we discuss dynamics.\n",
    "\n",
    "For now we will focus on the case of log utility, so that the equilibrium is determined by [(26.15)](#equation-equilibrium-quantity)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6eebdf3",
   "metadata": {},
   "source": [
    "### Evolution of capital\n",
    "\n",
    "The discussion above shows how equilibrium $ k_{t+1} $ is obtained given $ w_t $.\n",
    "\n",
    "From [(26.9)](#equation-wage-one) we can translate this into $ k_{t+1} $ as a function of $ k_t $\n",
    "\n",
    "In particular, since $ w_t = (1-\\alpha)k_t^\\alpha $, we have\n",
    "\n",
    "\n",
    "<a id='equation-law-of-motion-capital'></a>\n",
    "$$\n",
    "k_{t+1} = \\frac{\\beta}{1+\\beta} (1-\\alpha)(k_t)^{\\alpha} \\tag{26.16}\n",
    "$$\n",
    "\n",
    "If we iterate on this equation, we get a sequence for capital stock.\n",
    "\n",
    "Let’s plot the 45-degree diagram of these dynamics, which we write as\n",
    "\n",
    "$$\n",
    "k_{t+1} = g(k_t)\n",
    "    \\quad \\text{where }\n",
    "    g(k) := \\frac{\\beta}{1+\\beta} (1-\\alpha)(k)^{\\alpha}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c632c8b7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def k_update(k, α, β):\n",
    "    return β * (1 - α) * k**α /  (1 + β)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d904c9b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "α, β = 0.5, 0.9\n",
    "kmin, kmax = 0, 0.1\n",
    "n = 1000\n",
    "k_grid = np.linspace(kmin, kmax, n)\n",
    "k_grid_next = k_update(k_grid,α,β)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "ymin, ymax = np.min(k_grid_next), np.max(k_grid_next)\n",
    "\n",
    "ax.plot(k_grid, k_grid_next,  lw=2, alpha=0.6, label='$g$')\n",
    "ax.plot(k_grid, k_grid, 'k-', lw=1, alpha=0.7, label='$45^{\\circ}$')\n",
    "\n",
    "\n",
    "ax.legend(loc='upper left', frameon=False, fontsize=12)\n",
    "ax.set_xlabel('$k_t$', fontsize=12)\n",
    "ax.set_ylabel('$k_{t+1}$', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49faaf91",
   "metadata": {},
   "source": [
    "### Steady state (log case)\n",
    "\n",
    "The diagram shows that the model has a unique positive steady state, which we\n",
    "denote by $ k^* $.\n",
    "\n",
    "We can solve for $ k^* $ by setting $ k^* = g(k^*) $, or\n",
    "\n",
    "\n",
    "<a id='equation-steady-state-1'></a>\n",
    "$$\n",
    "k^* = \\frac{\\beta (1-\\alpha) (k^*)^{\\alpha}}{(1+\\beta)} \\tag{26.17}\n",
    "$$\n",
    "\n",
    "Solving this equation yields\n",
    "\n",
    "\n",
    "<a id='equation-steady-state-2'></a>\n",
    "$$\n",
    "k^* = \\left (\\frac{\\beta (1-\\alpha)}{1+\\beta} \\right )^{1/(1-\\alpha)} \\tag{26.18}\n",
    "$$\n",
    "\n",
    "We can get the steady state interest rate from [(26.10)](#equation-interest-rate-one), which yields\n",
    "\n",
    "$$\n",
    "R^* = \\alpha (k^*)^{\\alpha - 1} \n",
    "        = \\frac{\\alpha}{1 - \\alpha} \\frac{1 + \\beta}{\\beta}\n",
    "$$\n",
    "\n",
    "In Python we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53f4da8d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "k_star = ((β * (1 - α))/(1 + β))**(1/(1-α))\n",
    "R_star = (α/(1 - α)) * ((1 + β) / β)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac11ed6b",
   "metadata": {},
   "source": [
    "### Time series\n",
    "\n",
    "The 45-degree diagram above shows that time series of capital with positive initial conditions converge to this steady state.\n",
    "\n",
    "Let’s plot some time series that visualize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "633f8d60",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ts_length = 25\n",
    "k_series = np.empty(ts_length)\n",
    "k_series[0] = 0.02\n",
    "for t in range(ts_length - 1):\n",
    "    k_series[t+1] = k_update(k_series[t], α, β)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(k_series, label=\"capital series\")\n",
    "ax.plot(range(ts_length), np.full(ts_length, k_star), 'k--', label=\"$k^*$\")\n",
    "ax.set_ylim(0, 0.1)\n",
    "ax.set_ylabel(\"capital\")\n",
    "ax.set_xlabel(\"$t$\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35efa375",
   "metadata": {},
   "source": [
    "If you experiment with different positive initial conditions, you will see that the series always converges to $ k^* $.\n",
    "\n",
    "Below we also plot the gross interest rate over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a751694c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "R_series = α * k_series**(α - 1)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(R_series, label=\"gross interest rate\")\n",
    "ax.plot(range(ts_length), np.full(ts_length, R_star), 'k--', label=\"$R^*$\")\n",
    "ax.set_ylim(0, 4)\n",
    "ax.set_ylabel(\"gross interest rate\")\n",
    "ax.set_xlabel(\"$t$\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a4fd129",
   "metadata": {},
   "source": [
    "The interest rate reflects the marginal product of capital, which is high when capital stock is low."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d23076a",
   "metadata": {},
   "source": [
    "## CRRA preferences\n",
    "\n",
    "Previously, in our examples, we looked at the case of log utility.\n",
    "\n",
    "Log utility is a rather special case of CRRA utility with $ \\gamma \\to 1 $.\n",
    "\n",
    "In this section, we are going to assume that $ u(c) = \\frac{ c^{1-\n",
    "\\gamma}-1}{1-\\gamma} $, where $ \\gamma >0, \\gamma\\neq 1 $.\n",
    "\n",
    "This function is called the CRRA utility function.\n",
    "\n",
    "In other respects, the model is the same.\n",
    "\n",
    "Below we define the utility function in Python and construct a `namedtuple` to store the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6298fe56",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def crra(c, γ):\n",
    "    return c**(1 - γ) / (1 - γ)\n",
    "\n",
    "Model = namedtuple('Model', ['α',        # Cobb-Douglas parameter\n",
    "                             'β',        # discount factor\n",
    "                             'γ']        # parameter in CRRA utility\n",
    "                   )\n",
    "\n",
    "def create_olg_model(α=0.4, β=0.9, γ=0.5):\n",
    "    return Model(α=α, β=β, γ=γ)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "197473c7",
   "metadata": {},
   "source": [
    "Let’s also redefine the capital demand function to work with this `namedtuple`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59e5660f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def capital_demand(R, model):\n",
    "    return (α/R)**(1/(1-model.α)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fd5df32",
   "metadata": {},
   "source": [
    "### Supply\n",
    "\n",
    "For households, the Euler equation becomes\n",
    "\n",
    "\n",
    "<a id='equation-euler-crra'></a>\n",
    "$$\n",
    "(w_t - s_t)^{-\\gamma} = \\beta R^{1-\\gamma}_{t+1}  (s_t)^{-\\gamma} \\tag{26.19}\n",
    "$$\n",
    "\n",
    "Solving for savings, we have\n",
    "\n",
    "\n",
    "<a id='equation-saving-crra'></a>\n",
    "$$\n",
    "s_t \n",
    "    = s(w_t, R_{t+1}) \n",
    "    = w_t \\left [ \n",
    "        1 + \\beta^{-1/\\gamma} R_{t+1}^{(\\gamma-1)/\\gamma} \n",
    "      \\right ]^{-1} \\tag{26.20}\n",
    "$$\n",
    "\n",
    "Notice how, unlike the log case, savings now depends on the interest rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14c6113b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def savings_crra(w, R, model):\n",
    "    α, β, γ = model\n",
    "    return w / (1 + β**(-1/γ) * R**((γ-1)/γ)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f92ee63",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "model = create_olg_model()\n",
    "w = 2.0\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(R_vals, capital_demand(R_vals, model), \n",
    "        label=\"aggregate demand\")\n",
    "ax.plot(R_vals, savings_crra(w, R_vals, model), \n",
    "        label=\"aggregate supply\")\n",
    "\n",
    "ax.set_xlabel(\"$R_{t+1}$\")\n",
    "ax.set_ylabel(\"$k_{t+1}$\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbce094",
   "metadata": {},
   "source": [
    "### Equilibrium\n",
    "\n",
    "Equating aggregate demand for capital  (see [(26.11)](#equation-aggregate-demand-capital-olg))\n",
    "with our new aggregate supply function yields equilibrium capital.\n",
    "\n",
    "Thus, we set\n",
    "\n",
    "\n",
    "<a id='equation-equilibrium-crra-2'></a>\n",
    "$$\n",
    "w_t \\left [ 1 + \\beta^{-1/\\gamma} R_{t+1}^{(\\gamma-1)/\\gamma} \\right ]^{-1} \n",
    "    = \\left (\\frac{R_{t+1}}{\\alpha} \\right )^{1/(\\alpha - 1)} \\tag{26.21}\n",
    "$$\n",
    "\n",
    "This expression is quite complex and we cannot solve for $ R_{t+1} $ analytically.\n",
    "\n",
    "Combining [(26.10)](#equation-interest-rate-one) and [(26.21)](#equation-equilibrium-crra-2) yields\n",
    "\n",
    "\n",
    "<a id='equation-law-of-motion-capital-crra'></a>\n",
    "$$\n",
    "k_{t+1} = \\left [ 1 + \\beta^{-1/\\gamma} (\\alpha k^{\\alpha - 1}_{t+1})^{(\\gamma-1)/\\gamma} \\right ]^{-1} (1-\\alpha)(k_t)^{\\alpha} \\tag{26.22}\n",
    "$$\n",
    "\n",
    "Again, with this equation and $ k_t $ as given, we cannot solve for $ k_{t+1} $ by pencil and paper.\n",
    "\n",
    "In the exercise below, you will be asked to solve these equations numerically."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecc7a220",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68debadd",
   "metadata": {},
   "source": [
    "## Exercise 26.1\n",
    "\n",
    "Solve for the dynamics of equilibrium capital stock in the CRRA case numerically using [(26.22)](#equation-law-of-motion-capital-crra).\n",
    "\n",
    "Visualize the dynamics using a 45-degree diagram."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc73f831",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 26.1](https://intro.quantecon.org/#olg_ex1)\n",
    "\n",
    "To solve for $ k_{t+1} $ given $ k_t $ we use [Newton’s method](https://python.quantecon.org/newton_method.html).\n",
    "\n",
    "Let\n",
    "\n",
    "\n",
    "<a id='equation-crra-newton-1'></a>\n",
    "$$\n",
    "f(k_{t+1}, k_t)\n",
    "    =\n",
    "    k_{t+1} \n",
    "    \\left[ \n",
    "        1 + \\beta^{-1/\\gamma} \n",
    "        \\left ( \n",
    "            \\alpha k^{\\alpha-1}_{t+1} \n",
    "        \\right )^{(\\gamma-1)/\\gamma} \n",
    "    \\right] - (1-\\alpha) k^{\\alpha}_t =0 \\tag{26.23}\n",
    "$$\n",
    "\n",
    "If $ k_t $ is given then $ f $ is a function of unknown $ k_{t+1} $.\n",
    "\n",
    "Then we can use `scipy.optimize.newton` to solve $ f(k_{t+1}, k_t)=0 $ for $ k_{t+1} $.\n",
    "\n",
    "First let’s define $ f $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c65f5e3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def f(k_prime, k, model):\n",
    "    α, β, γ = model.α, model.β, model.γ\n",
    "    z = (1 - α) * k**α\n",
    "    a = α**(1-1/γ)\n",
    "    b = k_prime**((α * γ - α + 1) / γ)\n",
    "    p = k_prime + k_prime * β**(-1/γ) * a * b\n",
    "    return p - z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e394582",
   "metadata": {},
   "source": [
    "Now let’s define a function that finds the value of $ k_{t+1} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbe8975c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def k_update(k, model):\n",
    "    return optimize.newton(lambda k_prime: f(k_prime, k, model), 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "909be9b4",
   "metadata": {},
   "source": [
    "Finally, here is the 45-degree diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "948f53e8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "kmin, kmax = 0, 0.5\n",
    "n = 1000\n",
    "k_grid = np.linspace(kmin, kmax, n)\n",
    "k_grid_next = np.empty_like(k_grid)\n",
    "\n",
    "for i in range(n):\n",
    "    k_grid_next[i] = k_update(k_grid[i], model)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "ymin, ymax = np.min(k_grid_next), np.max(k_grid_next)\n",
    "\n",
    "ax.plot(k_grid, k_grid_next,  lw=2, alpha=0.6, label='$g$')\n",
    "ax.plot(k_grid, k_grid, 'k-', lw=1, alpha=0.7, label='$45^{\\circ}$')\n",
    "\n",
    "\n",
    "ax.legend(loc='upper left', frameon=False, fontsize=12)\n",
    "ax.set_xlabel('$k_t$', fontsize=12)\n",
    "ax.set_ylabel('$k_{t+1}$', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8227dadc",
   "metadata": {},
   "source": [
    "## Exercise 26.2\n",
    "\n",
    "The 45-degree diagram from the last exercise shows that there is a unique\n",
    "positive steady state.\n",
    "\n",
    "The positive steady state can be obtained by setting  $ k_{t+1} = k_t = k^* $ in [(26.22)](#equation-law-of-motion-capital-crra), which yields\n",
    "\n",
    "$$\n",
    "k^* = \n",
    "    \\frac{(1-\\alpha)(k^*)^{\\alpha}}\n",
    "    {1 + \\beta^{-1/\\gamma} (\\alpha (k^*)^{\\alpha-1})^{(\\gamma-1)/\\gamma}}\n",
    "$$\n",
    "\n",
    "Unlike the log preference case, the CRRA utility steady state $ k^* $\n",
    "cannot be obtained analytically.\n",
    "\n",
    "Instead, we solve for $ k^* $ using Newton’s method."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91c52134",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 26.2](https://intro.quantecon.org/#olg_ex2)\n",
    "\n",
    "We introduce a function $ h $ such that\n",
    "positive steady state is the root of $ h $.\n",
    "\n",
    "\n",
    "<a id='equation-crra-newton-2'></a>\n",
    "$$\n",
    "h(k^*) = k^*  \n",
    "    \\left [ \n",
    "        1 + \\beta^{-1/\\gamma} (\\alpha (k^*)^{\\alpha-1})^{(\\gamma-1)/\\gamma} \n",
    "    \\right ] - (1-\\alpha)(k^*)^{\\alpha} \\tag{26.24}\n",
    "$$\n",
    "\n",
    "Here it is in Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75b020ac",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def h(k_star, model):\n",
    "    α, β, γ = model.α, model.β, model.γ\n",
    "    z = (1 - α) * k_star**α\n",
    "    R1 = α ** (1-1/γ)\n",
    "    R2 = k_star**((α * γ - α + 1) / γ)\n",
    "    p = k_star + k_star * β**(-1/γ) * R1 * R2\n",
    "    return p - z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc8b09ff",
   "metadata": {},
   "source": [
    "Let’s apply Newton’s method to find the root:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff592d74",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "k_star = optimize.newton(h, 0.2, args=(model,))\n",
    "print(f\"k_star = {k_star}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06e0d5c1",
   "metadata": {},
   "source": [
    "## Exercise 26.3\n",
    "\n",
    "Generate three time paths for capital, from\n",
    "three distinct initial conditions, under the parameterization listed above.\n",
    "\n",
    "Use initial conditions for $ k_0 $ of $ 0.001, 1.2, 2.6 $ and time series length 10."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15a04654",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 26.3](https://intro.quantecon.org/#olg_ex3)\n",
    "\n",
    "Let’s define the constants and three distinct intital conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e642f9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ts_length = 10\n",
    "k0 = np.array([0.001, 1.2, 2.6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511f106f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def simulate_ts(model, k0_values, ts_length):\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ts = np.zeros(ts_length)\n",
    "\n",
    "    # simulate and plot time series\n",
    "    for k_init in k0_values:\n",
    "        ts[0] = k_init\n",
    "        for t in range(1, ts_length):\n",
    "            ts[t] = k_update(ts[t-1], model)\n",
    "        ax.plot(np.arange(ts_length), ts, '-o', ms=4, alpha=0.6,\n",
    "                label=r'$k_0=%g$' %k_init)\n",
    "    ax.plot(np.arange(ts_length), np.full(ts_length, k_star),\n",
    "            alpha=0.6, color='red', label=r'$k^*$')\n",
    "    ax.legend(fontsize=10)\n",
    "\n",
    "    ax.set_xlabel(r'$t$', fontsize=14)\n",
    "    ax.set_ylabel(r'$k_t$', fontsize=14)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c612d65c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "simulate_ts(model, k0, ts_length)"
   ]
  }
 ],
 "metadata": {
  "date": 1736815200.213997,
  "filename": "olg.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "The Overlapping Generations Model"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}