cass_koopmans_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "610e7993",
   "metadata": {},
   "source": [
    "\n",
    "<a id='cass-koopmans-1'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "690339a9",
   "metadata": {},
   "source": [
    "# Cass-Koopmans Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84cd315d",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture and [Cass-Koopmans Competitive Equilibrium](https://python.quantecon.org/cass_koopmans_2.html) describe a model that Tjalling Koopmans [[Koopmans, 1965](https://python.quantecon.org/zreferences.html#id91)]\n",
    "and David Cass [[Cass, 1965](https://python.quantecon.org/zreferences.html#id92)] used to analyze optimal growth.\n",
    "\n",
    "The model extends  the model of Robert Solow\n",
    "described in [an earlier lecture](https://python-programming.quantecon.org/python_oop.html).\n",
    "\n",
    "It does so by making  saving rate be a  decision, instead of a hard-wired constant.\n",
    "\n",
    "(Solow assumed a constant saving rate determined outside the model.)\n",
    "\n",
    "We describe two versions of the model, a planning problem  in this lecture, and a competitive equilibrium   in this lecture  [Cass-Koopmans Competitive Equilibrium](https://python.quantecon.org/cass_koopmans_2.html).\n",
    "\n",
    "Together, the two lectures  illustrate what is, in fact, a\n",
    "more general connection between a **planned economy** and a decentralized economy\n",
    "organized as a **competitive equilibrium**.\n",
    "\n",
    "This lecture is devoted to the planned economy version.\n",
    "\n",
    "In the planned economy, there are\n",
    "\n",
    "- no prices  \n",
    "- no budget constraints  \n",
    "\n",
    "\n",
    "Instead there is a dictator that tells people\n",
    "\n",
    "- what to produce  \n",
    "- what to invest in physical capital  \n",
    "- who is to consume what  and when  \n",
    "\n",
    "\n",
    "The lecture uses important ideas including\n",
    "\n",
    "- A min-max problem for solving a planning problem.  \n",
    "- A **shooting algorithm** for solving difference equations subject\n",
    "  to initial and terminal conditions.  \n",
    "- A **turnpike** property of  optimal paths for\n",
    "  long but finite-horizon economies.  \n",
    "- A **stable manifold** and a **phase plane**  \n",
    "\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d144b26",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from numba import jit, float64\n",
    "from numba.experimental import jitclass\n",
    "import numpy as np\n",
    "from quantecon.optimize import brentq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85f07af0",
   "metadata": {},
   "source": [
    "## The Model\n",
    "\n",
    "Time is discrete and takes values $ t = 0, 1 , \\ldots, T $ where $ T $ is  finite.\n",
    "\n",
    "(We’ll eventually study a limiting case in which  $ T = + \\infty $)\n",
    "\n",
    "A single good can either be consumed or invested in physical capital.\n",
    "\n",
    "The consumption good is not durable and depreciates completely if not\n",
    "consumed immediately.\n",
    "\n",
    "The capital good is durable but depreciates.\n",
    "\n",
    "We let $ C_t $ be the total consumption of a nondurable  consumption good at time $ t $.\n",
    "\n",
    "Let $ K_t $ be the stock of physical capital at time $ t $.\n",
    "\n",
    "Let $ \\vec{C} $ = $ \\{C_0,\\dots, C_T\\} $ and\n",
    "$ \\vec{K} $ = $ \\{K_0,\\dots,K_{T+1}\\} $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40ccc024",
   "metadata": {},
   "source": [
    "### Digression: Aggregation Theory\n",
    "\n",
    "We use a concept of a representative consumer to be thought of as follows.\n",
    "\n",
    "There is a unit mass of identical consumers indexed by $ \\omega \\in [0,1] $.\n",
    "\n",
    "Consumption of consumer $ \\omega $ is $ c(\\omega) $.\n",
    "\n",
    "Aggregate consumption is\n",
    "\n",
    "$$\n",
    "C = \\int_0^1 c(\\omega) d \\omega\n",
    "$$\n",
    "\n",
    "Consider  a welfare problem that chooses an allocation $ \\{c(\\omega)\\} $ across consumers to maximize\n",
    "\n",
    "$$\n",
    "\\int_0^1 u(c(\\omega)) d \\omega\n",
    "$$\n",
    "\n",
    "where $ u(\\cdot) $ is a concave utility function with $ u' >0, u'' < 0 $ and  maximization is subject to\n",
    "\n",
    "\n",
    "<a id='equation-eq-feas200'></a>\n",
    "$$\n",
    "C = \\int_0^1 c(\\omega) d \\omega . \\tag{34.1}\n",
    "$$\n",
    "\n",
    "Form a Lagrangian $ L = \\int_0^1 u(c(\\omega)) d \\omega + \\lambda [C - \\int_0^1 c(\\omega) d \\omega ] $.\n",
    "\n",
    "Differentiate under the integral signs with respect to each $ \\omega $ to  obtain the first-order\n",
    "necessary conditions\n",
    "\n",
    "$$\n",
    "u'(c(\\omega)) = \\lambda.\n",
    "$$\n",
    "\n",
    "These conditions imply that $ c(\\omega) $ equals a constant $ c $ that is independent\n",
    "of $ \\omega $.\n",
    "\n",
    "To find $ c $, use feasibility constraint [(34.1)](#equation-eq-feas200) to conclude that\n",
    "\n",
    "$$\n",
    "c(\\omega) = c = C.\n",
    "$$\n",
    "\n",
    "This line of argument indicates the special *aggregation theory* that lies beneath outcomes in which a representative consumer\n",
    "consumes amount $ C $.\n",
    "\n",
    "It appears often in aggregate economics.\n",
    "\n",
    "We shall use this aggregation theory here and also  in  this lecture  [Cass-Koopmans Competitive Equilibrium](https://python.quantecon.org/cass_koopmans_2.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cdb8925",
   "metadata": {},
   "source": [
    "#### An  Economy\n",
    "\n",
    "A representative household is endowed with one unit of labor at each\n",
    "$ t $ and likes the consumption good at each $ t $.\n",
    "\n",
    "The representative household inelastically supplies a single unit of\n",
    "labor $ N_t $ at each $ t $, so that\n",
    "$ N_t =1 \\text{ for all } t \\in \\{0, 1, \\ldots,  T\\} $.\n",
    "\n",
    "The representative household has preferences over consumption bundles\n",
    "ordered by the utility functional:\n",
    "\n",
    "\n",
    "<a id='equation-utility-functional'></a>\n",
    "$$\n",
    "U(\\vec{C}) = \\sum_{t=0}^{T} \\beta^t \\frac{C_t^{1-\\gamma}}{1-\\gamma} \\tag{34.2}\n",
    "$$\n",
    "\n",
    "where $ \\beta \\in (0,1) $ is a discount factor and $ \\gamma >0 $\n",
    "governs the curvature of the one-period utility function.\n",
    "\n",
    "Larger $ \\gamma $’s imply more curvature.\n",
    "\n",
    "Note that\n",
    "\n",
    "\n",
    "<a id='equation-utility-oneperiod'></a>\n",
    "$$\n",
    "u(C_t) = \\frac{C_t^{1-\\gamma}}{1-\\gamma} \\tag{34.3}\n",
    "$$\n",
    "\n",
    "satisfies $ u'>0,u''<0 $.\n",
    "\n",
    "$ u' > 0 $ asserts that the consumer prefers more to less.\n",
    "\n",
    "$ u''< 0 $ asserts that marginal utility declines with increases\n",
    "in $ C_t $.\n",
    "\n",
    "We assume that $ K_0 > 0 $ is an  exogenous  initial\n",
    "capital stock.\n",
    "\n",
    "There is an economy-wide production function\n",
    "\n",
    "\n",
    "<a id='equation-production-function'></a>\n",
    "$$\n",
    "F(K_t,N_t) = A K_t^{\\alpha}N_t^{1-\\alpha} \\tag{34.4}\n",
    "$$\n",
    "\n",
    "with $ 0 < \\alpha<1 $, $ A > 0 $.\n",
    "\n",
    "A feasible allocation $ \\vec{C}, \\vec{K} $ satisfies\n",
    "\n",
    "\n",
    "<a id='equation-allocation'></a>\n",
    "$$\n",
    "C_t + K_{t+1} \\leq F(K_t,N_t) + (1-\\delta) K_t \\quad \\text{for all } t \\in \\{0, 1, \\ldots,  T\\} \\tag{34.5}\n",
    "$$\n",
    "\n",
    "where $ \\delta \\in (0,1) $ is a depreciation rate of capital."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2092f86d",
   "metadata": {},
   "source": [
    "## Planning Problem\n",
    "\n",
    "A planner chooses an allocation $ \\{\\vec{C},\\vec{K}\\} $ to\n",
    "maximize [(34.2)](#equation-utility-functional) subject to [(34.5)](#equation-allocation).\n",
    "\n",
    "Let $ \\vec{\\mu}=\\{\\mu_0,\\dots,\\mu_T\\} $ be a sequence of\n",
    "nonnegative **Lagrange multipliers**.\n",
    "\n",
    "To find an optimal allocation, form a Lagrangian\n",
    "\n",
    "\n",
    "<a id='equation-eq-lagrangian201'></a>\n",
    "$$\n",
    "\\mathcal{L}(\\vec{C} ,\\vec{K} ,\\vec{\\mu} ) =\n",
    "\\sum_{t=0}^T \\beta^t\\left\\{ u(C_t)+ \\mu_t\n",
    "\\left(F(K_t,1) + (1-\\delta) K_t- C_t - K_{t+1} \\right)\\right\\} \\tag{34.6}\n",
    "$$\n",
    "\n",
    "and  pose the following min-max problem:\n",
    "\n",
    "\n",
    "<a id='equation-min-max-prob'></a>\n",
    "$$\n",
    "\\min_{\\vec{\\mu}} \\max_{\\vec{C},\\vec{K}} \\mathcal{L}(\\vec{C},\\vec{K},\\vec{\\mu} ) \\tag{34.7}\n",
    "$$\n",
    "\n",
    "- **Extremization** means\n",
    "  maximization with respect to $ \\vec{C}, \\vec{K} $ and\n",
    "  minimization with respect to $ \\vec{\\mu} $.  \n",
    "- Our problem satisfies\n",
    "  conditions that assure that  second-order\n",
    "  conditions are satisfied at an allocation that satisfies the\n",
    "  first-order necessary conditions that we are about to compute.  \n",
    "\n",
    "\n",
    "Before computing first-order conditions, we present some handy formulas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a849a88",
   "metadata": {},
   "source": [
    "### Useful Properties of Linearly Homogeneous Production Function\n",
    "\n",
    "The following technicalities will help us.\n",
    "\n",
    "Notice that\n",
    "\n",
    "$$\n",
    "F(K_t,N_t) = A K_t^\\alpha N_t^{1-\\alpha} = N_t A\\left(\\frac{K_t}{N_t}\\right)^\\alpha\n",
    "$$\n",
    "\n",
    "Define the **output per-capita production function**\n",
    "\n",
    "$$\n",
    "\\frac{F(K_t,N_t)}{N_t} \\equiv f\\left(\\frac{K_t}{N_t}\\right) = A\\left(\\frac{K_t}{N_t}\\right)^\\alpha\n",
    "$$\n",
    "\n",
    "whose argument is **capital per-capita**.\n",
    "\n",
    "It is useful to recall the following calculations for the marginal product of capital\n",
    "\n",
    "\n",
    "<a id='equation-useful-calc1'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial F(K_t,N_t)}{\\partial K_t}\n",
    "& =\n",
    "\\frac{\\partial N_t f\\left( \\frac{K_t}{N_t}\\right)}{\\partial K_t}\n",
    "\\\\ &=\n",
    "N_t f'\\left(\\frac{K_t}{N_t}\\right)\\frac{1}{N_t} \\quad \\text{(Chain rule)}\n",
    "\\\\ &=\n",
    "f'\\left.\\left(\\frac{K_t}{N_t}\\right)\\right|_{N_t=1}\n",
    "\\\\ &= f'(K_t)\n",
    "\\end{aligned} \\tag{34.8}\n",
    "$$\n",
    "\n",
    "and the marginal product of labor\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial F(K_t,N_t)}{\\partial N_t}\n",
    "&=\n",
    "\\frac{\\partial N_t f\\left( \\frac{K_t}{N_t}\\right)}{\\partial N_t} \\quad \\text{(Product rule)}\n",
    "\\\\ &=\n",
    "f\\left(\\frac{K_t}{N_t}\\right){+} N_t f'\\left(\\frac{K_t}{N_t}\\right) \\frac{-K_t}{N_t^2} \\quad \\text{(Chain rule)}\n",
    "\\\\ &=\n",
    "f\\left(\\frac{K_t}{N_t}\\right){-}\\frac{K_t}{N_t}f'\\left.\\left(\\frac{K_t}{N_t}\\right)\\right|_{N_t=1}\n",
    "\\\\ &=\n",
    "f(K_t) - f'(K_t) K_t\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "(Here we are using that $ N_t = 1 $ for all $ t $, so that $ K_t = \\frac{K_t}{N_t} $.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "667656ba",
   "metadata": {},
   "source": [
    "### First-order necessary conditions\n",
    "\n",
    "We now compute **first-order necessary conditions** for extremization of  Lagrangian [(34.6)](#equation-eq-lagrangian201):\n",
    "\n",
    "\n",
    "<a id='equation-constraint1'></a>\n",
    "$$\n",
    "C_t: \\qquad u'(C_t)-\\mu_t=0 \\qquad \\text{for all} \\quad t= 0,1,\\dots,T \\tag{34.9}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-constraint2'></a>\n",
    "$$\n",
    "K_t: \\qquad \\beta \\mu_t\\left[(1-\\delta)+f'(K_t)\\right] - \\mu_{t-1}=0 \\qquad \\text{for all } \\quad t=1,2,\\dots,T \\tag{34.10}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-constraint3'></a>\n",
    "$$\n",
    "\\mu_t:\\qquad F(K_t,1)+ (1-\\delta) K_t  - C_t - K_{t+1}=0 \\qquad \\text{for all } \\quad t=0,1,\\dots,T \\tag{34.11}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-constraint4'></a>\n",
    "$$\n",
    "K_{T+1}: \\qquad -\\mu_T \\leq 0, \\ \\leq 0 \\text{ if } K_{T+1}=0; \\ =0 \\text{ if } K_{T+1}>0 \\tag{34.12}\n",
    "$$\n",
    "\n",
    "In computing  [(34.10)](#equation-constraint2) we recognize that $ K_t $ appears\n",
    "in both the time  $ t $ and time $ t-1 $ feasibility constraints [(34.5)](#equation-allocation).\n",
    "\n",
    "Restrictions [(34.12)](#equation-constraint4) come from differentiating with respect\n",
    "to $ K_{T+1} $ and applying the following **Karush-Kuhn-Tucker condition** (KKT)\n",
    "(see [Karush-Kuhn-Tucker conditions](https://en.wikipedia.org/wiki/Karush-Kuhn-Tucker_conditions)):\n",
    "\n",
    "\n",
    "<a id='equation-kkt'></a>\n",
    "$$\n",
    "\\mu_T K_{T+1}=0 \\tag{34.13}\n",
    "$$\n",
    "\n",
    "Combining [(34.9)](#equation-constraint1) and [(34.10)](#equation-constraint2) gives\n",
    "\n",
    "$$\n",
    "\\beta u'\\left(C_t\\right)\\left[(1-\\delta)+f'\\left(K_t\\right)\\right]-u'\\left(C_{t-1}\\right)=0\n",
    "\\quad \\text{ for all } t=1,2,\\dots, T+1\n",
    "$$\n",
    "\n",
    "which can be rearranged to become\n",
    "\n",
    "\n",
    "<a id='equation-l12'></a>\n",
    "$$\n",
    "\\beta u'\\left(C_{t+1}\\right)\\left[(1-\\delta)+f'\\left(K_{t+1}\\right)\\right]=\n",
    "u'\\left(C_{t}\\right) \\quad \\text{ for all } t=0,1,\\dots, T \\tag{34.14}\n",
    "$$\n",
    "\n",
    "Applying  the inverse  marginal utility  of consumption  function on both sides of the above\n",
    "equation gives\n",
    "\n",
    "$$\n",
    "C_{t+1} =u'^{-1}\\left(\\left(\\frac{\\beta}{u'(C_t)}[f'(K_{t+1}) +(1-\\delta)]\\right)^{-1}\\right)\n",
    "$$\n",
    "\n",
    "which for our utility function [(34.3)](#equation-utility-oneperiod) becomes the consumption **Euler\n",
    "equation**\n",
    "\n",
    "\n",
    "<a id='equation-eq-consn-euler'></a>\n",
    "$$\n",
    "\\begin{aligned} C_{t+1} =\\left(\\beta C_t^{\\gamma}[f'(K_{t+1}) +\n",
    "(1-\\delta)]\\right)^{1/\\gamma} \n",
    "%\\notag\\\\= C_t\\left(\\beta [f'(K_{t+1}) +\n",
    "%(1-\\delta)]\\right)^{1/\\gamma} \n",
    "\\end{aligned} \\tag{34.15}\n",
    "$$\n",
    "\n",
    "which we can  combine with the feasibility constraint [(34.5)](#equation-allocation) to get\n",
    "\n",
    "\n",
    "<a id='equation-eq-systemdynamics'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "C_{t+1} & = C_t\\left(\\beta [f'(F(K_t,1)+ (1-\\delta) K_t  - C_t) +\n",
    "(1-\\delta)]\\right)^{1/\\gamma}  \\\\\n",
    "K_{t+1}  & = F(K_t,1)+ (1-\\delta) K_t  - C_t .\n",
    "\\end{aligned} \\tag{34.16}\n",
    "$$\n",
    "\n",
    "This is a pair of non-linear first-order difference equations that map $ C_t, K_t $ into $ C_{t+1}, K_{t+1} $ and that  an optimal sequence $ \\vec C , \\vec K $ must satisfy.\n",
    "\n",
    "It must also satisfy the initial condition that $ K_0 $ is given and $ K_{T+1} = 0 $.\n",
    "\n",
    "Below we define a `jitclass` that stores parameters and functions\n",
    "that define our economy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f4baec5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "planning_data = [\n",
    "    ('γ', float64),    # Coefficient of relative risk aversion\n",
    "    ('β', float64),    # Discount factor\n",
    "    ('δ', float64),    # Depreciation rate on capital\n",
    "    ('α', float64),    # Return to capital per capita\n",
    "    ('A', float64)     # Technology\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "435013e6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jitclass(planning_data)\n",
    "class PlanningProblem():\n",
    "\n",
    "    def __init__(self, γ=2, β=0.95, δ=0.02, α=0.33, A=1):\n",
    "\n",
    "        self.γ, self.β = γ, β\n",
    "        self.δ, self.α, self.A = δ, α, A\n",
    "\n",
    "    def u(self, c):\n",
    "        '''\n",
    "        Utility function\n",
    "        ASIDE: If you have a utility function that is hard to solve by hand\n",
    "        you can use automatic or symbolic differentiation\n",
    "        See https://github.com/HIPS/autograd\n",
    "        '''\n",
    "        γ = self.γ\n",
    "\n",
    "        return c ** (1 - γ) / (1 - γ) if γ!= 1 else np.log(c)\n",
    "\n",
    "    def u_prime(self, c):\n",
    "        'Derivative of utility'\n",
    "        γ = self.γ\n",
    "\n",
    "        return c ** (-γ)\n",
    "\n",
    "    def u_prime_inv(self, c):\n",
    "        'Inverse of derivative of utility'\n",
    "        γ = self.γ\n",
    "\n",
    "        return c ** (-1 / γ)\n",
    "\n",
    "    def f(self, k):\n",
    "        'Production function'\n",
    "        α, A = self.α, self.A\n",
    "\n",
    "        return A * k ** α\n",
    "\n",
    "    def f_prime(self, k):\n",
    "        'Derivative of production function'\n",
    "        α, A = self.α, self.A\n",
    "\n",
    "        return α * A * k ** (α - 1)\n",
    "\n",
    "    def f_prime_inv(self, k):\n",
    "        'Inverse of derivative of production function'\n",
    "        α, A = self.α, self.A\n",
    "\n",
    "        return (k / (A * α)) ** (1 / (α - 1))\n",
    "\n",
    "    def next_k_c(self, k, c):\n",
    "        ''''\n",
    "        Given the current capital Kt and an arbitrary feasible\n",
    "        consumption choice Ct, computes Kt+1 by state transition law\n",
    "        and optimal Ct+1 by Euler equation.\n",
    "        '''\n",
    "        β, δ = self.β, self.δ\n",
    "        u_prime, u_prime_inv = self.u_prime, self.u_prime_inv\n",
    "        f, f_prime = self.f, self.f_prime\n",
    "\n",
    "        k_next = f(k) + (1 - δ) * k - c\n",
    "        c_next = u_prime_inv(u_prime(c) / (β * (f_prime(k_next) + (1 - δ))))\n",
    "\n",
    "        return k_next, c_next"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07a052b6",
   "metadata": {},
   "source": [
    "We can construct an economy with the Python code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa80db21",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "pp = PlanningProblem()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9538c35a",
   "metadata": {},
   "source": [
    "## Shooting Algorithm\n",
    "\n",
    "We use  **shooting** to compute an optimal allocation\n",
    "$ \\vec{C}, \\vec{K} $ and an associated Lagrange multiplier sequence\n",
    "$ \\vec{\\mu} $.\n",
    "\n",
    "First-order necessary conditions\n",
    "[(34.9)](#equation-constraint1), [(34.10)](#equation-constraint2), and\n",
    "[(34.11)](#equation-constraint3)  for the planning problem form a system of **difference equations** with\n",
    "two boundary conditions:\n",
    "\n",
    "- $ K_0 $ is a given **initial condition** for capital  \n",
    "- $ K_{T+1} =0 $ is a **terminal condition** for capital that we\n",
    "  deduced from the first-order necessary condition for $ K_{T+1} $\n",
    "  the KKT condition [(34.13)](#equation-kkt)  \n",
    "\n",
    "\n",
    "We have no initial condition for the Lagrange multiplier\n",
    "$ \\mu_0 $.\n",
    "\n",
    "If we did, our job would be easy:\n",
    "\n",
    "- Given $ \\mu_0 $ and $ k_0 $, we could compute $ c_0 $ from\n",
    "  equation [(34.9)](#equation-constraint1) and then $ k_1 $ from equation\n",
    "  [(34.11)](#equation-constraint3) and $ \\mu_1 $ from equation\n",
    "  [(34.10)](#equation-constraint2).  \n",
    "- We could continue in this way to compute the remaining elements of\n",
    "  $ \\vec{C}, \\vec{K}, \\vec{\\mu} $.  \n",
    "\n",
    "\n",
    "However, we woujld not be assured that the Kuhn-Tucker condition [(34.13)](#equation-kkt) would be satisfied.\n",
    "\n",
    "Furthermore,  we don’t have an initial condition for $ \\mu_0 $.\n",
    "\n",
    "So this won’t work.\n",
    "\n",
    "Indeed, part of our task is to compute the **optimal** value of $ \\mu_0 $.\n",
    "\n",
    "To compute $ \\mu_0 $ and the other objects we want, a  simple modification of the above procedure will work.\n",
    "\n",
    "It is  called the **shooting algorithm**.\n",
    "\n",
    "It is  an instance of a **guess and verify**\n",
    "algorithm that consists of the following steps:\n",
    "\n",
    "- Guess an initial Lagrange multiplier $ \\mu_0 $.  \n",
    "- Apply the **simple algorithm** described above.  \n",
    "- Compute $ K_{T+1} $ and check whether it\n",
    "  equals zero.  \n",
    "- If  $ K_{T+1} =0 $, we have solved the problem.  \n",
    "- If $ K_{T+1} > 0 $, lower $ \\mu_0 $ and try again.  \n",
    "- If $ K_{T+1} < 0 $, raise $ \\mu_0 $ and try again.  \n",
    "\n",
    "\n",
    "The following Python code implements the shooting algorithm for the\n",
    "planning problem.\n",
    "\n",
    "(Actually, we modified the preceding  algorithm slightly by starting with a guess for\n",
    "$ c_0 $ instead of $ \\mu_0 $ in the following code.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d887224",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def shooting(pp, c0, k0, T=10):\n",
    "    '''\n",
    "    Given the initial condition of capital k0 and an initial guess\n",
    "    of consumption c0, computes the whole paths of c and k\n",
    "    using the state transition law and Euler equation for T periods.\n",
    "    '''\n",
    "    if c0 > pp.f(k0) + (1 - pp.δ) * k0:\n",
    "        print(\"initial consumption is not feasible\")\n",
    "\n",
    "        return None\n",
    "\n",
    "    # initialize vectors of c and k\n",
    "    c_vec = np.empty(T+1)\n",
    "    k_vec = np.empty(T+2)\n",
    "\n",
    "    c_vec[0] = c0\n",
    "    k_vec[0] = k0\n",
    "\n",
    "    for t in range(T):\n",
    "        k_vec[t+1], c_vec[t+1] = pp.next_k_c(k_vec[t], c_vec[t])\n",
    "\n",
    "    k_vec[T+1] = pp.f(k_vec[T]) + (1 - pp.δ) * k_vec[T] - c_vec[T]\n",
    "\n",
    "    return c_vec, k_vec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "391f8f3e",
   "metadata": {},
   "source": [
    "We’ll start with an incorrect guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd37e38b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "paths = shooting(pp, 0.2, 0.3, T=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ebbbff",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(14, 5))\n",
    "\n",
    "colors = ['blue', 'red']\n",
    "titles = ['Consumption', 'Capital']\n",
    "ylabels = ['$c_t$', '$k_t$']\n",
    "\n",
    "T = paths[0].size - 1\n",
    "for i in range(2):\n",
    "    axs[i].plot(paths[i], c=colors[i])\n",
    "    axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])\n",
    "\n",
    "axs[1].scatter(T+1, 0, s=80)\n",
    "axs[1].axvline(T+1, color='k', ls='--', lw=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35de7cf0",
   "metadata": {},
   "source": [
    "Evidently, our initial guess for $ \\mu_0 $ is too high, so\n",
    "initial consumption too low.\n",
    "\n",
    "We know this because we miss our $ K_{T+1}=0 $ target on the high\n",
    "side.\n",
    "\n",
    "Now we automate things with a search-for-a-good $ \\mu_0 $\n",
    "algorithm that stops when we hit the target $ K_{t+1} = 0 $.\n",
    "\n",
    "We use a **bisection method**.\n",
    "\n",
    "We make an initial guess for $ C_0 $ (we can eliminate\n",
    "$ \\mu_0 $ because $ C_0 $ is an exact function of\n",
    "$ \\mu_0 $).\n",
    "\n",
    "We know that the lowest $ C_0 $ can ever be is $ 0 $ and that the\n",
    "largest it can be is initial output $ f(K_0) $.\n",
    "\n",
    "Guess $ C_0 $  and shoot forward to $ T+1 $.\n",
    "\n",
    "If $ K_{T+1}>0 $, we take it to be our new **lower** bound\n",
    "on $ C_0 $.\n",
    "\n",
    "If $ K_{T+1}<0 $, we take  it to be our new **upper** bound.\n",
    "\n",
    "Make a new guess for $ C_0 $ that is  halfway between our new\n",
    "upper and lower bounds.\n",
    "\n",
    "Shoot forward again, iterating on these steps until we converge.\n",
    "\n",
    "When $ K_{T+1} $ gets close enough to $ 0 $ (i.e., within an error\n",
    "tolerance bounds), we stop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da73a675",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def bisection(pp, c0, k0, T=10, tol=1e-4, max_iter=500, k_ter=0, verbose=True):\n",
    "\n",
    "    # initial boundaries for guess c0\n",
    "    c0_upper = pp.f(k0)\n",
    "    c0_lower = 0\n",
    "\n",
    "    i = 0\n",
    "    while True:\n",
    "        c_vec, k_vec = shooting(pp, c0, k0, T)\n",
    "        error = k_vec[-1] - k_ter\n",
    "\n",
    "        # check if the terminal condition is satisfied\n",
    "        if np.abs(error) < tol:\n",
    "            if verbose:\n",
    "                print('Converged successfully on iteration ', i+1)\n",
    "            return c_vec, k_vec\n",
    "\n",
    "        i += 1\n",
    "        if i == max_iter:\n",
    "            if verbose:\n",
    "                print('Convergence failed.')\n",
    "            return c_vec, k_vec\n",
    "\n",
    "        # if iteration continues, updates boundaries and guess of c0\n",
    "        if error > 0:\n",
    "            c0_lower = c0\n",
    "        else:\n",
    "            c0_upper = c0\n",
    "\n",
    "        c0 = (c0_lower + c0_upper) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9b40c4f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_paths(pp, c0, k0, T_arr, k_ter=0, k_ss=None, axs=None):\n",
    "\n",
    "    if axs is None:\n",
    "        fix, axs = plt.subplots(1, 3, figsize=(16, 4))\n",
    "    ylabels = ['$c_t$', '$k_t$', '$\\mu_t$']\n",
    "    titles = ['Consumption', 'Capital', 'Lagrange Multiplier']\n",
    "\n",
    "    c_paths = []\n",
    "    k_paths = []\n",
    "    for T in T_arr:\n",
    "        c_vec, k_vec = bisection(pp, c0, k0, T, k_ter=k_ter, verbose=False)\n",
    "        c_paths.append(c_vec)\n",
    "        k_paths.append(k_vec)\n",
    "\n",
    "        μ_vec = pp.u_prime(c_vec)\n",
    "        paths = [c_vec, k_vec, μ_vec]\n",
    "\n",
    "        for i in range(3):\n",
    "            axs[i].plot(paths[i])\n",
    "            axs[i].set(xlabel='t', ylabel=ylabels[i], title=titles[i])\n",
    "\n",
    "        # Plot steady state value of capital\n",
    "        if k_ss is not None:\n",
    "            axs[1].axhline(k_ss, c='k', ls='--', lw=1)\n",
    "\n",
    "        axs[1].axvline(T+1, c='k', ls='--', lw=1)\n",
    "        axs[1].scatter(T+1, paths[1][-1], s=80)\n",
    "\n",
    "    return c_paths, k_paths"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db384bda",
   "metadata": {},
   "source": [
    "Now we can solve the model and plot the paths of consumption, capital, and Lagrange multiplier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc35705c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_paths(pp, 0.3, 0.3, [10]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff0bb21d",
   "metadata": {},
   "source": [
    "## Setting Initial Capital to Steady State Capital\n",
    "\n",
    "When  $ T \\rightarrow +\\infty $, the optimal allocation converges to\n",
    "steady state values of $ C_t $ and $ K_t $.\n",
    "\n",
    "It is instructive to set $ K_0 $ equal\n",
    "to the $ \\lim_{T \\rightarrow + \\infty } K_t $, which we’ll call  steady state capital.\n",
    "\n",
    "In a steady state $ K_{t+1} = K_t=\\bar{K} $ for all very\n",
    "large $ t $.\n",
    "\n",
    "Evalauating  feasibility constraint [(34.5)](#equation-allocation) at $ \\bar K $ gives\n",
    "\n",
    "\n",
    "<a id='equation-feasibility-constraint'></a>\n",
    "$$\n",
    "f(\\bar{K})-\\delta \\bar{K} = \\bar{C} \\tag{34.17}\n",
    "$$\n",
    "\n",
    "Substituting $ K_t = \\bar K $ and $ C_t=\\bar C $ for\n",
    "all $ t $ into [(34.14)](#equation-l12) gives\n",
    "\n",
    "$$\n",
    "1=\\beta \\frac{u'(\\bar{C})}{u'(\\bar{C})}[f'(\\bar{K})+(1-\\delta)]\n",
    "$$\n",
    "\n",
    "Defining $ \\beta = \\frac{1}{1+\\rho} $, and cancelling gives\n",
    "\n",
    "$$\n",
    "1+\\rho = 1[f'(\\bar{K}) + (1-\\delta)]\n",
    "$$\n",
    "\n",
    "Simplifying gives\n",
    "\n",
    "$$\n",
    "f'(\\bar{K}) = \\rho +\\delta\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\bar{K} = f'^{-1}(\\rho+\\delta)\n",
    "$$\n",
    "\n",
    "For   production function [(34.4)](#equation-production-function), this becomes\n",
    "\n",
    "$$\n",
    "\\alpha \\bar{K}^{\\alpha-1} = \\rho + \\delta\n",
    "$$\n",
    "\n",
    "As an example, after setting $ \\alpha= .33 $,\n",
    "$ \\rho = 1/\\beta-1 =1/(19/20)-1 = 20/19-19/19 = 1/19 $, $ \\delta = 1/50 $,\n",
    "we get\n",
    "\n",
    "$$\n",
    "\\bar{K} = \\left(\\frac{\\frac{33}{100}}{\\frac{1}{50}+\\frac{1}{19}}\\right)^{\\frac{67}{100}} \\approx 9.57583\n",
    "$$\n",
    "\n",
    "Let’s verify this with Python and then use this steady state\n",
    "$ \\bar K $ as our initial capital stock $ K_0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcdc2a8f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ρ = 1 / pp.β - 1\n",
    "k_ss = pp.f_prime_inv(ρ+pp.δ)\n",
    "\n",
    "print(f'steady state for capital is: {k_ss}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "911c40d4",
   "metadata": {},
   "source": [
    "Now we plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bcb9424",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_paths(pp, 0.3, k_ss, [150], k_ss=k_ss);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05350cc2",
   "metadata": {},
   "source": [
    "Evidently,  with a large value of\n",
    "$ T $, $ K_t $ stays near $ K_0 $ until $ t $ approaches $ T $ closely.\n",
    "\n",
    "Let’s see what the planner does when we set\n",
    "$ K_0 $ below $ \\bar K $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcb85aa0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_paths(pp, 0.3, k_ss/3, [150], k_ss=k_ss);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2387e7f",
   "metadata": {},
   "source": [
    "Notice how the planner pushes capital toward the steady state, stays\n",
    "near there for a while, then pushes $ K_t $ toward the terminal\n",
    "value $ K_{T+1} =0 $ when $ t $ closely approaches $ T $.\n",
    "\n",
    "The following graphs compare optimal outcomes as we vary $ T $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cfb08f8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_paths(pp, 0.3, k_ss/3, [150, 75, 50, 25], k_ss=k_ss);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42601fac",
   "metadata": {},
   "source": [
    "## A Turnpike Property\n",
    "\n",
    "The following calculation indicates that when  $ T $ is very large,\n",
    "the optimal capital stock stays close to\n",
    "its steady state value most of the time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b485d594",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_paths(pp, 0.3, k_ss/3, [250, 150, 50, 25], k_ss=k_ss);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a1e00f",
   "metadata": {},
   "source": [
    "In the above graphs, different colors  are associated with\n",
    "different horizons $ T $.\n",
    "\n",
    "Notice that as the horizon increases, the planner keeps $ K_t $\n",
    "closer to the steady state value $ \\bar K $ for longer.\n",
    "\n",
    "This pattern reflects a **turnpike** property of the steady state.\n",
    "\n",
    "A rule of thumb for the planner is\n",
    "\n",
    "- from $ K_0 $, push $ K_t $ toward\n",
    "  the steady state and stay close to the steady state until time approaches $ T $.  \n",
    "\n",
    "\n",
    "The planner accomplishes this by adjusting the saving rate $ \\frac{f(K_t) - C_t}{f(K_t)} $\n",
    "over time.\n",
    "\n",
    "Let’s calculate and  plot the saving rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a69664",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def saving_rate(pp, c_path, k_path):\n",
    "    'Given paths of c and k, computes the path of saving rate.'\n",
    "    production = pp.f(k_path[:-1])\n",
    "\n",
    "    return (production - c_path) / production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2066818",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_saving_rate(pp, c0, k0, T_arr, k_ter=0, k_ss=None, s_ss=None):\n",
    "\n",
    "    fix, axs = plt.subplots(2, 2, figsize=(12, 9))\n",
    "\n",
    "    c_paths, k_paths = plot_paths(pp, c0, k0, T_arr, k_ter=k_ter, k_ss=k_ss, axs=axs.flatten())\n",
    "\n",
    "    for i, T in enumerate(T_arr):\n",
    "        s_path = saving_rate(pp, c_paths[i], k_paths[i])\n",
    "        axs[1, 1].plot(s_path)\n",
    "\n",
    "    axs[1, 1].set(xlabel='t', ylabel='$s_t$', title='Saving rate')\n",
    "\n",
    "    if s_ss is not None:\n",
    "        axs[1, 1].hlines(s_ss, 0, np.max(T_arr), linestyle='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdf73caa",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_saving_rate(pp, 0.3, k_ss/3, [250, 150, 75, 50], k_ss=k_ss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd4279ce",
   "metadata": {},
   "source": [
    "## A Limiting Infinite Horizon Economy\n",
    "\n",
    "We want to set $ T = +\\infty $.\n",
    "\n",
    "The appropriate thing to do is to replace terminal condition\n",
    "[(34.12)](#equation-constraint4) with\n",
    "\n",
    "$$\n",
    "\\lim_{T \\rightarrow +\\infty} \\beta^T u'(C_T) K_{T+1} = 0 ,\n",
    "$$\n",
    "\n",
    "a condition that will be satisfied by a path that converges to an\n",
    "optimal steady state.\n",
    "\n",
    "We can approximate the optimal path by starting from an arbitrary initial\n",
    "$ K_0 $ and shooting towards the optimal steady state\n",
    "$ K $ at a large but finite $ T+1 $.\n",
    "\n",
    "In the following code, we do this for a large $ T $ and plot consumption, capital, and the\n",
    "saving rate.\n",
    "\n",
    "We know that in the steady state that the saving rate is constant\n",
    "and that $ \\bar s= \\frac{f(\\bar K)-\\bar C}{f(\\bar K)} $.\n",
    "\n",
    "From [(34.17)](#equation-feasibility-constraint) the steady state saving rate equals\n",
    "\n",
    "$$\n",
    "\\bar s =\\frac{ \\delta \\bar{K}}{f(\\bar K)}\n",
    "$$\n",
    "\n",
    "The steady state saving rate $ \\bar S = \\bar s f(\\bar K) $ is\n",
    "the amount required to offset capital depreciation each period.\n",
    "\n",
    "We first study optimal capital paths that start below the steady\n",
    "state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4262a9f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# steady state of saving rate\n",
    "s_ss = pp.δ * k_ss / pp.f(k_ss)\n",
    "\n",
    "plot_saving_rate(pp, 0.3, k_ss/3, [130], k_ter=k_ss, k_ss=k_ss, s_ss=s_ss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c27eda89",
   "metadata": {},
   "source": [
    "Since $ K_0<\\bar K $, $ f'(K_0)>\\rho +\\delta $.\n",
    "\n",
    "The planner chooses a positive saving rate that is higher than  the steady state\n",
    "saving rate.\n",
    "\n",
    "Note that $ f''(K)<0 $, so as $ K $ rises, $ f'(K) $ declines.\n",
    "\n",
    "The planner slowly lowers the saving rate until reaching a steady\n",
    "state in which $ f'(K)=\\rho +\\delta $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2213354c",
   "metadata": {},
   "source": [
    "## Stable Manifold and Phase Diagram\n",
    "\n",
    "We now describe a classic diagram that describes an optimal $ (K_{t+1}, C_t) $ path.\n",
    "\n",
    "The diagram has $ K $ on the ordinate axis and $ C $ on the coordinate axis.\n",
    "\n",
    "Given an arbitrary and fixed  $ K $, a fixed point $ C $ of the consumption Euler equation  [(34.15)](#equation-eq-consn-euler)\n",
    "satisfies\n",
    "\n",
    "$$\n",
    "C=C\\left(\\beta\\left[f^{\\prime}\\left(f\\left(K\\right)+\\left(1-\\delta\\right)K-C\\right)+\\left(1-\\delta\\right)\\right]\\right)^{1/\\gamma}\n",
    "$$\n",
    "\n",
    "which implies\n",
    "\n",
    "\n",
    "<a id='equation-eq-tildec'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "C &=f\\left(K\\right)+\\left(1-\\delta\\right)K-f^{\\prime-1}\\left(\\frac{1}{\\beta}-\\left(1-\\delta\\right)\\right)  \\\\\n",
    " &\\equiv \\tilde{C} \\left(K\\right)\n",
    "\\end{aligned} \\tag{34.18}\n",
    "$$\n",
    "\n",
    "A positive fixed point  $ C = \\tilde C(K) $ exists only if $ f\\left(K\\right)+\\left(1-\\delta\\right)K-f^{\\prime-1}\\left(\\frac{1}{\\beta}-\\left(1-\\delta\\right)\\right)>0 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "994e317f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def C_tilde(K, pp):\n",
    "\n",
    "    return pp.f(K) + (1 - pp.δ) * K - pp.f_prime_inv(1 / pp.β - 1 + pp.δ)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f80e6c88",
   "metadata": {},
   "source": [
    "Next note that given a time-invariant  arbitrary $ C $,  a fixed point $ K $ of the feasibility condition  [(34.5)](#equation-allocation) solves the following equation\n",
    "\n",
    "$$\n",
    "K = f(K) + (1 - \\delta K) - C .\n",
    "$$\n",
    "\n",
    "A fixed point of the above equation is described by  a function\n",
    "\n",
    "\n",
    "<a id='equation-eq-tildek'></a>\n",
    "$$\n",
    "K = \\tilde K(C) \\tag{34.19}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37fb7584",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def K_diff(K, C, pp):\n",
    "    return pp.f(K) - pp.δ * K - C\n",
    "\n",
    "@jit\n",
    "def K_tilde(C, pp):\n",
    "\n",
    "    res = brentq(K_diff, 1e-6, 100, args=(C, pp))\n",
    "\n",
    "    return res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c714daf7",
   "metadata": {},
   "source": [
    "A  steady state $ \\left(K_s, C_s\\right) $ is a pair $ (K,C) $ that  satisfies both equations [(34.18)](#equation-eq-tildec) and [(34.19)](#equation-eq-tildek).\n",
    "\n",
    "It is thus the intersection of the two curves $ \\tilde{C} $ and $ \\tilde{K} $ that we’ll plot in Figure Fig. 34.1 below.\n",
    "\n",
    "We can compute $ K_s $ by solving the equation $ K_s = \\tilde{K}\\left(\\tilde{C}\\left(K_s\\right)\\right) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3367f4a8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def K_tilde_diff(K, pp):\n",
    "\n",
    "    K_out = K_tilde(C_tilde(K, pp), pp)\n",
    "\n",
    "    return K - K_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d543b4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "res = brentq(K_tilde_diff, 8, 10, args=(pp,))\n",
    "\n",
    "Ks = res.root\n",
    "Cs = C_tilde(Ks, pp)\n",
    "\n",
    "Ks, Cs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2e2a3b7",
   "metadata": {},
   "source": [
    "We can use the shooting algorithm to  compute  trajectories that approach $ \\left(K_s, C_s\\right) $.\n",
    "\n",
    "For a given $ K $, let’s compute $ \\vec{C} $ and $ \\vec{K} $ for a large $ T $ , e.g., $ =200 $.\n",
    "\n",
    "We compute  $ C_0 $ by the bisection algorithm that assures that  $ K_T=K_s $.\n",
    "\n",
    "Let’s compute  two trajectories towards $ \\left(K_s, C_s\\right) $ that  start from different sides of $ K_s $: $ \\bar{K}_0=1e-3<K_s<\\bar{K}_1=15 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eed258a2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "c_vec1, k_vec1 = bisection(pp, 5, 15, T=200, k_ter=Ks)\n",
    "c_vec2, k_vec2 = bisection(pp, 1e-3, 1e-3, T=200, k_ter=Ks)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a3f899c",
   "metadata": {},
   "source": [
    "The following code generates Figure Fig. 34.1, which is patterned on a graph that appears  on  page 411 of [[Intriligator, 2002](https://python.quantecon.org/zreferences.html#id3)].\n",
    "\n",
    "Figure Fig. 34.1 is a classic “phase plane” with  “state” variable $ K $ on the ordinate axis and “co-state” variable $ C $ on the coordinate axis.\n",
    "\n",
    "Figure Fig. 34.1 plots   three curves:\n",
    "\n",
    "- the blue line  graphs $ C = \\tilde C (K) $ of fixed points described by equation [(34.18)](#equation-eq-tildec).  \n",
    "- the red line graphs $ K = \\tilde K(C) $ of fixed points described by equation [(34.19)](#equation-eq-tildek)  \n",
    "- the green line graphs the stable traced out by paths that converge to the steady state starting from an arbitrary $ K_0 $ at time $ 0 $.  \n",
    "  - for a given $ K_0 $, the shooting algorithm sets $ C_0 $ to the coordinate on the green line in order to initiate a path that converges to the optimal steady state  \n",
    "  - the arrows on the green line show the direction in which  dynamics [(34.16)](#equation-eq-systemdynamics) push successive $ (K_{t+1}, C_t) $ pairs.  \n",
    "\n",
    "\n",
    "In addition to the three curves, Figure Fig. 34.1 plots  arrows that point where the  dynamics [(34.16)](#equation-eq-systemdynamics) drive the system  when, for a given $ K_0 $, $ C_0 $ is  not on the stable manifold depicted in the green line.\n",
    "\n",
    "- If $ C_0 $ is set below the green line for a given $ K_0 $, too much capital is accumulated  \n",
    "- If $ C_0 $ is set above the green line for a given $ K_0 $, too little capital is accumulated  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57a2450d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(7, 5))\n",
    "\n",
    "K_range = np.arange(1e-1, 15, 0.1)\n",
    "C_range = np.arange(1e-1, 2.3, 0.1)\n",
    "\n",
    "# C tilde\n",
    "ax.plot(K_range, [C_tilde(Ks, pp) for Ks in K_range], color='b')\n",
    "ax.text(11.8, 4, r'$C=\\tilde{C}(K)$', color='b')\n",
    "\n",
    "# K tilde\n",
    "ax.plot([K_tilde(Cs, pp) for Cs in C_range], C_range, color='r')\n",
    "ax.text(2, 1.5, r'$K=\\tilde{K}(C)$', color='r')\n",
    "\n",
    "# stable branch\n",
    "ax.plot(k_vec1[:-1], c_vec1, color='g')\n",
    "ax.plot(k_vec2[:-1], c_vec2, color='g')\n",
    "ax.quiver(k_vec1[5], c_vec1[5],\n",
    "          k_vec1[6]-k_vec1[5], c_vec1[6]-c_vec1[5],\n",
    "          color='g')\n",
    "ax.quiver(k_vec2[5], c_vec2[5],\n",
    "          k_vec2[6]-k_vec2[5], c_vec2[6]-c_vec2[5],\n",
    "          color='g')\n",
    "ax.text(12, 2.5, r'stable branch', color='g')\n",
    "\n",
    "# (Ks, Cs)\n",
    "ax.scatter(Ks, Cs)\n",
    "ax.text(Ks-1.2, Cs+0.2, '$(K_s, C_s)$')\n",
    "\n",
    "# arrows\n",
    "K_range = np.linspace(1e-3, 15, 20)\n",
    "C_range = np.linspace(1e-3, 7.5, 20)\n",
    "K_mesh, C_mesh = np.meshgrid(K_range, C_range)\n",
    "\n",
    "next_K, next_C = pp.next_k_c(K_mesh, C_mesh)\n",
    "ax.quiver(K_range, C_range, next_K-K_mesh, next_C-C_mesh)\n",
    "\n",
    "# infeasible consumption area\n",
    "ax.text(0.5, 5, \"infeasible\\n consumption\")\n",
    "\n",
    "ax.set_ylim([0, 7.5])\n",
    "ax.set_xlim([0, 15])\n",
    "\n",
    "ax.set_xlabel('$K$')\n",
    "ax.set_ylabel('$C$')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42c3b221",
   "metadata": {},
   "source": [
    "## Concluding Remarks\n",
    "\n",
    "In [Cass-Koopmans Competitive Equilibrium](https://python.quantecon.org/cass_koopmans_2.html),  we study a decentralized version of an economy with exactly the same\n",
    "technology and preference structure as deployed here.\n",
    "\n",
    "In that lecture, we replace the  planner of this lecture with Adam Smith’s **invisible hand**.\n",
    "\n",
    "In place of quantity choices made by the planner, there are market prices that are set by a *deus ex machina* from outside the model, a so-called invisible hand.\n",
    "\n",
    "Equilibrium market prices must reconcile distinct decisions that are made independently\n",
    "by a representative household and a representative firm.\n",
    "\n",
    "The relationship between a command economy like the one studied in this lecture and a market economy like that\n",
    "studied in [Cass-Koopmans Competitive Equilibrium](https://python.quantecon.org/cass_koopmans_2.html) is a foundational topic in general equilibrium theory and welfare economics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cec5300d",
   "metadata": {},
   "source": [
    "### Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa685fd",
   "metadata": {},
   "source": [
    "### Exercise 34.1\n",
    "\n",
    "- Plot the optimal consumption, capital, and saving paths when the\n",
    "  initial capital level begins at 1.5 times the steady state level\n",
    "  as we shoot towards the steady state at $ T=130 $.  \n",
    "- Why does the saving rate respond as it does?  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30188f2a",
   "metadata": {},
   "source": [
    "### Solution to[ Exercise 34.1](https://python.quantecon.org/#ck1_ex1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6d9a273",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_saving_rate(pp, 0.3, k_ss*1.5, [130], k_ter=k_ss, k_ss=k_ss, s_ss=s_ss)"
   ]
  }
 ],
 "metadata": {
  "date": 1736552526.713921,
  "filename": "cass_koopmans_1.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Cass-Koopmans Model"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}