Class_08_RCK.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ramsey-Cass-Koopmans (RCK) model\n",
    "\n",
    "The RCK model is an extension of the Solow-Swan model, pioneered by [Ramsey (1928)](http://www.jstor.org/stable/2224098) and extended by [Cass (1965)](http://www.jstor.org/stable/2295827) and [Koopmans (1965)](http://cowles.yale.edu/sites/default/files/files/pub/d01/d0163.pdf). Known also as the Neoclassical Growth Model, it is a \"workhorse\" of modern macroeconomics and forms the basis of the majority of modern works in economic growth and business cycles literature.\n",
    "\n",
    "While in the OLG model there were an infinite number of agents (since in every new period a new generation entered the economy) and they had a finite time horizon, in the Ramsey model there is a finite number of agents with an infinite time horizon.\n",
    "\n",
    "This can be justified by the following argument. In the (basic) OLG model parents are \"selfish\" and leave no assets to their children. However in the real world we see that at least some individuals leave sizeable bequests and support their children financially throughout their lives. This implies that parents care for their children's \"utility\". \n",
    "\n",
    "In the RCK model it is assumed that parents care about their children as much as they care about themselves. Therefore, each household is in fact a \"dynasty\" of individuals with an infinite planning horizon, and the welfare function of a household is constructed as follows:\n",
    "\n",
    "\\begin{align}\n",
    "U = \\sum_{t=0}^{\\infty}\\beta^{t}u\\left(c_{t}\\right)\n",
    "\\end{align}\n",
    "\n",
    "where $\\beta$ is a discount factor and $u\\left(\\cdot\\right)$ is an instantaneous utility function (felicity function).\n",
    "\n",
    "Usually the instantaneous utility function is represented by a Constant Relative Risk Aversion (CRRA) function:\n",
    "\n",
    "\\begin{align}\n",
    "u\\left(c_{t}\\right) = \\frac{c_{t}^{1-\\sigma}-1}{1-\\sigma}\n",
    "\\end{align}\n",
    "\n",
    "where $\\sigma>0$ regulates the degree of risk aversion.\n",
    "\n",
    "We will first consider a RCK model with simplifying assumptions that yield an analytical solution, and then consider a general case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Firms' problem\n",
    "\n",
    "The problem of the firms will be the same for both simplified and general cases. The firms want to hire optimal quantities of labor and capital to maximize their profits:\n",
    "\n",
    "\\begin{align}\n",
    "\\max\\quad\\Pi_{t}=K_{t}^{\\alpha}\\left(A_{t}L_{t}\\right)^{1-\\alpha}-w_{t}L_{t}-r_{t}^{k}K_{t}\n",
    "\\end{align}\n",
    "\n",
    "First order conditions (FOCs):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\Pi_{t}}{\\partial L_{t}} & =\\left(1-\\alpha\\right)K_{t}^{\\alpha}A_{t}^{1-\\alpha}L_{t}^{-\\alpha}-w_{t}=0\\\\\n",
    "\\frac{\\partial\\Pi_{t}}{\\partial K_{t}} & =\\alpha K_{t}^{\\alpha-1}\\left(A_{t}L_{t}\\right)^{1-\\alpha}-r_{t}^{k}=0\n",
    "\\end{align}\n",
    "\n",
    "Simplify the above expressions:\n",
    "\n",
    "\\begin{align}\n",
    "w_{t} & =\\left(1-\\alpha\\right)A_{t}\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha}\\\\\n",
    "r_{t}^{k} & =\\alpha\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha-1}\n",
    "\\end{align}\n",
    "\n",
    "The above equations determine both the wage $w$ and the capital rental rate $r^{k}$ as functions of capital (per effective labor). The real interest rate is equal to the capital rental rate, less the capital depreciation rate:\n",
    "\n",
    "\\begin{align}\n",
    "r_{t} & =r_{t}^{k}-\\delta\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make graphs appear within the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Import numerical computations library\n",
    "import numpy as np\n",
    "\n",
    "# Import plotting library\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Import minimization function\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "# Import functions from AppliedMacroPlots.py\n",
    "from AppliedMacroPlots import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simplified model\n",
    "\n",
    "### Households' problem\n",
    "\n",
    "Let us now make the simplifying assumptions. There is no population and technology growth, and without loss of generality we can normalize $A=1$, so that capital per effective labor and capital per worker are equivalent. We will also assume total depreciation of capital ($\\delta=1$) and that the instantaneous utility function is logarithmic.\n",
    "\n",
    "The representative \"dynasty\" solves the following problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max\\quad & U=\\sum_{t=0}^{\\infty}\\beta^{t}\\log c_{t}\\\\\n",
    "\\text{subject to}\\quad & c_{t}+a_{t+1}=w_{t}+\\left(1+r_{t}\\right)a_{t}\n",
    "\\end{align}\n",
    "\n",
    "Write down the Lagrangian:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L}=\\sum_{t=0}^{\\infty}\\beta^{t}\\log c_{t}+\\sum_{t=0}^{\\infty}\\lambda_{t}\\left[w_{t}+\\left(1+r_{t}\\right)a_{t}-c_{t}-a_{t+1}\\right]\n",
    "\\end{align}\n",
    "\n",
    "Note that we have an infinity of time periods and Lagrange multipliers. Rewrite the Lagrangian so that it is easier to take derivatives:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L} & =\\ldots+\\beta^{t}\\log c_{t}+\\ldots\\\\\n",
    " & \\quad+\\ldots+\\lambda_{t}\\left[w_{t}+\\left(1+r_{t}\\right)a_{t}-c_{t}-a_{t+1}\\right]\\\\\n",
    " & \\quad+\\lambda_{t+1}\\left[w_{t+1}+\\left(1+r_{t+1}\\right)a_{t+1}-c_{t+1}-a_{t+2}\\right]+\\ldots\n",
    "\\end{align}\n",
    "\n",
    "In each period $t$ we can choose period $t$ consumption $c_{t}$ and end-of-period $t$ assets $a_{t+1}$.\n",
    "\n",
    "First order conditions (FOCs):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial c_{t}} & =\\beta^{t}\\frac{1}{c_{t}}+\\lambda_{t}\\left[-1\\right]=0\\\\\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial a_{t+1}} & =\\lambda_{t}\\left[-1\\right]+\\lambda_{t+1}\\left[\\left(1+r_{t+1}\\right)\\right]=0\n",
    "\\end{align}\n",
    "\n",
    "Simplify and rewrite:\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda_{t} & =\\beta^{t}\\frac{1}{c_{t}}\\\\\n",
    "\\lambda_{t} & =\\left(1+r_{t+1}\\right)\\lambda_{t+1}\n",
    "\\end{align}\n",
    "\n",
    "Join conditions:\n",
    "\n",
    "\\begin{align}\n",
    "\\beta^{t}\\frac{1}{c_{t}} & =\\left(1+r_{t+1}\\right)\\beta^{t+1}\\frac{1}{c_{t+1}}\\\\\n",
    "\\frac{1}{c_{t}} & =\\left(1+r_{t+1}\\right)\\beta\\frac{1}{c_{t+1}}\\\\\n",
    "c_{t+1} & =\\beta\\left(1+r_{t+1}\\right)c_{t}\n",
    "\\end{align}\n",
    "\n",
    "We have obtained the standard Euler equation.\n",
    "\n",
    "The optimal behavior of the households can then be restated as:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t+1} & =\\beta\\left(1+r_{t+1}\\right)c_{t}\\\\\n",
    "c_{t}+a_{t+1} & =w_{t}+\\left(1+r_{t}\\right)a_{t}\n",
    "\\end{align}\n",
    "\n",
    "And after some algebra the optimal consumpion in period $t$ can be expressed as:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}=\\frac{1}{1-\\beta}\\left[\\left(1+r_{t}\\right)a_{t}+\\sum_{i=0}^{\\infty}\\frac{w_{t+i}}{1+\\bar{r}_{t+i}}\\right]\n",
    "\\end{align}\n",
    "\n",
    "where:\n",
    "\n",
    "\\begin{align}\n",
    "1+\\bar{r}_{t+i}=\\begin{cases}\n",
    "1 & \\text{for }i=0\\\\\n",
    "\\left(1+r_{t+1}\\right)\\cdot\\ldots\\cdot\\left(1+r_{t+i}\\right) & \\text{for }i=1,2,\\ldots\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "As you can see, the optimal behavior of households requires them to formulate accurate forecasts of prices $w$ and $r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Equilibrium\n",
    "\n",
    "The above solution of the households' problem is still relatively general, as we haven't made full use of our simplifying assumptions yet. Consider now the general equilibrium, where all markets clear. \n",
    "\n",
    "Since there is no government and no international trade, domestic capital is the only asset which can be in positive net supply:\n",
    "\n",
    "\\begin{align}\n",
    "a_{t} = k_{t} \\quad \\text{for all } t\n",
    "\\end{align}\n",
    "\n",
    "And the budget constraint of the households can be rewritten as:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}+a_{t+1} & =w_{t}+\\left(1+r_{t}\\right)a_{t} \\\\\n",
    "c_{t}+k_{t+1}&=w_{t}+\\left(1+r_{t}\\right)k_{t}\n",
    "\\end{align}\n",
    "\n",
    "We can also use the expressions for prices (recall that $A=1$ and $\\delta=1$):\n",
    "\n",
    "\\begin{align}\n",
    "w_{t} & =\\left(1-\\alpha\\right)A_{t}\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha}=\\left(1-\\alpha\\right)k_{t}^{\\alpha}\\\\\n",
    "r_{t}^{k} & =\\alpha\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha-1}=\\alpha k_{t}^{\\alpha-1}\\\\\n",
    "r_{t} & =r_{t}^{k}-\\delta=\\alpha k_{t}^{\\alpha-1}-1\n",
    "\\end{align}\n",
    "\n",
    "And plug them into the budget constraint:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}+k_{t+1} & =\\left(1-\\alpha\\right)k_{t}^{\\alpha}+\\left(1+\\alpha k_{t}^{\\alpha-1}-1\\right)k_{t}\\\\\n",
    "c_{t}+k_{t+1} & =\\left(1-\\alpha\\right)k_{t}^{\\alpha}+\\alpha k_{t}^{\\alpha}\\\\\n",
    "c_{t}+k_{t+1} & =k_{t}^{\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "### Solution\n",
    "\n",
    "Finally, let us guess-and-verify that this economy behaves as a Solow-Swan economy, and the households save a constant fraction $s$ of their income:\n",
    "\n",
    "\\begin{align}\n",
    "k_{t+1} & =sy_{t}=s\\cdot k_{t}^{\\alpha}\\\\\n",
    "c_{t} & =\\left(1-s\\right)y_{t}=\\left(1-s\\right)\\cdot k_{t}^{\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "Use the Euler equation:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t+1} & =\\beta\\left(1+r_{t+1}\\right)c_{t}\\\\\n",
    "c_{t+1} & =\\beta\\alpha k_{t+1}^{\\alpha-1}c_{t}\\\\\n",
    "\\left(1-s\\right)\\cdot k_{t+1}^{\\alpha} & =\\beta\\alpha k_{t+1}^{\\alpha-1}\\left(1-s\\right)\\cdot k_{t}^{\\alpha}\\\\\n",
    "k_{t+1} & =\\beta\\alpha\\cdot k_{t}^{\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "Indeed, the saving rate is a constant:\n",
    "\n",
    "\\begin{align}\n",
    "s=\\alpha\\beta\n",
    "\\end{align}\n",
    "\n",
    "Note that $\\partial s / \\partial \\beta > 0$ and that $s \\leq \\alpha = s_{GR}$, with equality when $\\beta=1$ (households care just as much about the future and the present). This implies that the RCK economy is never dynamically inefficient.\n",
    "\n",
    "Under our assumptions, we can obtain an analytical, closed form solution (no need to forecast prices!) of the RCK model:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t} & =\\left(1-\\alpha\\beta\\right)\\cdot k_{t}^{\\alpha}\\\\\n",
    "k_{t+1} & =\\alpha\\beta\\cdot k_{t}^{\\alpha}\n",
    "\\end{align}\n",
    "\n",
    "Finally, let us find the steady state of the model:\n",
    "\n",
    "\\begin{align}\n",
    "k^{*} & =\\alpha\\beta\\left(k^{*}\\right)^{\\alpha}\\\\\n",
    "\\left(k^{*}\\right)^{1-\\alpha} & =\\alpha\\beta\\\\\n",
    "k^{*} & =\\left(\\alpha\\beta\\right)^{1/\\left(1-\\alpha\\right)}\\\\\n",
    "c^{*} & =\\left(1-\\alpha\\beta\\right)\\cdot \\left(\\alpha\\beta\\right)^{\\alpha/\\left(1-\\alpha\\right)}\n",
    "\\end{align}\n",
    "\n",
    "The first plot below illustrates the process of convergence to the steady state over time. Note that since $\\delta=1$, the convergence is very fast.\n",
    "\n",
    "The second plot is a new type of plot, called a phase diagram, and will illustrate the convergence to the steady state in the $(k, c)$ space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "α = 0.33\n",
    "β = 0.96\n",
    "δ = 1\n",
    "\n",
    "# Steady state\n",
    "k_star = (α*β)**(1/(1-α))\n",
    "c_star = (1-α*β) * k_star**α\n",
    "\n",
    "# Number of periods\n",
    "T = 10+1\n",
    "\n",
    "k_t = np.zeros(T)\n",
    "c_t = np.zeros(T)\n",
    "\n",
    "# Initial level of capital per worker\n",
    "k_0 = k_star/20\n",
    "\n",
    "# Solving the model forward\n",
    "k_t[0] = k_0\n",
    "c_t[0] = (1-α*β) * k_0**α\n",
    "\n",
    "for t in range(T-1):\n",
    "    k_t[t+1] = α*β * k_t[t]**α\n",
    "    c_t[t+1] = (1-α*β) * k_t[t+1]**α\n",
    "\n",
    "# Store results for future use\n",
    "k_simple, c_simple = k_t, c_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convergence plot\n",
    "plt.plot(k_t, lw=2, label='$k$')\n",
    "plt.plot(c_t, lw=2, label='$c$')\n",
    "\n",
    "plt.hlines(k_star, 0, T-1, color='C0', lw=1, linestyle='--')\n",
    "plt.hlines(c_star, 0, T-1, color='C1', lw=1, linestyle='--')\n",
    "\n",
    "plt.title('Convergence to steady state over time')\n",
    "plt.xlabel('Period')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.savefig('RCK_trans.pdf', transparent=True) #, bbox_inches='tight', pad_inches=0.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Phase diagram\n",
    "\n",
    "kk = np.linspace(0, 1.5*k_star, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star, 1000)\n",
    "\n",
    "plt.plot(kk, kk**α - δ*kk, lw=2, label='$k_{t+1}=k_{t}$')\n",
    "plt.plot(kk**0 * k_star, cc, lw=2, label='$c_{t+1}=c_{t}$')\n",
    "plt.plot(k_star, c_star, 'ro', label='Steady state')\n",
    "plt.plot(k_t, c_t, 'k-', label='Saddle path')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per worker $k$')\n",
    "plt.ylabel('Consumption per worker $c$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Phase diagram\n",
    "\n",
    "kk = np.linspace(0, 1.5*k_star, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star, 1000)\n",
    "\n",
    "plt.plot(kk, kk**α - δ*kk, lw=2, label='$k_{t+1}=k_{t}$')\n",
    "plt.plot(kk**0 * k_star, cc, lw=2, label='$c_{t+1}=c_{t}$')\n",
    "plt.plot(k_star, c_star, 'ro', label='Steady state')\n",
    "plt.plot(kk, (1-α*β)*kk**α, 'k-', label='Saddle path')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per worker $k$')\n",
    "plt.ylabel('Consumption per worker $c$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "\n",
    "plt.savefig('RCK_simple.pdf', transparent=True) #, bbox_inches='tight', pad_inches=0.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phase diagram shows a number of conditions. \n",
    "\n",
    "The first one, $k_{t+1}=k_{t}$, describes a set of combinations of $k$ and $c$ for which capital per worker does not change over time. We can find them via the budget/resource constraint:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}+k_{t+1} & =k_{t}^{\\alpha}\\\\\n",
    "c_{t} &= k_{t}^{\\alpha} - k_{t+1}\n",
    "\\end{align}\n",
    "\n",
    "If $k_{t+1}=k_{t}=k$:\n",
    "\n",
    "\\begin{align}\n",
    "c = k^{\\alpha} - k\n",
    "\\end{align}\n",
    "\n",
    "The second condition, $c_{t+1}=c_{t}$, describes a set of combinations of $k$ and $c$ for which consumption per worker does not change over time. We can find them via the Euler equation:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t+1} & =\\beta\\left(1+r_{t+1}\\right)c_{t}\\\\\n",
    "1 &= \\beta\\left(1+r_{t+1}\\right)\\\\\n",
    "r_{t+1} &= \\frac{1}{\\beta} - 1\\\\\n",
    "\\alpha k_{t+1}^{\\alpha-1}-1  &= \\frac{1}{\\beta} - 1 \\\\\n",
    "k_{t+1}^{\\alpha-1} &= \\frac{1}{\\alpha\\beta}\\\\\n",
    "k_{t+1} &= \\left(\\alpha\\beta\\right)^{1/\\left(1-\\alpha\\right)} = k^{*}\\\\\n",
    "\\end{align}\n",
    "\n",
    "That means that when the capital per worker is at its steady state level, $c_{t+1}=c_{t}$ for all $c$.\n",
    "\n",
    "Obviously, the intersection of $k_{t+1}=k_{t}$ and $c_{t+1}=c_{t}$ is by definition the steady state of the model. \n",
    "\n",
    "Finally, the \"saddle path\" is our model's solution. Importantly, there is only one path that leads from any possible level of capital per worker $k$ to the steady state. This guarantees that the model has a unique solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neoclassical Growth Model\n",
    "\n",
    "We now turn to analyzing the full RCK/NGM model. Population grows at rate $n$, while technology improves at rate $g$. The firms' problem is the same as before. What is a bit different is the households' problem.\n",
    "\n",
    "### Households' problem\n",
    "\n",
    "The representative \"dynasty\" solves the following problem:\n",
    "\n",
    "\\begin{align}\n",
    "\\max\\quad & U=\\sum_{t=0}^{\\infty}\\beta^{t}\\cdot\\frac{c_{t}^{1-\\sigma}-1}{1-\\sigma}\\\\\n",
    "\\text{subject to}\\quad & c_{t}+\\left(1+n\\right)a_{t+1}=w_{t}+\\left(1+r_{t}\\right)a_{t}\n",
    "\\end{align}\n",
    "\n",
    "Lagrangian:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L}=\\sum_{t=0}^{\\infty}\\beta^{t}\\cdot\\frac{c_{t}^{1-\\sigma}-1}{1-\\sigma}+\\sum_{t=0}^{\\infty}\\lambda_{t}\\left[w_{t}+\\left(1+r_{t}\\right)a_{t}-c_{t}-\\left(1+n\\right)a_{t+1}\\right]\n",
    "\\end{align}\n",
    "\n",
    "Expand the Lagrangian:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathcal{L} & =\\ldots+\\beta^{t}\\cdot\\frac{c_{t}^{1-\\sigma}-1}{1-\\sigma}+\\ldots\\\\\n",
    " & \\quad+\\ldots+\\lambda_{t}\\left[w_{t}+\\left(1+r_{t}\\right)a_{t}-c_{t}-\\left(1+n\\right)a_{t+1}\\right]\\\\\n",
    " & \\quad+\\lambda_{t+1}\\left[w_{t+1}+\\left(1+r_{t+1}\\right)a_{t+1}-c_{t+1}-\\left(1+n\\right)a_{t+2}\\right]+\\ldots\n",
    "\\end{align}\n",
    "\n",
    "First order conditions (FOCs):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial c_{t}} & =\\beta^{t}c_{t}^{-\\sigma}+\\lambda_{t}\\left[-1\\right]=0\\\\\n",
    "\\frac{\\partial\\mathcal{L}}{\\partial a_{t+1}} & =\\lambda_{t}\\left[-\\left(1+n\\right)\\right]+\\lambda_{t+1}\\left[\\left(1+r_{t+1}\\right)\\right]=0\n",
    "\\end{align}\n",
    "\n",
    "Simplify and rewrite:\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda_{t} & =\\beta^{t}c_{t}^{-\\sigma}\\\\\n",
    "\\left(1+n\\right)\\lambda_{t} & =\\left(1+r_{t+1}\\right)\\lambda_{t+1}\n",
    "\\end{align}\n",
    "\n",
    "Join conditions:\n",
    "\n",
    "\\begin{align}\n",
    "\\left(1+n\\right)\\beta^{t}c_{t}^{-\\sigma} & =\\left(1+r_{t+1}\\right)\\beta^{t+1}c_{t+1}^{-\\sigma}\\\\\n",
    "c_{t}^{-\\sigma} & =\\frac{1+r_{t+1}}{1+n}\\beta c_{t+1}^{-\\sigma}\\\\\n",
    "\\left(\\frac{c_{t+1}}{c_{t}}\\right)^{\\sigma} & =\\beta\\frac{1+r_{t+1}}{1+n}\\\\\n",
    "\\frac{c_{t+1}}{c_{t}} & =\\left(\\beta\\frac{1+r_{t+1}}{1+n}\\right)^{1/\\sigma}\\\\\n",
    "c_{t+1} & =\\left(\\beta\\frac{1+r_{t+1}}{1+n}\\right)^{1/\\sigma}c_{t}\n",
    "\\end{align}\n",
    "\n",
    "We have obtained the Euler equation. Note that if $n=0$ and $\\sigma=1$, we go back to our usual form of Euler equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Equilibrium\n",
    "\n",
    "Again, since there is no government and no international trade, domestic capital is the only asset which can be in positive net supply:\n",
    "\n",
    "\\begin{align}\n",
    "a_{t} = k_{t} \\quad \\text{for all } t\n",
    "\\end{align}\n",
    "\n",
    "And the budget constraint of the households can be rewritten as:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}+\\left(1+n\\right)a_{t+1} & =w_{t}+\\left(1+r_{t}\\right)a_{t} \\\\\n",
    "c_{t}+\\left(1+n\\right)k_{t+1}&=w_{t}+\\left(1+r_{t}\\right)k_{t}\n",
    "\\end{align}\n",
    "\n",
    "We can also use the expressions for prices (but this time technology improves over time, so we make use of per effective labor variables):\n",
    "\n",
    "\\begin{align}\n",
    "w_{t} & =\\left(1-\\alpha\\right)A_{t}\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha}=\\left(1-\\alpha\\right)A_{t}\\hat{k}_{t}^{\\alpha}\\\\\n",
    "r_{t}^{k} & =\\alpha\\left(\\frac{K_{t}}{A_{t}L_{t}}\\right)^{\\alpha-1}=\\alpha \\hat{k}_{t}^{\\alpha-1}\\\\\n",
    "r_{t} & =r_{t}^{k}-\\delta=\\alpha \\hat{k}_{t}^{\\alpha-1}-\\delta\n",
    "\\end{align}\n",
    "\n",
    "And plug them into the budget constraint:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t}+\\left(1+n\\right)k_{t+1} & =\\left(1-\\alpha\\right)A_{t}\\hat{k}_{t}^{\\alpha}+(1+\\alpha\\hat{k}_{t}^{\\alpha-1}-\\delta)k_{t}\\quad|\\quad:A_{t}\\\\\n",
    "\\hat{c}_{t}+\\left(1+n\\right)\\frac{k_{t+1}}{A_{t}} & =\\left(1-\\alpha\\right)\\hat{k}_{t}^{\\alpha}+(1+\\alpha\\hat{k}_{t}^{\\alpha-1}-\\delta)\\hat{k}_{t}\\\\\n",
    "\\hat{c}_{t}+\\left(1+n\\right)\\frac{A_{t+1}}{A_{t}}\\frac{k_{t+1}}{A_{t+1}} & =\\left(1-\\alpha\\right)\\hat{k}_{t}^{\\alpha}+\\alpha\\hat{k}_{t}^{\\alpha}+\\left(1-\\delta\\right)\\hat{k}_{t}\\\\\n",
    "\\left(1+n\\right)\\left(1+g\\right)\\hat{k}_{t+1} & =\\hat{k}_{t}^{\\alpha}-\\hat{c}_{t}+\\left(1-\\delta\\right)\\hat{k}_{t}\\\\\n",
    "\\hat{k}_{t+1} & =\\frac{\\hat{k}_{t}^{\\alpha}-\\hat{c}_{t}+\\left(1-\\delta\\right)\\hat{k}_{t}}{\\left(1+n\\right)\\left(1+g\\right)}\n",
    "\\end{align}\n",
    "\n",
    "The $\\hat{k}_{t}^{\\alpha}-\\hat{c}_{t}$ part represents gross investment, and is the analogue of $s\\cdot\\hat{k}_{t}^{\\alpha}$ in the Solow-Swan and OLG models.\n",
    "\n",
    "Plug in the interest rate into the Euler equation and rewrite it in per effective labor terms:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t+1} & =\\left(\\beta\\frac{1+r_{t+1}}{1+n}\\right)^{1/\\sigma}c_{t}\\quad|\\quad:A_{t}\\\\\n",
    "\\frac{c_{t+1}}{A_{t}} & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\frac{c_{t}}{A_{t}}\\\\\n",
    "\\frac{A_{t+1}}{A_{t}}\\frac{c_{t+1}}{A_{t+1}} & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\hat{c}_{t}\\\\\n",
    "\\left(1+g\\right)\\hat{c}_{t+1} & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\hat{c}_{t}\\\\\n",
    "\\hat{c}_{t+1} & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\frac{\\hat{c}_{t}}{\\left(1+g\\right)}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Steady state\n",
    "\n",
    "If $\\hat{c}_{t+1}=\\hat{c}_{t}$, then:\n",
    "\n",
    "\\begin{align}\n",
    "\\left(1+g\\right) & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\\\\n",
    "\\left(1+g\\right)^{\\sigma} & =\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\\\\n",
    "\\frac{\\left(1+g\\right)^{\\sigma}\\left(1+n\\right)}{\\beta} & =1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta\\\\\n",
    "\\alpha\\hat{k}_{t+1}^{\\alpha-1} & =\\frac{\\left(1+g\\right)^{\\sigma}\\left(1+n\\right)}{\\beta}-\\left(1-\\delta\\right)\\\\\n",
    "\\hat{k}_{t+1}^{\\alpha-1} & =\\frac{\\left(1+g\\right)^{\\sigma}\\left(1+n\\right)-\\beta\\left(1-\\delta\\right)}{\\alpha\\beta}\\\\\n",
    "\\hat{k}^{*} & =\\left[\\frac{\\alpha\\beta}{\\left(1+g\\right)^{\\sigma}\\left(1+n\\right)-\\beta\\left(1-\\delta\\right)}\\right]^{1/\\left(1-\\alpha\\right)}\n",
    "\\end{align}\n",
    "\n",
    "Note that if $\\beta=1$ and $\\sigma=1$, then we can express $\\hat{k}^{*}$ as:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{k}^{*}=\\left[\\frac{\\alpha}{\\left(1+g\\right)\\left(1+n\\right)-\\left(1-\\delta\\right)}\\right]^{1/\\left(1-\\alpha\\right)}=\\left[\\frac{\\alpha}{\\delta+n+g+ng}\\right]^{1/\\left(1-\\alpha\\right)}\n",
    "\\end{align}\n",
    "\n",
    "which is equivalent to the Solow-Swan model solution, given that the saving rate is at its golden rule level.\n",
    "\n",
    "If $\\hat{k}_{t+1}=\\hat{k}_{t}=\\hat{k}$, then:\n",
    "\n",
    "\\begin{align}\n",
    "\\left(1+n\\right)\\left(1+g\\right)\\hat{k}_{t+1} & =\\hat{k}_{t}^{\\alpha}-\\hat{c}_{t}+\\left(1-\\delta\\right)\\hat{k}_{t}\\\\\n",
    "\\left(1+n\\right)\\left(1+g\\right)\\hat{k} & =\\hat{k}^{\\alpha}-\\hat{c}+\\left(1-\\delta\\right)\\hat{k}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{c} & =\\hat{k}^{\\alpha}+\\left(1-\\delta\\right)\\hat{k}-\\left(1+n\\right)\\left(1+g\\right)\\hat{k}\\\\\n",
    "\\hat{c} & =\\hat{k}^{\\alpha}-\\left(\\delta+n+g+ng\\right)\\hat{k}\n",
    "\\end{align}\n",
    "\n",
    "And the steady state level of consumption per effective labor is given by:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{c}^{*} & =(\\hat{k}^{*})^{\\alpha}-\\left(\\delta+n+g+ng\\right)\\hat{k}^{*}\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transition (saddle) path\n",
    "\n",
    "Along the transition to the steady state, the model variables have to obey the following conditions:\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{k}_{t+1} & =\\frac{\\hat{k}_{t}^{\\alpha}-\\hat{c}_{t}+\\left(1-\\delta\\right)\\hat{k}_{t}}{\\left(1+n\\right)\\left(1+g\\right)}\\\\\n",
    "\\hat{c}_{t+1} & =\\left(\\beta\\frac{1+\\alpha\\hat{k}_{t+1}^{\\alpha-1}-\\delta}{1+n}\\right)^{1/\\sigma}\\frac{\\hat{c}_{t}}{\\left(1+g\\right)}\n",
    "\\end{align}\n",
    "\n",
    "Solving the RCK model is equivalent to finding, for some initial level of capital per effective labor $\\hat{k}_{0}$, such $\\hat{c}_{0}$ that the above forward equations bring $\\hat{k}$ and $\\hat{c}$ to the steady state.\n",
    "\n",
    "We will solve the model by the shooting algorithm. We will try different values of $\\hat{c}_{0}$ and see whether they are too high or too low."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define needed functions\n",
    "\n",
    "def ProductionFunction(k, params):\n",
    "    α = params['α']\n",
    "    return k**α\n",
    "\n",
    "def InterestRate(k, params):\n",
    "    α = params['α']\n",
    "    δ = params['δ']\n",
    "    return α*k**(α-1) - δ\n",
    "    \n",
    "def SteadyState(params):\n",
    "    α = params['α']\n",
    "    δ = params['δ']\n",
    "    n = params['n']\n",
    "    g = params['g']\n",
    "    β = params['β']\n",
    "    σ = params['σ']\n",
    "    \n",
    "    r_star = (1+n) * (1+g)**σ / β - 1\n",
    "    k_star = (α/(r_star+δ))**(1/(1-α))\n",
    "    c_star = ProductionFunction(k_star, params) - (δ+n+g+n*g)*k_star\n",
    "    \n",
    "    return k_star, c_star\n",
    "\n",
    "def ResourceConstraint(k, c, params):\n",
    "    α = params['α']\n",
    "    δ = params['δ']\n",
    "    n = params['n']\n",
    "    g = params['g']\n",
    "    \n",
    "    return (ProductionFunction(k, params) + (1-δ)*k - c)/((1+n)*(1+g))\n",
    "\n",
    "def EulerEquation(k, c, params):\n",
    "    n = params['n']\n",
    "    g = params['g']\n",
    "    β = params['β']\n",
    "    σ = params['σ']\n",
    "    \n",
    "    k_next = ResourceConstraint(k, c, params)\n",
    "    \n",
    "    if k_next > 0:\n",
    "        r = InterestRate(k_next, params)\n",
    "        c_next = (β*(1+r)/(1+n))**(1/σ) * c / (1+g)\n",
    "        return c_next\n",
    "    else:\n",
    "        return 0\n",
    "\n",
    "def Constant_k(k, params):\n",
    "    δ = params['δ']\n",
    "    n = params['n']\n",
    "    g = params['g']\n",
    "    \n",
    "    return ProductionFunction(k, params) - (δ+n+g+n*g)*k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "params = {'α':0.33, 'δ':0.1, 'n':0.01, 'g':0.02, 'σ':2, 'β':0.96}\n",
    "\n",
    "# Steady state\n",
    "print(SteadyState(params))\n",
    "k_star, c_star = SteadyState(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forward equations and convergence criterion\n",
    "\n",
    "def Path(c_0, k_0, params, T=100):\n",
    "    \n",
    "    T += 1\n",
    "    \n",
    "    k_t = np.zeros(T)\n",
    "    c_t = np.zeros(T)\n",
    "    \n",
    "    k_t[0] = k_0\n",
    "    c_t[0] = c_0\n",
    "    \n",
    "    for t in range(T-1):\n",
    "        k_t[t+1] = ResourceConstraint(k_t[t], c_t[t], params)\n",
    "        if k_t[t+1] > 0:\n",
    "            c_t[t+1] = EulerEquation(k_t[t], c_t[t], params)\n",
    "        else:\n",
    "            k_t[t+1] = 0\n",
    "            c_t[t+1] = 0\n",
    "            \n",
    "    return k_t, c_t\n",
    "\n",
    "def Path_crit(c_0, k_0, params, T=100):\n",
    "    \n",
    "    k_t, c_t = Path(c_0, k_0, params, T)\n",
    "    k_star, c_star = SteadyState(params)\n",
    "    \n",
    "    ss_diff = np.sqrt((k_t-k_star)**2 + (c_t-c_star)**2)\n",
    "    \n",
    "    return np.min(ss_diff) + ss_diff[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we will make attempts to target the steady state, but in the first case the initial guess will be too low and in the second too high."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_0 = k_star / 20\n",
    "c_0 = 0.34\n",
    "\n",
    "k_t, c_t = Path(c_0, k_0, params)\n",
    "\n",
    "plt.plot(k_t, c_t, 'ko-')\n",
    "plt.plot(k_star, c_star, 'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_0 = k_star / 20\n",
    "c_0 = 0.35\n",
    "\n",
    "k_t, c_t = Path(c_0, k_0, params)\n",
    "\n",
    "plt.plot(k_t, c_t, 'ko-')\n",
    "plt.plot(k_star, c_star, 'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Confirm that convergence criterion attains minimum in the chosen range\n",
    "\n",
    "cc_0 = np.linspace(0.34, 0.345)\n",
    "\n",
    "for c_0 in cc_0:\n",
    "    plt.plot(c_0, Path_crit(c_0, k_0, params), 'bo')\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the function minimum, starting from an initial guess\n",
    "result = minimize(Path_crit, 0.34, args=(k_0, params, 100), method='Nelder-Mead')\n",
    "print(result)\n",
    "\n",
    "c_0 = result.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Phase diagram\n",
    "k_t, c_t = Path(c_0, k_0, params, 50)\n",
    "\n",
    "kk = np.linspace(0, 1.5*k_star, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star, 1000)\n",
    "\n",
    "plt.plot(kk, Constant_k(kk, params), lw=2, label='$\\hat{k}_{t+1}=\\hat{k}_{t}$')\n",
    "plt.plot(kk**0 * k_star, cc, lw=2, label='$\\hat{c}_{t+1}=\\hat{c}_{t}$')\n",
    "plt.plot(k_star, c_star, 'ro', label='Steady state')\n",
    "plt.plot(k_t, c_t, 'k-', label='Saddle path')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per effective labor $\\hat{k}$')\n",
    "plt.ylabel('Consumption per effective labor $\\hat{c}$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now verify that our code runs correctly, by solving again the simplified model, for which we have an analytic solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "params_simple = {'α':0.33, 'δ':1, 'n':0, 'g':0, 'σ':1, 'β':0.96}\n",
    "\n",
    "# Steady state\n",
    "print(SteadyState(params_simple))\n",
    "k_star_simple, c_star_simple = SteadyState(params_simple)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the function minimum, starting from an initial guess\n",
    "result_simple = minimize(Path_crit, 0.1, \n",
    "                         args=(k_star_simple / 20, params_simple, 100), \n",
    "                         method='Nelder-Mead')\n",
    "print(result_simple)\n",
    "\n",
    "c_0_simple = result_simple.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate transition path\n",
    "k_check, c_check = Path(c_0_simple, k_star_simple / 20, params_simple, 10)\n",
    "\n",
    "kk = np.linspace(0, 1.5*k_star_simple, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star_simple, 1000)\n",
    "\n",
    "plt.plot(kk, Constant_k(kk, params_simple), lw=2, label='$k_{t+1}=k_{t}$')\n",
    "plt.plot(kk**0 * k_star_simple, cc, lw=2, label='$c_{t+1}=c_{t}$')\n",
    "plt.plot(k_star_simple, c_star_simple, 'ro', label='Steady state')\n",
    "plt.plot(k_check, c_check, 'k-', lw=4, label='Saddle path')\n",
    "plt.plot(k_simple, c_simple, 'y-', lw=2, label='Analytic solution')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per worker $k$')\n",
    "plt.ylabel('Consumption per worker $c$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.plot(k_simple, (c_simple-c_check)**2, lw=2)\n",
    "plt.title('Squared errors')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems that our numerical procedure can recover the analytical solution.\n",
    "\n",
    "Let us now take a look at the behavior of capital, consumption and saving rate over time in our full model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_t = 1 - c_t/ProductionFunction(k_t, params)\n",
    "s_star = 1-c_star/ProductionFunction(k_star, params)\n",
    "\n",
    "plt.plot(k_t, lw=2)\n",
    "plt.hlines(k_star, 0, 50, lw=1, linestyle='--')\n",
    "plt.title('Capital per effective labor $\\hat{k}$')\n",
    "plt.xlabel('Period')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(c_t, lw=2)\n",
    "plt.hlines(c_star, 0, 50, lw=1, linestyle='--')\n",
    "plt.title('Consumption per effective labor $\\hat{c}$')\n",
    "plt.xlabel('Period')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(100*s_t, lw=2)\n",
    "plt.hlines(100*s_star, 0, 50, lw=1, linestyle='--')\n",
    "plt.title('Saving rate $s$ (%)')\n",
    "plt.xlabel('Period')\n",
    "plt.ylim(20, 30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually the RCK model generates decreasing saving rate along the transition path. The condition for when the saving rate is decreasing is summarized by:\n",
    "\n",
    "\\begin{align}\n",
    "s^{*} < 1/\\sigma\n",
    "\\end{align}\n",
    "\n",
    "while the BGP saving rate is a function of model parameters. The saving rates of countries rarely exceed 40%, and empirical estimates for $\\sigma$ usually place it around 2 (although with high uncertainty).\n",
    "\n",
    "Below I present a set of parameters for which the saving rate is increasing over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "params_dec = {'α':0.33, 'δ':0.1, 'n':0.01, 'g':0.02, 'σ':5, 'β':0.96}\n",
    "\n",
    "# Steady state\n",
    "print(SteadyState(params_dec))\n",
    "k_star_dec, c_star_dec = SteadyState(params_dec)\n",
    "\n",
    "# Find the function minimum, starting from an initial guess\n",
    "result_dec = minimize(Path_crit, 0, \n",
    "                      args=(k_star_dec / 20, params_dec, 100), \n",
    "                      method='Nelder-Mead')\n",
    "print('\\n', result_dec)\n",
    "\n",
    "# Generate transition path\n",
    "c_0_dec = result_dec.x\n",
    "k_dec, c_dec = Path(c_0_dec, k_star_dec / 20, params_dec, 100)\n",
    "\n",
    "# Check if steady state was reached\n",
    "\n",
    "kk = np.linspace(0, 1.5*k_star_dec, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star_dec, 1000)\n",
    "\n",
    "plt.plot(kk, Constant_k(kk, params_dec), lw=2, label='$\\hat{k}_{t+1}=\\hat{k}_{t}$')\n",
    "plt.plot(kk**0 * k_star_dec, cc, lw=2, label='$\\hat{c}_{t+1}=\\hat{c}_{t}$')\n",
    "plt.plot(k_star_dec, c_star_dec, 'ro', label='Steady state')\n",
    "plt.plot(k_dec, c_dec, 'k-', label='Saddle path')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per effective labor $\\hat{k}$')\n",
    "plt.ylabel('Consumption per effective labor $\\hat{c}$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# Saving rate\n",
    "s_dec = 1 - c_dec/ProductionFunction(k_dec, params_dec)\n",
    "s_star_dec = 1-c_star_dec/ProductionFunction(k_star_dec, params_dec)\n",
    "\n",
    "plt.plot(100*s_dec[:51], lw=2)\n",
    "plt.hlines(100*s_star_dec, 0, 50, lw=1, linestyle='--')\n",
    "plt.title('Saving rate $s$ (%)')\n",
    "plt.xlabel('Period')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This last exercise also gives us insight on the shape of the transition path as a function of $\\sigma$. \n",
    "\n",
    "As can be seen below, transition path for higher $\\sigma$ consumers runs close to the $\\hat{k}_{t+1}=\\hat{k}_{t}$ condition.\n",
    "\n",
    "Also, higher $\\sigma$ consumers end up with lower consumption per effective worker in the steady state, provided that $g>0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "params_low = {'α':0.33, 'δ':0.1, 'n':0.01, 'g':0.02, 'σ':0.2, 'β':0.96}\n",
    "\n",
    "# Steady state\n",
    "print(SteadyState(params_low))\n",
    "k_star_low, c_star_low = SteadyState(params_low)\n",
    "\n",
    "# Find the function minimum, starting from an initial guess\n",
    "result_low = minimize(Path_crit, 0, \n",
    "                      args=(k_star_dec / 20, params_low, 50), \n",
    "                      method='Nelder-Mead')\n",
    "print('\\n', result_low)\n",
    "\n",
    "# Generate transition path\n",
    "c_0_low = result_low.x\n",
    "k_low, c_low = Path(c_0_low, k_star_dec / 20, params_low, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kk = np.linspace(0, 1.5*k_star_low, 1000)\n",
    "cc = np.linspace(0, 1.5*c_star_low, 1000)\n",
    "\n",
    "plt.plot(kk, Constant_k(kk, params_low), lw=2)\n",
    "\n",
    "plt.plot(kk**0 * k_star_low, cc, 'r--', lw=2)\n",
    "plt.plot(k_star_low, c_star_low, 'ro', label='Steady state ($\\sigma_{L}$)')\n",
    "plt.plot(k_low, c_low, 'r-', label='Saddle path ($\\sigma_{L}$)')\n",
    "\n",
    "plt.plot(kk**0 * k_star_dec, cc, 'b--', lw=2)\n",
    "plt.plot(k_star_dec, c_star_dec, 'bo', label='Steady state ($\\sigma_{L}$)')\n",
    "plt.plot(k_dec, c_dec, 'b-', label='Saddle path ($\\sigma_{L}$)')\n",
    "\n",
    "plt.title('Phase diagram of the RCK model')\n",
    "plt.xlabel('Capital per effective labor $\\hat{k}$')\n",
    "plt.ylabel('Consumption per effective labor $\\hat{c}$')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "# plt.xlim(0, 4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare consumption paths\n",
    "\n",
    "plt.plot(c_t, 'r')\n",
    "plt.plot(c_dec[:51], 'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this context, risk aversion can be related to the standard deviation of consumption paths: consumers with high $\\sigma$ start with higher and end up with lower consumption that consumers with low $\\sigma$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extension: Stone-Geary preferences\n",
    "\n",
    "The standard RCK model predicts that for the \"usual\" parameter values saving rates will decline with the level of GDP per worker. This also counterfactually implies that countries with lowest levels of GDP per worker should also have highest saving / investment rates.\n",
    "\n",
    "To correct for this issue, we may introduce a modification to the utility function that requires that an individual has to consume at least a certain minimum to survive, and only the excess consumption yields utility. This is captured by the Stone-Geary function and can be expressed as:\n",
    "\n",
    "\\begin{align}\n",
    "u\\left(c_{t}\\right)=\\frac{\\left(c_{t}-\\bar{c}\\right)^{1-\\sigma}-1}{1-\\sigma}\n",
    "\\end{align}\n",
    "\n",
    "and the corresponding Euler equation reads:\n",
    "\n",
    "\\begin{align}\n",
    "c_{t+1}= \\left(\\beta \\frac{1+r_{t+1}}{1+n}\\right)^{1/\\sigma} \\left(c_{t}-\\bar{c}\\right) + \\bar{c}\n",
    "\\end{align}\n",
    "\n",
    "In economies where $y\\approx\\bar{c}$ (output per person is close to the survival requirement) households will save very little. Only when output per person becomes sufficiently high, those households will start saving more (even if $\\sigma<1/s^{*}$). This generates an \"inverse-U\" convergence pattern: the least developed countries exhibit relatively low growth rates, while the middle income countries grow the fastest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EulerEquation_SG(k, c, params):\n",
    "    n = params['n']\n",
    "    g = params['g']\n",
    "    β = params['β']\n",
    "    σ = params['σ']\n",
    "    c_min = params['c_min']\n",
    "    \n",
    "    k_next = ResourceConstraint(k, c, params)\n",
    "    \n",
    "    if k_next > 0:\n",
    "        r = InterestRate(k_next, params)\n",
    "        c_next = (β*(1+r)/(1+n))**(1/σ) * (c-c_min) / (1+g) + c_min\n",
    "        return c_next\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forward equations and convergence criterion\n",
    "\n",
    "def Path_SG(c_0, k_0, params, T=100):\n",
    "    \n",
    "    T += 1\n",
    "    \n",
    "    k_t = np.zeros(T)\n",
    "    c_t = np.zeros(T)\n",
    "    \n",
    "    k_t[0] = k_0\n",
    "    c_t[0] = c_0\n",
    "    \n",
    "    for t in range(T-1):\n",
    "        k_t[t+1] = ResourceConstraint(k_t[t], c_t[t], params)\n",
    "        if k_t[t+1] > 0:\n",
    "            c_t[t+1] = EulerEquation_SG(k_t[t], c_t[t], params)\n",
    "        else:\n",
    "            k_t[t+1] = 0\n",
    "            c_t[t+1] = 0\n",
    "            \n",
    "    return k_t, c_t\n",
    "\n",
    "def Path_crit_SG(c_0, k_0, params, T=100):\n",
    "    \n",
    "    k_t, c_t = Path_SG(c_0, k_0, params, T)\n",
    "    k_star, c_star = SteadyState(params)\n",
    "    \n",
    "    ss_diff = np.sqrt((k_t-k_star)**2 + (c_t-c_star)**2)\n",
    "    \n",
    "    return np.min(ss_diff) + ss_diff[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for Stone-Geary\n",
    "params_SG = {'α':0.33, 'δ':0.1, 'n':0.01, 'g':0, 'σ':2, 'β':0.96, 'c_min':0.49}\n",
    "\n",
    "# Steady state\n",
    "print(SteadyState(params_SG))\n",
    "k_star_SG, c_star_SG = SteadyState(params_SG)\n",
    "\n",
    "# Initial level of capital\n",
    "k_0_SG = k_star_SG / 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the function minimum, starting from an initial guess\n",
    "result_SG = minimize(Path_crit_SG, 0.5,\n",
    "                     args=(k_0_SG, params_SG, 100),\n",
    "                     method='Nelder-Mead')\n",
    "print(result_SG)\n",
    "\n",
    "c_0_SG = result_SG.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_t_SG, c_t_SG = Path_SG(c_0_SG, k_0_SG, params_SG, 100)\n",
    "\n",
    "y_t_SG = ProductionFunction(k_t_SG, params_SG)\n",
    "\n",
    "s_t_SG = 1 - c_t_SG/y_t_SG\n",
    "\n",
    "g_y_t_SG = 100*(y_t_SG[1:]/y_t_SG[:-1]-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(k_t_SG, lw=2)\n",
    "plt.hlines(k_star_SG, 0, 100, lw=1, linestyle='--')\n",
    "plt.title('Capital per person $k$')\n",
    "plt.xlabel('Period')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(c_t_SG, lw=2)\n",
    "plt.hlines(c_star_SG, 0, 100, lw=1, linestyle='--')\n",
    "plt.title('Consumption per person $c$')\n",
    "plt.xlabel('Period')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(100*s_t_SG, lw=2)\n",
    "plt.hlines(100*s_t_SG[-1], 0, 100, lw=1, linestyle='--')\n",
    "plt.title('Saving rate $s$ (%)')\n",
    "plt.xlabel('Period')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(g_y_t_SG, lw=2)\n",
    "plt.title('Growth rate of GDP per person $g_{y}$ (%)')\n",
    "plt.xlabel('Period')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.log(1+y_t_SG), 100*s_t_SG, lw=3)\n",
    "plt.title('Saving rate vs income for Stone-Geary utility function')\n",
    "plt.xlabel('Logarithm of GDP per person')\n",
    "plt.ylabel('Saving rate (%)')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(np.log(1+y_t_SG[:-1]), g_y_t_SG, lw=3)\n",
    "plt.title('Saving rate vs income for Stone-Geary utility function')\n",
    "plt.xlabel('Logarithm of GDP per person')\n",
    "plt.ylabel('Annual growth rate (%)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}