{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#*Numerical Methods for Economists:* The Ramsey-Cass-Koopmans Model\n", "\n", "In the previous lab, we learned how to simulate a deterministic, continuous-time version of the Solow growth model where the savings rate of the economy was taken to be exogenous. In this lab, we will learn how to solve and simulate a deteministic, continuous time version of the neoclassical optimal growth model, often referred to as the Ramsey-Cass-Koopmans model, where the savings rate is endogenously chosen optimally by a representative household. \n", "\n", "A preview of what we will cover in this lab:\n", "\n", "* **Task 1:** Coding the optimal growth model in Python\n", "* **Task 2:** Calibrating the model\n", "* **Task 3:** Analyzing the phase diagram\n", "* **Task 4:** Solving and simulating the model\n", "* **Task 5:** Impulse response functions\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# These libaries will be doing the bulk of the heavy lifting!\n", "import numpy as np\n", "import pandas as pd\n", "import sympy as sp\n", "from scipy import linalg\n", "\n", "# these models contain all of the code\n", "import ramsey\n", "import pwt\n", "\n", "# matplotlib is for graphics\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# this function will allow for pretty LaTex printing within the IPython notebook\n", "from sympy.interactive import printing\n", "printing.init_printing(use_latex=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Task 1: Coding the optimal growth model in Python\n", "\n", "In this task you will learn how to program a continuous time version of the optimal growth model in Python. Before we starting writing code, we should stop and think about the \"primitives\" (i.e., basic bulding blocks) of an optimal growth model. \n", "\n", "## Coding $f(k)$ and $f'(k)$:\n", "In the previous lab we began coding up the Solow model by writing a Python function representing the intensive form of a [Cobb-Douglas](http://en.wikipedia.org/wiki/Cobb%E2%80%93Douglas_production_function) production function. In this lab we will assume that $f$ has the more general constant elasticity of substitution (CES) form\n", "\n", "$$f(k(t)) = \\bigg[\\alpha k(t)^\\gamma + (1 - \\alpha)\\bigg]^{\\frac{1}{\\gamma}} \\tag{1.1a}$$\n", "\n", "where $0 < \\alpha < 1$ and $-\\infty < \\gamma < 1$. The parameter $\\gamma = \\frac{\\sigma - 1}{\\sigma}$ where $\\sigma$ is the elasticity of substitution between capital and effective labor. Recall that the CES production technology nests the Cobb-Douglas function\n", "\n", "$$ f(k(t)) = k(t)^{\\alpha} \\tag{1.1b}$$\n", "\n", "as a special case when $\\sigma=1$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1:\n", "\n", "Use the information above to complete the [Python function](http://www.greenteapress.com/thinkpython/html/thinkpython004.html) `ces_output` defined below. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ces_output(t, k, params):\n", " \"\"\"\n", " Constant elasticity of substitution (CES) production function.\n", "\n", " Arguments:\n", "\n", " t: (array) Time.\n", " k: (array) Capital (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " y: (array-like) Output (per person/ effective person)\n", "\n", " \"\"\"\n", " # extract params\n", " alpha = params['alpha']\n", " sigma = params['sigma']\n", " gamma = (sigma - 1) / sigma\n", " \n", " # nest Cobb-Douglas as special case\n", " if gamma == 0:\n", " y = # INSERT CODE HERE!!\n", " else:\n", " y = # INSERT CODE HERE!!\n", " \n", " return y" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution is as follows\n", "def ces_output(t, k, params):\n", " \"\"\"\n", " Constant elasticity of substitution (CES) production function.\n", "\n", " Arguments:\n", "\n", " t: (array) Time.\n", " k: (array) Capital (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " y: (array-like) Output (per person/ effective person)\n", "\n", " \"\"\"\n", " # extract params\n", " alpha = params['alpha']\n", " sigma = params['sigma']\n", " gamma = (sigma - 1) / sigma\n", " \n", " # nest Cobb-Douglas as special case\n", " if gamma == 0:\n", " y = k**alpha\n", " else:\n", " y = (alpha * k**gamma + (1 - alpha))**(1 / gamma)\n", " \n", " return y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also need to code a function defining the marginal product of capital per effective person. From your homework assignment, you should recall that \n", "\n", "$$ f'(k) = \\alpha k(t)^{\\gamma-1}\\bigg[\\alpha k(t)^{\\gamma} + (1-\\alpha)\\bigg]^{\\frac{1}{\\gamma} - 1} \\tag{1.2a}$$\n", "\n", "for $\\sigma \\neq 1$ and \n", "\n", "$$f'(k) = \\alpha k(t)^{\\alpha-1} \\tag{1.2b}$$\n", "\n", "for the Cobb-Douglas special case when $\\sigma = 1$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2:\n", "\n", "Use the information above to complete the [Python function](http://www.greenteapress.com/thinkpython/html/thinkpython004.html) `marginal_product_capital` defined below. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def marginal_product_capital(t, k, params):\n", " \"\"\"\n", " Marginal product of capital for CES production function.\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " k: (array-like) Capital (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " mpk: (array-like) Derivative of output with respect to capital, k.\n", "\n", " \"\"\"\n", " # extract params\n", " alpha = params['alpha']\n", " sigma = params['sigma']\n", " gamma = (sigma - 1) / sigma\n", " \n", " # nest Cobb-Douglas as special case\n", " if gamma == 0:\n", " mpk = # INSERT CODE HERE!!\n", " else:\n", " mpk = # INSERT CODE HERE!!\n", " \n", " return mpk\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution!\n", "def marginal_product_capital(t, k, params):\n", " \"\"\"\n", " Marginal product of capital for CES production function.\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " k: (array-like) Capital (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " mpk: (array-like) Derivative of output with respect to capital, k.\n", "\n", " \"\"\"\n", " # extract params\n", " alpha = params['alpha']\n", " sigma = params['sigma']\n", " gamma = (sigma - 1) / sigma\n", " \n", " # nest Cobb-Douglas as special case\n", " if gamma == 0:\n", " mpk = alpha * k**(alpha - 1)\n", " else:\n", " mpk = alpha * k**(gamma - 1) * (alpha * k**gamma + (1 - alpha))**((1 / gamma) - 1)\n", " \n", " return mpk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coding the equilibrium system of ODEs:\n", "\n", "From the lecture you should recall that the equilibrium of the continuous-time, optimal growth model is described by the following two-dimensional system of ordinary differential equations (ODEs).\n", "\n", "\\begin{align}\n", " \t\\dot{k} =& f(k(t)) - (n+g+\\delta)k(t) - c(t),\\ k(0) = k_0 \\tag{1.3}\\\\\n", "\t\\dot{c} =& \\left[\\frac{f'(k(t)) - \\delta - \\rho - \\theta g}{\\theta}\\right]c(t),\\ 0 < \\lim_{t\\rightarrow\\infty} |c(t)| < \\infty \\tag{1.4}\n", "\\end{align}\n", "\n", "Equation 1.3 is the equation of motion for capital per effective person and equation 1.4 is the consumption Euler equation. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3:\n", "\n", "Use the information above to complete the [Python functions](http://www.greenteapress.com/thinkpython/html/thinkpython004.html) `equation_motion_capital` and `consumption_euler_equation` defined below. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def equation_motion_capital(t, vec, params):\n", " \"\"\"\n", " Equation of motion for capital (per person/ effective person).\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " vec: (array-like) Vector of endogenous variables [k, c] where k is \n", " capital (per person/effective person) and c is consumption \n", " (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " k_dot: (array-like) Rate of change in the stock of captial \n", " (per person/effective person).\n", "\n", " \"\"\"\n", " # extract params\n", " n = params['n']\n", " g = params['g']\n", " delta = params['delta']\n", " \n", " # unpack the vector of endogenous variables\n", " k, c = vec\n", " \n", " # formula for the equation of motion for capital \n", " k_dot = # INSERT CODE HERE!!\n", " \n", " return k_dot\n", "\n", "def consumption_euler_equation(t, vec, params):\n", " \"\"\"\n", " Equation of motion for consumption (per person/ effective person).\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " vec: (array-like) Vector of endogenous variables [k, c] where k is \n", " capital (per person/effective person) and c is consumption \n", " (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " c_dot: (array-like) Rate of change in consumption (per person/effective \n", " person).\n", "\n", " \"\"\"\n", " # extract params\n", " g = params['g']\n", " delta = params['delta']\n", " rho = params['rho']\n", " theta = params['theta']\n", " \n", " # unpack the vector of endogenous variables\n", " k, c = vec\n", " \n", " # marginal product of capital\n", " mpk = # INSERT CODE HERE!!\n", "\n", " # formula for the equation of motion for consumption\n", " c_dot = # INSERT CODE HERE!!\n", " \n", " return c_dot" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solutions!\n", "def equation_motion_capital(t, vec, params):\n", " \"\"\"\n", " Equation of motion for capital (per person/ effective person).\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " vec: (array-like) Vector of endogenous variables [k, c] where k is \n", " capital (per person/effective person) and c is consumption \n", " (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " k_dot: (array-like) Rate of change in the stock of captial \n", " (per person/effective person).\n", "\n", " \"\"\"\n", " # extract params\n", " n = params['n']\n", " g = params['g']\n", " delta = params['delta']\n", " \n", " # unpack the vector of endogenous variables\n", " k, c = vec\n", " \n", " # formula for the equation of motion for capital \n", " k_dot = ces_output(t, k, params) - (n + g + delta) * k - c\n", " \n", " return k_dot\n", "\n", "def consumption_euler_equation(t, vec, params):\n", " \"\"\"\n", " Equation of motion for consumption (per person/ effective person).\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " vec: (array-like) Vector of endogenous variables [k, c] where k is \n", " capital (per person/effective person) and c is consumption \n", " (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " c_dot: (array-like) Rate of change in consumption (per person/effective \n", " person).\n", "\n", " \"\"\"\n", " # extract params\n", " g = params['g']\n", " delta = params['delta']\n", " rho = params['rho']\n", " theta = params['theta']\n", " \n", " # unpack the vector of endogenous variables\n", " k, c = vec\n", " \n", " # marginal product of capital\n", " mpk = marginal_product_capital(t, k, params)\n", "\n", " # formula for the equation of motion for consumption\n", " c_dot = ((mpk - delta - rho - theta * g) / theta) * c\n", " \n", " return c_dot" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Coding the instantaneous utility function, $u(c)$:\n", "\n", "Recall that the consumption Euler equation as defined above assumes that the representative household hold has constant relative risk aversion (CRRA) utility.\n", "\n", "$$u(c) = \\frac{c^{1 - \\theta} - 1}{1 - \\theta} \\tag{1.5a}$$\n", "\n", "where utility is a function of consumption per effective person, $c$. The parameter, $0 \\lt \\theta$, plays a dual role:\n", "\n", "* $\\theta$ is the agent's coefficient of relative risk aversion\n", "* $\\frac{1}{\\theta}$ is the agent's elasticity of inter-temporal substitution\n", "\n", "Things to remember about this utility function:\n", "\n", "When $0 \\lt \\theta$ but \"close\" to zero, consumer has almost linear (i.e., risk neutral) utility. As $\\theta$ increases, households becomes more risk averse (i.e., the curvature of the utility function increases!). Note that this increase in curvature becomes most apparent in the behavior of the utility function for levels of consumption less than 1. As $\\theta$ increase the \"penalty\" for consuming less than 1 unit becomes increasingly severe!\n", "\n", "Why is there a '-1' in the numerator come from? It has the effect of making total utility positive or negative depending on whether $c$ is greater than or less than one. Also \n", "\n", "$$\\lim_{\\theta \\rightarrow 1}\\ \\frac{c^{1 - \\theta} - 1}{1 - \\theta}\\ =\\ \\ln\\ c. \\tag{1.5b}$$\n", "\n", "It is also **not** necessary! Leaving off the '-1' will still yield CRRA preferences." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def crra_utility(c, params):\n", " \"\"\"\n", " The optimal growth model assumes that the representative individual has CRRA utility. \n", " utility. Individual utility is a function of consumption per effective person, c. The \n", " parameter, 0 < theta, plays a dual role: \n", " \n", " 1) theta is the agent's coefficient of relative risk aversion\n", " 2) 1 / theta is the agent's elasticity of inter-temporal substitution\n", " \n", " N.B.: If theta = 1, then CRRA utility is equivalent to log utility! \n", " \n", " Arguments:\n", " \n", " c: (array) Consumption per effective person\n", " \n", " params: (dict) Dictionary of model parameters.\n", " \n", " Returns: \n", " \n", " utility: (array) An array of values respresenting utility from consumption per \n", " effective worker.\n", " \n", " \"\"\"\n", " # extract the params\n", " theta = params['theta']\n", " \n", " if theta == 1:\n", " utility = np.log(c)\n", " else:\n", " utility = (c**(1 - theta) - 1) / (1 - theta)\n", " \n", " return utility" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot this utility function for different values of the elasticity of intertemporal substitution..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create a grid of points for plotting\n", "gridmax, gridsize = 6, 1000\n", "grid = np.linspace(0, gridmax, gridsize)\n", "\n", "# Create a new figure\n", "plt.figure(figsize=(8,6))\n", "\n", "# create an array of some values for theta\n", "theta_values = np.logspace(-2, 2, 5)\n", "\n", "# Plot utility for varying values of theta\n", "for theta in theta_values:\n", " params = {'theta':theta}\n", " plt.plot(grid, crra_utility(grid, params), '-', label=r'$\\theta=%.2f$' %theta)\n", "\n", "# Don't forget to label your axes!\n", "plt.xlabel('$c$', fontsize=15, family='serif')\n", "plt.ylabel('$u(c)$', rotation='horizontal', fontsize=15, family='serif')\n", "plt.axis([0, 6, -10, 5])\n", "\n", "# Add a title to the plot\n", "plt.title(r'CRRA Utility for different values of $\\theta$', fontsize=20, family='serif')\n", "\n", "# Add a legend\n", "plt.legend(loc=2, frameon=False)\n", "\n", "# Turn on grid lines\n", "plt.grid()\n", "\n", "# save and display figure!\n", "plt.savefig('graphics/CRRA-utility.png')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Coding the Jacobian matrix:\n", "The next piece of the model that we need to construct is the [Jacobian](http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant) matrix of partial derivatives. \n", "\n", "$$ \n", "\\mathcal{J} =\n", "\\begin{bmatrix}\n", "\\frac{\\partial \\dot{k}}{\\partial k} & \\frac{\\partial \\dot{k}}{\\partial c}\\\\\n", "\\frac{\\partial \\dot{c}}{\\partial k} & \\frac{\\partial \\dot{c}}{\\partial c}\n", "\\end{bmatrix} =\n", "\\begin{bmatrix}\n", "f'(k) - (n + g + \\delta) & -1\\\\\n", "\\frac{1}{\\theta}f''(k) & \\frac{f'(k) - \\delta - \\rho - \\theta g}{\\theta}\n", "\\end{bmatrix} \n", "\\tag{1.5}\n", "$$\n", "\n", "The Jacobian matrix plays several important roles in the optimal growth model (and many other economic models!):\n", "\n", "1. We can use $\\mathcal{J}$ to construct a linear approximation of the full non-linear model via a first-order Taylor expansion around the model's long-run steady-state.\n", "2. We can use the [eigenvalues and eigenvectors](http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of $\\mathcal{J}$ to assess the model's stability properties.\n", "3. $\\mathcal{J}$ is also an important input into finite-difference methods for solving boundary value problems (BVPs). \n", "\n", "Although the above matrix of partial derivatives looks fairly straightforward to compute by hand, you may recall from your homework that deriving an analytic expression for $f''$ with CES production is tedious and error prone. In general, a much better strategy is to get the computer to take complicated and error-prone derivatives for you! There are three general approaches that one can use in order to get a computer to take derivatives: [numeric differentiation](http://en.wikipedia.org/wiki/Numerical_differentiation), [symbolic differentiation](http://en.wikipedia.org/wiki/Symbolic_differentiation), or [automatic differentiation](http://en.wikipedia.org/wiki/Automatic_differentiation). \n", "\n", "While numerical differentiation is often very fast, routines that are simple to implement may fail to achieve a high level of accuracy. Automatic differentiation is the gold-standard but open-source libraries implmenting the techniques are rare (although there are several such libraries currently under development in Python). We will take a middle of the road approach and use symbolic differentiation via [SymPy](http://sympy.org/en/index.html) to compute the model's Jacobian." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# declare symbolic variables\n", "c, k = sp.symbols('c,k')\n", "\n", "# declare symbolic parameters\n", "n, g, alpha, delta, rho, sigma, theta = sp.symbols('n,g,alpha,delta,rho,sigma,theta')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We always want to re-use code whenever possible. We have already coded functions for output, marginal product of capital, the equation of motion for capital, and the consumption Euler equation. By being slightly clever, we can use Sympy to differentiate our existing functions!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a dictionary of symbolic parameters (recall that params is an input to our functions)\n", "symbolic_params = {'n':n, 'g':g, 'alpha':alpha, 'delta':delta, 'rho':rho, \n", " 'sigma':sigma, 'theta':theta}\n", "\n", "# wrap the model equations so that they are functions of k and c only!\n", "k_dot = lambda k, c: equation_motion_capital(0, [k, c], symbolic_params)\n", "c_dot = lambda k, c: consumption_euler_equation(0, [k, c], symbolic_params)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a SymPy matrix containing the wrapped functions\n", "F = sp.Matrix([k_dot(k, c), c_dot(k, c)])\n", "\n", "# now compute the Jacobian is trivial!\n", "symbolic_jacobian = F.jacobian([k, c])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# try to simplify as much as possible\n", "symbolic_jacobian.simplify()\n", "\n", "# display the symbolic jacobian for the general CES case\n", "symbolic_jacobian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to differentiate the Cobb-Douglas special case, we need to re-wrap the model equations with a new dictionary of parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set sigma = 1.0 in order to specify the Cobb-Douglas special case\n", "symbolic_params = {'n':n, 'g':g, 'alpha':alpha, 'delta':delta, 'rho':rho, \n", " 'sigma':1.0, 'theta':theta}\n", "\n", "k_dot = lambda k, c: equation_motion_capital(0, [k, c], symbolic_params)\n", "c_dot = lambda k, c: consumption_euler_equation(0, [k, c], symbolic_params)\n", "F = sp.Matrix([k_dot(k, c), c_dot(k, c)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# now computing the Jacobian \n", "symbolic_jacobian = F.jacobian([k, c])\n", "\n", "# try to simplify as much as possible\n", "symbolic_jacobian.simplify()\n", "\n", "# display the symbolic jacobian for the Cobb-Douglas case\n", "symbolic_jacobian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can actually combine all of the above into a single Python function as follows..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def ramsey_jacobian(t, vec, params):\n", " \"\"\"\n", " The Jacobian for the optimal growth model is computed using symbolic \n", " differentiation and then evaluated numerically.\n", "\n", " Arguments:\n", "\n", " t: (array-like) Time.\n", " vec: (array-like) Vector of endogenous variables [k, c] where k is \n", " capital (per person/effective person) and c is consumption \n", " (per person/effective person).\n", " params: (dict) Dictionary of parameter values.\n", " \n", " Returns:\n", "\n", " jac: (array-like) Jacobian matrix of partial derivatives.\n", "\n", " \"\"\"\n", " # declare two symbolic variables\n", " k, c = sp.symbols('k,c')\n", " \n", " # represent the system of equations as a SymPy matrix\n", " k_dot = lambda k, c: equation_motion_capital(t, [k, c], params)\n", " c_dot = lambda k, c: consumption_euler_equation(t, [k, c], params)\n", " F = sp.Matrix([k_dot(k, c), c_dot(k, c)])\n", "\n", " # now computing the Jacobian is trivial!\n", " symbolic_jacobian = F.jacobian([k, c])\n", " \n", " # lambdify the function in order to evaluate it numerically\n", " numeric_jacobian = sp.lambdify(args=[k, c], expr=symbolic_jacobian, modules='numpy') \n", " \n", " # output num_jac evaluated at vec\n", " evaluated_jacobian = np.array(numeric_jacobian(vec[0], vec[1]))\n", " \n", " return evaluated_jacobian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Coding the steady state values of $k$ and $c$\n", "Because the Ramsey Model is still relatively simply, we can write down an analytic expressions for the steady state value of capital per effective person and consumption per effective person in terms of the structural parameters of the model. \n", "\n", "In general, the steady state value for capital per effective person is the value of $k$ that solves\n", "\n", "$$ f'(k) = \\delta + \\rho + \\theta g. \\tag{1.6}$$\n", "\n", "For the CES case you should be able to show that \n", "\n", "$$k^* = (1 - \\alpha)^{\\frac{1}{\\gamma}}\\Bigg[\\left(\\frac{\\alpha}{\\delta + \\rho + \\theta g}\\right)^{\\frac{\\gamma}{\\gamma - 1}} - \\left(\\frac{\\alpha}{1 - \\alpha}\\right)\\Bigg]^{-\\frac{1}{\\gamma}} \\tag{1.7a}$$\n", "\n", "when $\\sigma \\ne 1$, and\n", "\n", "$$ k^* = \\left(\\frac{\\alpha}{\\delta+\\rho+\\theta g}\\right)^{\\frac{1}{1-\\alpha}} \\tag{1.7b}$$\n", "\n", "when $\\sigma = 1$. Once the steady state value for capital per effective person is known the steady state value of consumption per effective person can be computed using the equation of motion for capital per effective person.\n", "\n", "$$ c^* = f(k^*) - (n + g + \\delta)k^* \\tag{1.8}$$\n", "\n", "Now we need to code these expressions as Python functions..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def analytic_k_star(params): \n", " \"\"\"Steady-state level of capital (per person/effective person).\"\"\"\n", " # extract params\n", " n = params['n']\n", " g = params['g']\n", " alpha = params['alpha']\n", " delta = params['delta']\n", " rho = params['rho']\n", " theta = params['theta']\n", " sigma = params['sigma']\n", " gamma = (sigma - 1) / sigma\n", " \n", " # nest Cobb-Douglas as special case\n", " if gamma == 0:\n", " k_star = (alpha / (delta + rho + theta * g))**(1 / (1 - alpha))\n", " else:\n", " k_star = (1 - alpha)**(1 / gamma) * ((alpha / (delta + rho + theta * g))**(gamma / (gamma - 1)) - alpha)**(-1 / gamma)\n", " \n", " return k_star\n", "\n", "def analytic_c_star(params): \n", " \"\"\"Steady-state level of consumption (per person/effective person).\"\"\"\n", " # extract params\n", " n = params['n']\n", " g = params['g']\n", " delta = params['delta']\n", " \n", " # compute k_star\n", " k_star = analytic_k_star(params)\n", " \n", " # compute c_star\n", " c_star = ces_output(0, k_star, params) - (n + g + delta) * k_star\n", " \n", " return c_star" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Choosing baseline parameters:\n", "Finally, to complete the description of the model, we need to provide some \"reasonable\" values for the model parameters in the form of a [Python dictionary](http://www.greenteapress.com/thinkpython/html/thinkpython012.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# start by defining a dictionary of parameter values\n", "baseline_params = {'alpha':0.33, 'n':0.025, 'g':0.025, 'rho':0.04, 'theta':2.5, \n", " 'delta':0.1, 'sigma':1.0, 'A0':1.0, 'L0':1.0}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating an instance of the `ramsey.Model` class:\n", "Now we are ready to create an instance of the `ramsey.Model` class representing the optimal growth model..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# we create an instance of the optimgal growth model as follows...\n", "model = ramsey.Model(output=ces_output, \n", " mpk=marginal_product_capital, \n", " k_dot=equation_motion_capital, \n", " c_dot=consumption_euler_equation,\n", " utility=crra_utility,\n", " jacobian=ramsey_jacobian, \n", " params=baseline_params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adding the steady state expressions to the model:\n", "\n", "Now we add the steady state functions to the model using by passing a dictionary containing the steady state expressions to the `set_functions` method of the model's `steady_state` attribute." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create a dictionary containing the steady state expressions\n", "steady_state_funcs = {'k_star':analytic_k_star, 'c_star':analytic_c_star}\n", "\n", "# pass it as an arg to the set_functions method\n", "model.steady_state.set_functions(func_dict=steady_state_funcs)\n", "\n", "# compute the steady state values!\n", "model.steady_state.set_values()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# display the values as follows\n", "model.steady_state.values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Task 2: Calibration:\n", "\n", "In this task I am going to walk you through a simple strategy for calibrating an optimal growth model with CES production using data from the [Penn World Tables (PWT)](http://www.rug.nl/research/ggdc/data/penn-world-table).\n", "\n", "## Labor force growth rate, $n$:\n", "To obtain a measure of the labor force growth rate for country $i$, I regress the logarithm of total employed persons, $\\verb|emp|$, in country $i$ from the PWT on a constant and a linear time trend.\n", "\n", "$$ \\ln\\ \\verb|emp|_i = \\beta_0 + \\beta_1 \\verb|t| + \\epsilon_i \\tag{2.1}$$\n", "\n", "The estimated coefficient $\\beta_1$ is then used as my estimate for the $n$. To estimate the initial number of employed persons, $L_0$, I use $e^{\\beta_0}$ as the estimate for $L_0$.\n", "\n", "## Capital depreciation rate, $\\delta$: \n", "The PWT provides estimated depreciation rates that vary across both time and countries. As an estimate for the rate of capital depreciation for country $i$ I use the time-average of $\\verb|delta_k|_{it}$ as provided by the PWT. \n", "\n", "## Elasticity of substitution between capital and effective labor, $\\sigma$, and $\\alpha$: \n", "Setting `sigma = 1.0` will calibrate a model with Cobb-Douglas production. In which case $\\alpha$ is estimated using data on a country's capital share. Capital's share for country $i$ in year $t$, $\\alpha_{it}$ is computed as $\\alpha_{it} = 1 - \\verb|labsh|_{it}$, where $\\verb|labsh|_{it}$ is the labor share of output/income for country $i$ in year $t$ provided by the PWT. I then use the time-average of $\\alpha_{it}$ as the estimate of capital's share for country $i$.\n", "\n", "If `sigma != 1.0` then both $\\sigma$ and $\\alpha$ are jointly estimated using non-linear least squares to solve the following optimization problem.\n", "\n", "$$\\min_{\\sigma, \\alpha}\\ \\sum_{t=1}^{T}\\left[\\ln Y_{i,t} - \\left(\\frac{\\sigma}{\\sigma - 1}\\right)\\ln\\bigg(\\alpha K_{i,t}^{\\frac{\\sigma-1}{\\sigma}} + (1 - \\alpha)(A_{i,t}L_{i,t})^{\\frac{\\sigma-1}{\\sigma}}\\bigg)\\right]. \\tag{2.2}$$\n", "\n", "Estimating the parameters of a CES production function is far from trivial:\n", "\n", "1. Objective function has lots of local minima (which makes finding the global minimum hard!).\n", "2. Objective function has a discontinuity at $\\sigma=1.0$.\n", "3. Estimation results may be sensitive to initial guess.\n", "\n", "The estimation procedure...\n", "\n", "## Growth rate of technology, $g$:\n", "To obtain a measure of the growth rate of technology for country $i$, I first adjust the total factor productivity measure reported by the PWT, $\\verb|rtfpna|$ (which excludes the human capital contribution to TFP), and then regress the logarithm of this adjusted measure of TFP, $\\verb|atfpna|$, for country $i$ on a constant and a linear time trend.\n", "\n", "$$\\ln\\ \\verb|atfpna|_i = \\beta_0 + \\beta_1 \\verb|t| + \\epsilon_i \\tag{2.3}$$\n", "\n", "The estimated coefficient $\\beta_1$ is then used as my estimate for the $g$. To estimate the initial level of technology, $A_0$, I use $e^{\\beta_0}$ as the estimate for $A_0$.\n", "\n", "## Elasticity of intertemporal substitution, $\\frac{1}{\\theta}$:\n", "The inverse elasticity of intertemporal substitution is calibrated in so that the model's balanced growth path capital-output ratio matches the average capital-output ratio observed in the data. Put another way, I choose $\\theta$ to solve:\n", "\n", "$$ \\frac{k^*(\\alpha, \\delta, \\rho, \\sigma, \\theta, g)}{y^*(\\alpha, \\delta, \\rho, \\sigma, \\theta, g)} - \\frac{1}{T}\\sum_{t=1}^{T} \\frac{K_{i,t}}{Y_{i,t}} = 0 \\tag{2.4}$$\n", "\n", "## Discount rate, $\\rho$: \n", "The discount rate is a free parameter and must be specified by the user. The calibrated value of $\\theta$ may be sensitive to the choice of $\\rho$.\n", "\n", "All of this is being done \"behind the scenes\". All you need to in order to calibrate the model is pick an [ISO 3 country code](http://en.wikipedia.org/wiki/ISO_3166-1_alpha-3)! Here is an example of the syntax for calibrating an optimal growth model to UK data for the years 1950-2011 (i.e., the entire sample).\n", "\n", " ramsey.calibrate_ces(model=model, iso3_code='GBR', bounds=[1950, 2011], sigma0=0.5, alpha0=0.5, theta0=1.5, rho=0.04)\n", " \n", "The keyword arguments are defined as follows: \n", "\n", "* `model` is an instance `ramsey.Model` class; \n", "* `iso3_code` is a valid ISO 3 country code; \n", "* `bounds` is a list containing the start and end years for the subset of data you wish to use in the calibration; \n", "* `sigma0` is the initial guess of the optimal value of $\\sigma$ (setting `sigma=1.0` will calibrate a model with Cobb-Douglas production); \n", "* `alpha0` is the initial guess for $\\alpha$ (value is ignored if `sigma=1.0`); \n", "* `theta0` is the initial guess of the \"true\" value of $\\theta$;\n", "* `rho` is the user specifed value for the discount rate. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# always check the docstring!\n", "ramsey.calibrate_ces?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# calibrate the model to the UK using only data from 1950-2011\n", "ramsey.calibrate_ces(model, iso3_code='GBR', bounds=[1950, 2011], sigma0=0.5, alpha0=0.5, theta0=2.5, rho=0.04)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# disply the calibrated parameters\n", "print model.params" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# display the steady state values\n", "print model.steady_state.values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Task 3: Analyzing the phase diagram\n", "\n", "From the lectures you should recall that when the current level of capital per effectiver worker is less than its long-run steady state value, then \n", "\n", "$$f'(k) > \\delta + \\rho + \\theta g$$\n", "\n", "and consumption per effective worker is rising. Similarly, when the current level of capital per effectiver worker is greater than its long-run steady state value, then\n", "\n", "$$f'(k) < \\delta + \\rho + \\theta g$$\n", "\n", "and consumption per effective worker is falling. The locus of points for which consumption per effective worker is constant (i.e., the set of points for which $\\dot{c}=0$) is called the $\\dot{c}=0$ locus.\n", "\n", "\n", " that the phase diagram for the optimal growth model combines the locus of points for which consumption per effective worker is constant (i.e., the set of points for which $\\dot{c}=0$) with the locus of points for which capital per effective worker is constant (i.e., the set of points for which $\\dot{k}=0$).\n", "\n", "\n", "\n", "The syntax for plotting a basic phase diagram uses the `plot_phase_diagram` method of the `Model` class as follows:\n", "\n", " model.plot_phase_diagram(gridmax=50, N=1000, arrows = True)\n", " \n", "the `gridmax` argument specifies the maximum value of the grid of values for capital per effective worker to include in the plot, `N` is the number of grid points to include in the plot, `arrows` specifies whether you wish to include directional arrows in the plot. All of this information (and much more!) is contained in the docstring for the `plot_phase_diagram` method." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# always check the docstring of a function before using it!\n", "model.plot_phase_diagram?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# generate a simple phase diagram plot\n", "plt.figure(figsize=(8,6))\n", "\n", "model.plot_phase_diagram(gridmax=150, N=1000, arrows=True)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exploring parameter shocks\n", "\n", "The syntax for plotting the response of the phase diagram to shocks to the various exogenous parameters of the model is as follows:\n", "\n", " model.plot_phase_diagram(gridmax=100, N=1000, arrows = True, param='delta', shock=0.9, reset=True)\n", " \n", "The `gridmax`, `N`, and `arrows` arguments are the as above. The `param` argument is the model parameter that you wish to shock, `shock` specifies the magnitude of the *multiplicative* shock that you wish to apply to the parameter, and `reset` specifies whether you wish to reset the model parameters to their pre-shock values.\n", "\n", "Thw next cell contains demonstrates how to plot the result of a 50% decrease in the growth rate of technology, $g$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# suppose there is a 50% negative shock to g\n", "plt.figure(figsize=(8,6))\n", "model.plot_phase_diagram(gridmax=150, N=1000, arrows=True, param='g', shock=0.5, reset=True)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that the location of the $\\dot{c}=0$ locus is determined by the steady state value of $k$ which solves\n", "\n", "$$f'(k^*) = \\delta + \\rho + \\theta g.$$\n", "\n", "Differentiating the above expression with respect to $g$ yields\n", "\n", "$$f''(k^*)\\frac{\\partial k^*}{\\partial g} = \\theta \\implies \\frac{\\partial k^*}{\\partial g} = \\frac{\\theta}{f''(k^*)} < 0$$\n", "\n", "is negative because $f''(k) < 0$ for all $k$. Therefore a fall in $g$ actually increase $k^*$ and thus shifts the $\\dot{c}=0$ locus to the right. To see why the $\\dot{k}=0$ locus shifts recall that the $\\dot{k}=0$ locus is the set of values for $c$ which lead to actual investment equaling break even investment.\n", "\n", "$$ 0 = f(k) - c - (n+g+\\delta)k \\implies f(k) - c = (n + g + \\delta)k$$\n", "\n", "A fall in $g$ clearly reduces break even invesment. In order for actual investment to fall, consumption per effective worker must rise. Thus a fall in $g$ raises the $\\dot{k}=0$ locus upwards for every level of $k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4: \n", "Describe how each of the following shocks affects both the $\\dot{c}=0$ locus and the $\\dot{k}=0$ locus. What happens the the balanced growth path values of $c$ and $k$?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### part a) \n", "A rise in the elasticity of substitution between consumption today and consumption tomorrow. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# INSERT CODE HERE!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The code in the cell below uses the [IPython magic command](http://ipython.org/ipython-doc/dev/interactive/tutorial.html) `%run` with the `-i` flag to Python script `exercise_4a.py` in IPython's namespace instead of an empty one. This is useful really useful if you are experimenting with code written in a text editor which depends on variables defined interactively. For more information about the `%run` magic command...read the docstring!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# check the docstring!\n", "%run?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution (100% increase in elasticity of substitution)\n", "%run -i exercise_4a.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### part b)\n", "A fall in the elasticity of substitution between capital and effective labor. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# INSERT CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution (10% fall)\n", "%run -i exercise_4b.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### part c)\n", "An increase in the depreciation rate of physical capital." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# INSERT CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution (50% increase)\n", "%run -i exercise_4c.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Task 4: Simulating and solving the model \n", "Below is the familiar phase diagram for optimal growth model. Because the optimal growth model is saddle point stable (I will formally define saddle point stability in a minute) there exist two invariant manifolds: a stable manifold, $M_S$ and an unstable manifold, $M_U$. These manifolds are called invariant because any path that begins on $M_S (M_U)$ will remain on $M_S (M_U)$. $M_S$ is called the stable manifold because any path that begins on the stable manifold will eventually converge to the steady state; $M_U$ is called the unstable manifold because any path the begins on $M_U$ will diverge away from the steady state. In order to solve the optimal growth model we need to compute its stable manifold, $M_S$.\n", "\n", "
\n", " \"Forward\n", "
\n", "\n", "**Note that what Romer calls the \"saddle-path\" is the same thing that I am calling the stable-manifold, $M_S$!**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulating the optimal growth model\n", "Before discussing the various methods for solving the optimal growth model it is necessary to first learn how to simulate trajectories (or time-paths) of the model. Just as we did in the previous lab on the Solow model, we will simulate trajectories for the optimal growth model using the `integrate` method. The following is a basic example of the syntax used to simulate a time path for the economy given any initial condition. \n", "\n", " model.integrate(t0=0, y0=[k0, c0], h=1.0, T=10, integrator='lsoda', **kwargs)\n", " \n", "The arguments are defined as follows: \n", "\n", "* `t0`: is the initial condition for the independent variable (which without loss of generality we can always set to zero for the optimal growth model); \n", "* `y0`: is an array-like container of initial conditions for the endogenous variables of the model (in this case we need to specify *both* $k_0$ and $c_0$; \n", "* `h`: is a floating point number that defines the step-size used by the numerical solver when approximating the time-path of the economy;\n", "* `T`: is an integer representing the desired length of the time-path;\n", "* `integrator`: is a valid finite-difference method for solving systems of ordinary differential equations.\n", "* `**kwargs`: is an optional dictionary of keyword arguments that allow the user to directly tune the integrator as needed. Note that valid keyword arguments are integrator specific. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# always check the docstring!\n", "model.integrate?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# initial conditions for c and k\n", "c0 = model.steady_state.values['c_star'] / 4\n", "k0 = model.steady_state.values['k_star'] / 4\n", "\n", "# combine into a vector (ordering is important!)\n", "init_vec = np.array([k0, c0])\n", "\n", "# simulate a time path of the economy (column ordering is [t, k, c]!)\n", "model.integrate(t0=0, y0=init_vec, h=1.0, T=10, integrator='lsoda')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot a trajectory of the economy in phase space as follows..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=135, N=1000, arrows=True)\n", "\n", "# simulate the model\n", "traj = model.integrate(t0=0, y0=init_vec, h=1.0, T=50, integrator='lsoda')\n", "\n", "# plot the trajectory \n", "model.plot_trajectory(traj, color='r')\n", "\n", "# demarcate the initial condition\n", "plt.vlines(k0, 0, c0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "plt.hlines(c0, 0, k0, linestyle='dashed', color='k')\n", "plt.yticks([c0], ['$c_0$'])\n", "\n", "# change the plot title\n", "plt.title('A time path of the economy in phase space...', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a case of very low initial consumption. Here both $\\dot{c}$ and $\\dot{k}$ are initially positive, but small. From the consumption Euler equation, $\\dot{c}$ is directly proportional to $c$: when $c$ is small, $\\dot{c}$ is also small. Therefore $c$ remains low and the economy evenually crosses the $\\dot{c}=0$ locus. After crossing the $\\dot{c}=0$ locus, $\\dot{c}$ becomes negative while $\\dot{k}$ remains positive. The economy moves down and to the right.\n", "\n", "The code in the cell below fixes the initial condition for capital per effective worker and then plots trajectories for various initial conditions for consumption per effective worker..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=130, N=1000, arrows=False)\n", "\n", "# new initial condition for k\n", "k0 = model.steady_state.values['k_star'] / 2\n", "\n", "# set of initial conditions for c\n", "N = 10\n", "c_lower = model.steady_state.values['c_star'] / 4\n", "c_upper = 2 * model.steady_state.values['c_star']\n", "init_vals = np.linspace(c_lower, c_upper, N)\n", "\n", "# color scheme\n", "color_map = mpl.cm.hot(np.linspace(0, 1, N))\n", "\n", "for i, c0 in enumerate(init_vals):\n", " \n", " # simulate the model\n", " traj = model.integrate(t0=0, y0=[k0, c0], h=0.1, T=20, integrator='lsoda')\n", "\n", " # plot the trajectory \n", " model.plot_trajectory(traj, color=color_map[i])\n", "\n", "# demarcate the initial condition\n", "plt.axvline(k0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "\n", "# change the plot title\n", "plt.title('Time paths of the economy in phase space...', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can easily repeat the above experiment for an initial condition for capital per effective worker which is greater than its steady state value..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=130, N=1000, arrows=False)\n", "\n", "# new initial condition for k\n", "k0 = 2 * model.steady_state.values['k_star']\n", "\n", "# set of initial conditions for c\n", "N = 10\n", "c_lower = model.steady_state.values['c_star'] / 4\n", "c_upper = 2 * model.steady_state.values['c_star']\n", "init_vals = np.linspace(c_lower, c_upper, N)\n", "\n", "# color scheme\n", "color_map = mpl.cm.hot(np.linspace(0, 1, N))\n", "\n", "for i, c0 in enumerate(init_vals):\n", " \n", " # simulate the model\n", " traj = model.integrate(t0=0, y0=[k0, c0], h=0.1, T=20, integrator='lsoda')\n", "\n", " # plot the trajectory \n", " model.plot_trajectory(traj, color=color_map[i])\n", "\n", "# demarcate the initial condition\n", "plt.axvline(k0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "\n", "# change the plot title\n", "plt.title('Time paths of the economy in phase space...', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise:\n", "\n", "Play around with initial condition for consumption per effective worker, $c_0$, and see if you can find the a value for $c_0$ consistent with a trajectory that \"hits\" the steady at $T=50$. If you are successful, try increasing the length of the trajectory. Can you \"hit\" the steady state at $T=100$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=130, N=1000, arrows=False)\n", "\n", "# you need to modify c0 and T to try and \"hit\" the steady state!\n", "k0 = 4 * model.steady_state.values['k_star']\n", "c0 = # INSERT YOUR CODE HERE!\n", "T = # INSERT YOUR CODE HERE!\n", "traj = model.integrate(t0=0, y0=[k0, c0], h=0.1, T=T, integrator='dopri5')\n", "\n", "# plot the trajectory \n", "model.plot_trajectory(traj, color='r')\n", "\n", "# demarcate the initial condition\n", "plt.vlines(k0, 0, c0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "plt.hlines(c0, 0, k0, linestyle='dashed', color='k')\n", "plt.yticks([c0], ['$c_0$'])\n", "\n", "# change the plot title\n", "plt.title('Can you hit the steady state?', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving the model via linearization\n", "A frequently used qualitative technique for analyzing systems of non-linear differential equations is to draw inferences from linear approximations to the system. Such linear approximations are derived from the Taylor expansion of the system around an equilibrium point. If the model has multiple equilibria, each equilibrium requires a *separate* linear approximation. While linear (or higher order polynomial) approximations derived using a Taylor expansion can give the exact value of the function at the point of approximation, the approximation errors may grow rapidly as we move farther away from the point of expansion. Same ideas apply to a linear approximation of arbitrary non-linear system. At the point of expansion (typically some steady-state equilibrium point of the system), the linear approximation pinpoints exactly the same equilibrium of the non-linear system. For a sufficiently small neighborhood of the equilibrium point, the behavior of the linear approximation should be roughly the same as the behavior of the non-linear system. \n", "\n", "To construct a linear approximation of the optimal growth model, begin by noting that time enters only indirectly into the analysis. Thus we can write\n", "\n", "$$\\dot{k}(t) = \\dot{k}(k, c) \\\\ \\dot{c}(t) = \\dot{c}(k, c)$$ \n", "\n", "We now wish to approximate this system of equations using a first-order Taylor expansion around the long-run equilibrium point $(k^*, c^*)$.\n", "\n", "\\begin{align}\n", "\\dot{k}(k, c) \\approx& \\dot{k}(k^*, c^*) + \\frac{\\partial \\dot{k}}{\\partial k}\\bigg|_{k=k^*, c=c^*}(k - k^*) + \\frac{\\partial \\dot{k}}{\\partial c}\\bigg|_{k=k^*, c=c^*}(c - c^*) \\\\\n", "\\dot{c}(k, c) \\approx& \\dot{c}(k^*, c^*) + \\frac{\\partial \\dot{c}}{\\partial k}\\bigg|_{k=k^*, c=c^*}(k - k^*) + \\frac{\\partial \\dot{c}}{\\partial c}\\bigg|_{k=k^*, c=c^*}(c - c^*) \n", "\\end{align}\n", "\n", "Note that $\\dot{k}(k^*, c^*) = \\dot{c}(k^*, c^*) =0$, and define the following change of variables.\n", "\n", "$$\\tilde{k}(t) = k(t) - k^* \\\\ \\tilde{c}(t) = c(t) - c^*$$ \n", "\n", "These definitions imply that $\\dot{\\tilde{k}}(t) = \\dot{k}(t)$ and $\\dot{\\tilde{c}}(t) = \\dot{c}(t)$. Using this change of variables we can re-write the original non-linear system as follows.\n", "\n", "\\begin{align}\n", "\\begin{bmatrix}\n", "\t\\dot{\\tilde{k}}(t) \\\\ \n", "\t\\dot{\\tilde{c}}(t)\n", "\\end{bmatrix}\n", "=&\n", "\\begin{bmatrix}\n", "\t\\frac{\\partial \\dot{k}}{\\partial k} & \\frac{\\partial \\dot{k}}{\\partial c} \\\\\n", "\t\\frac{\\partial \\dot{c}}{\\partial k} & \\frac{\\partial \\dot{c}}{\\partial c} \\\\\n", "\\end{bmatrix}\n", "\\Bigg|_{k=k^*, c=c^*}\n", "\\begin{bmatrix}\n", "\t\\tilde{k}(t) \\\\\n", "\t\\tilde{c}(t)\n", "\\end{bmatrix} \\tag{4.2.1} \\\\\n", "\\tilde{k}(0) =& k_0 - k^* \\tag{4.2.2}\\\\\n", "\\lim_{t\\rightarrow\\infty} |\\tilde{c}(t)| =& 0 \\tag{4.2.3}\n", "\\end{align}\n", "\n", "where $\\tilde{k}(t)$ is the pre-determined (i.e., state) variable and $\\tilde{c}(t)$ is a free (i.e., control) variable. Note that we have transformed our original non-linear boundary value problem into a linear boundary value problem with constant coefficients which, assuming certain determinacy conditions are satisfied, can be solved explicitly.\n", "\n", "The next step in the linearization is to compute the Jacobian matrix of partial derivatives and evaluate it at the long-run equilibrium. \n", "\n", "\\begin{align}\n", "\\mathcal{J}\\big|_{k=k^*, c=c^*}=& \n", "\\begin{bmatrix}\n", "\t\\frac{\\partial \\dot{k}}{\\partial k} & \\frac{\\partial \\dot{k}}{\\partial c} \\\\\n", "\t\\frac{\\partial \\dot{c}}{\\partial k} & \\frac{\\partial \\dot{c}}{\\partial c} \\\\\n", "\\end{bmatrix}\n", "\\Bigg|_{k=k^*, c=c^*}\n", "=\n", "\\begin{bmatrix}\n", "\tf'(k^*) - (n + g + \\delta) & -1 \\\\\n", "\t\\frac{c^*}{\\theta}f''(k^*) & \\frac{f'(k^*) - \\delta - \\rho - \\theta g}{\\theta} \n", "\\end{bmatrix}\n", "=\n", "\\begin{bmatrix}\n", "\t\\beta & -1 \\\\\n", "\t\\frac{c^*}{\\theta}f''(k^*) & 0\n", "\\end{bmatrix}\n", "\\end{align}\n", "\n", "where $\\beta = f'(k^*) - (n + g + \\delta)$. Substituting this result into 4.2.1 yields:\n", "\n", "\\begin{align}\n", "\\begin{bmatrix}\n", "\t\\dot{\\tilde{k}}(t) \\\\ \n", "\t\\dot{\\tilde{c}}(t)\n", "\\end{bmatrix}\n", "=&\n", "\\mathcal{J}\\big|_{k=k^*, c=c^*}\n", "\\begin{bmatrix}\n", "\t\\tilde{k}(t) \\\\\n", "\t\\tilde{c}(t)\n", "\\end{bmatrix} \\tag{4.2.4} \\\\\n", "\\tilde{k}(0) =& k_0 - k^* \\\\\n", "\\lim_{t\\rightarrow\\infty} |\\tilde{c}(t)| =& 0\n", "\\end{align}\n", "\n", "Now we need to solve this system. To solve linear systems of differential equations with constant coefficients we need to diagonalize the matrix $\\mathcal{J}$ by decomposing it into the product of two special matrices \n", "\n", "$$ P = \\begin{bmatrix}P_{00} & P_{01} \\\\ P_{10} & P_{11} \\\\ \\end{bmatrix},\\ \\Lambda=\\begin{bmatrix}\\lambda_{0} & 0 \\\\ 0 & \\lambda_1 \\\\ \\end{bmatrix} \\tag{4.2.5}$$ \n", "\n", "such that \n", "\n", "\\begin{equation}\n", "\\mathcal{J}\\big|_{k=k^*, c=c^*} = P^{-1}\\Lambda P \\tag{4.2.6}\n", "\\end{equation} \n", "\n", "where the matrix $\\Lambda$ is a diagonal matrix of whose entries are the eigenvalues of $\\mathcal{J}$ and the rows of the matrix $P$ are the normalized left-eigenvectors of $\\mathcal{J}$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# store the model's steady state values as a list\n", "steady_state = [model.steady_state.values['k_star'], model.steady_state.values['c_star']]\n", "\n", "# evaluate the jacobian of the system at the models steady state \n", "evaluated_jacobian = ramsey_jacobian(0, steady_state, model.params)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# display the result\n", "print evaluated_jacobian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computing eigenvalues and eigenvectors by hand is an incredibly tedious and error-prone process. Fortunately there exist a large number of high-quality algorithms for computing eigvalues and eigenvectors. We will use the `scipy.linalg.eig` routine. You can learn more about this routine by checking out the docstring..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# always read the docstring! \n", "linalg.eig?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Computing the eigenvalues and eigenvectors can be done as follows..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# linalg.eig returns a tuple of the form (values, vectors)\n", "eig_vals, eig_vecs = linalg.eig(evaluated_jacobian, left=True, right=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# create the diagonal matrix Lambda using the eigenvalues\n", "Lambda = np.diag(eig_vals)\n", "\n", "# create the matrix P where rows are eigenvectors!\n", "P = eig_vecs.T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Confirm that the matrix diagonalization formula that I provided above is in fact correct..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# compute P^-1\n", "P_inverse = linalg.inv(P)\n", "\n", "# compare this result...\n", "print P_inverse.dot(Lambda).dot(P)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# ...with this one!\n", "print evaluated_jacobian" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Need to discuss determinacy issues. In general, systems of linear equations can have zero, one, or an infinite number of solutions. We want 4.2.5 to have a unique solution. In order for the system to have a unique solution, we need...\n", "\n", "1. Number of eigenvalues with *negative* real parts (i.e., stable eigenvalues) equals the number of pre-determined variables;\n", "2. Number of eigenvalues with *positive* real parts (i.e., unstable eigenvalues) equals the number of free or control variables. \n", "\n", "The optimal growth model has one pre-determined variable (i.e., capital per effective worker) and one free or control variable (i.e., consumption per effective worker). Thus in order for the model to have a unique solution we need one eigenvalue with negative real part and one eigenvalue with positive real part.\n", "\n", "Let's check the eigenvalues to make sure that we have a unique solution!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# mathematically speaking, the saddle path is unique!\n", "eig_vals" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Without loss of generality we can re-order the rows and columns of $P$ and $\\Lambda$ so that the first row of $P$ corresponds to the stable eigenvector and $\\lambda_0$ is the stable eigenvalue." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# swap rows of P \n", "P = P[::-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# note first row contains the stable eigenvector!\n", "P" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# swap rows and colums of Lambda\n", "Lambda = Lambda[::-1, ::-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# now lambda_0 has negative real part\n", "Lambda" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why can we swap rows and columns? [Elementary row/column operations](http://en.wikipedia.org/wiki/Elementary_row_operations#Operations) don't impact diagonalization procedure!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# compute P^-1\n", "P_inverse = linalg.inv(P)\n", "\n", "# still equals Jacobian!\n", "print P_inverse.dot(Lambda).dot(P)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The whole point of diagonalizing the Jacobian is that doing so allows us to re-write equation 4.2.5 as follows.\n", "\n", "\\begin{align}\n", "\\begin{bmatrix}\n", "\t\\dot{\\tilde{k}}(t) \\\\ \n", "\t\\dot{\\tilde{c}}(t)\n", "\\end{bmatrix}\n", "=&\n", "P^{-1}\\Lambda P\n", "\\begin{bmatrix}\n", "\t\\tilde{k}(t) \\\\\n", "\t\\tilde{c}(t)\n", "\\end{bmatrix} \\tag{4.2.8}\n", "\\end{align}\n", "\n", "Pre-multiplication by $P$ yields...\n", "\n", "\\begin{align}\n", "\\left(P\n", "\\begin{bmatrix}\n", "\t\\dot{\\tilde{k}}(t) \\\\ \n", "\t\\dot{\\tilde{c}}(t)\n", "\\end{bmatrix}\\right)\n", "=\n", "\\Lambda \\left(P\n", "\\begin{bmatrix}\n", "\t\\tilde{k}(t) \\\\\n", "\t\\tilde{c}(t)\n", "\\end{bmatrix}\\right) \\tag{4.2.9}\n", "\\end{align}\n", "\n", "Now we define *another* change of variables...\n", "\n", "$$\n", "\\begin{bmatrix}\n", "\t\\hat{k} \\\\ \n", "\t\\hat{c}\n", "\\end{bmatrix}\n", "=\n", "P\n", "\\begin{bmatrix}\n", "\t\\tilde{k} \\\\ \n", "\t\\tilde{c}\n", "\\end{bmatrix} \\tag{4.2.10}\n", "$$\n", "\n", "...in order to further reduce 4.2.5 to the following *separable* system of differential equations...\n", "\n", "\\begin{align}\n", "\\begin{bmatrix}\n", "\t\\dot{\\hat{k}}(t) \\\\ \n", "\t\\dot{\\hat{c}}(t)\n", "\\end{bmatrix}\n", "=&\n", "\\Lambda\n", "\\begin{bmatrix}\n", "\t\\hat{k}(t) \\\\\n", "\t\\hat{c}(t)\n", "\\end{bmatrix} \\tag{4.2.11}\\\\\n", "\\hat{k}(0) =& \\hat{k}(0) \\\\\n", "\\lim_{t\\rightarrow\\infty} |\\hat{c}(t)| =& 0 \n", "\\end{align}\n", "\n", "Why is system described by 4.2.11 *separable*? The differential equation governing the evolution of $\\hat{k}$ doesn't depend on the differential equation governing the evolution of $\\hat{c}$! This means that each of the differential equations can be solved separately:\n", "\n", "\\begin{align}\n", " \\hat{k}(t) =& \\hat{k}(0)e^{\\lambda_0 t} \\tag{4.2.12} \\\\\n", " \\hat{c}(t) =& \\hat{c}(0)e^{\\lambda_1 t} \\tag{4.2.13} \\\\\n", "\\end{align}\n", "\n", "where\n", "\n", "\\begin{align}\n", " \\hat{k}(0) =& [P_{00}\\tilde{k}(0) + P_{01}\\tilde{c}(0)] \\\\\n", " \\hat{c}(0) =& [P_{10}\\tilde{k}(0) + P_{11}\\tilde{c}(0)] \\\\\n", "\\end{align}\n", "\n", "Now $\\lambda_1 > 0$ implies that $\\hat{c}(t)$ is exploding as $t\\rightarrow \\infty$! But recall that we need to satisfy $\\lim_{t\\rightarrow\\infty} |\\hat{c}(t)|$. Solution is to choose $\\tilde{c}(0)$ such that $\\hat{c}(0)=0$ which forces $\\hat{c}(t)=0$.\n", "\n", "$$ P_{10}\\tilde{k}(0) + P_{11}\\tilde{c}(0) = 0 \\implies \\tilde{c}(0) = -\\frac{P_{10}}{P_{11}}\\tilde{k}(0) \\implies \\tilde{c}(t) = -\\frac{P_{10}}{P_{11}}\\tilde{k}(t) \\tag{4.2.14} $$\n", "\n", "Substituting our result for $\\tilde{c}_0$ (and a bit of algebra!) allows us to recover the solution to 4.2.5.\n", "\n", "\\begin{align}\n", " \\tilde{k}(t) =& \\tilde{k}(0)e^{\\lambda_0 t} \\tag{4.2.15}\\\\\n", " \\tilde{c}(t) =& -\\frac{P_{10}}{P_{11}}\\tilde{k}(t) \\tag{4.2.16}\\\\\n", "\\end{align}\n", "\n", "Finally, reversing our original change of variables yields the desired linear approximations of the functions $k(t)$ and $c(t)$ that solve the original system of differential equations.\n", "\n", "\\begin{align}\n", " k(t) =& k^* + (k(0) - k^*)e^{\\lambda_0 t} \\tag{4.2.17} \\\\\n", " c(t) =& c^* - \\frac{P_{10}}{P_{11}}(k(t) - k^*) \\tag{4.2.18} \\\\\n", "\\end{align}\n", "\n", "Although linearization is really a *qualitative* technique for analyzing it is also the most widely used *quantitative*, numerical method in all of economics. Anyone interested in pursuing economic research (particularly macroeconomic research) should have a solid understanding of this linearization (in particular where it can go horribly wrong!) \n", "\n", "The syntax for computing a linear approximation of the saddle path/stable manifold is as follows:\n", "\n", " model.solve_linearization(k0=k0, ti=t)\n", " \n", "You need to provide an initial condition for capital per effectiver worker, `k0`, as well as a grid of time points, `ti` at which you want to compute the solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# specify an initial condition and a set of time points\n", "k0 = model.steady_state.values['k_star'] / 2\n", "t = np.linspace(0, 100, 10)\n", "\n", "# generate a linearized approximation to the lower branch of the stable manifold\n", "ms_lower_linear = model.solve_linearization(k0, t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# display the linear approximation\n", "print ms_lower_linear" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot a phase diagram\n", "model.plot_phase_diagram(50)\n", "\n", "# plot the lower half of the linearized stable manifold\n", "model.plot_trajectory(ms_lower_linear, color='r')\n", "\n", "# demarcate the initial condition\n", "c0 = ms_lower_linear[0,2]\n", "plt.vlines(k0, 0, c0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "plt.hlines(c0, 0, k0, linestyle='dashed', color='k')\n", "plt.yticks([c0], ['$c_0$'])\n", "\n", "# change the title\n", "plt.title('Linear approximation of saddle-path/stable manifold', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise:\n", "Compute and plot a linear approximation to the upper branch of the saddle-path/stable manifold." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# specify an initial condition and a set of time points\n", "k0 = # INSERT YOUR CODE HERE!\n", "t = # INSERT YOUR CODE HERE!\n", "\n", "# generate a linearized approximation to the lower branch of the stable manifold\n", "ms_upper_linear = # INSERT YOUR CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# specify an initial condition and a set of time points\n", "k0 = 2 * model.steady_state.values['k_star']\n", "t = np.linspace(0, 100, 10)\n", "\n", "# generate a linearized approximation to the lower branch of the stable manifold\n", "ms_upper_linear = model.solve_linearization(k0, t)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot a phase diagram\n", "model.plot_phase_diagram(50)\n", "\n", "# plot the lower half of the linearized stable manifold\n", "model.plot_trajectory(ms_lower_linear, color='r')\n", "\n", "# plot the lower half of the linearized stable manifold\n", "# INSERT YOUR CODE HERE!\n", "\n", "# change the title\n", "plt.title('Linear approximation of saddle-path/stable manifold', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot a phase diagram\n", "model.plot_phase_diagram(50)\n", "\n", "# plot the lower half of the linearized stable manifold\n", "model.plot_trajectory(ms_lower_linear, color='r')\n", "\n", "# plot the lower half of the linearized stable manifold\n", "model.plot_trajectory(ms_upper_linear, color='r')\n", "\n", "# change the title\n", "plt.title('Linear approximation of saddle-path/stable manifold', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving the model via forward-shooting\n", "In the previous section, you found a linear approximation to the stable manifold $M_S$ that was only valid in a neighborhood of the steady state. In this task, we will see how to solve for the full non-linear stable manifold using a \"forward shooting\" approach. The forward shooting method for finding $M_S$ proceeds by guessing a feasible value for $c_0$ and then using some IVP scheme to generate the implied solution trajectories for the variables $k(t)$ and $c(t)$. If the initial choice of $c_0$ is too small, then the solution trajectory eventually crosses the $\\dot{c}=0$ locus and $c(t)$ begins to fall. Similarly, if the choice of $c_0$ is too large, then our path eventually crosses the $\\dot{k}=0$ curve at which point $k(t)$ will start to fall. These observations motivate the forward shooting scheme described in pseudo-code below for an initial condition $k(0)=k_0 < k^*$\n", "\n", "The pseudo-code below basically translates the above logic into something closer to an algorithm that we could implement in order to solve for an approximation of the stable manifold.\n", "\n", " # Bracket the true c(0) by setting cL=0 and cH = y + (1 - delta)k. \n", " cL = 0\n", " cH = f(k0) + (1 - delta) * k0\n", " \n", " # Guess that c0 = 0.5(cH + cL) and specify a stopping tolerance, tol > 0.\n", " c0 = 0.5 * (cH + cL)\n", " \n", " while True:\n", " Solve the model as an IVP with k(0)=k0 and c(0)=c0\n", " if c is falling:\n", " # check for convergence\n", " if |c(T) - c^*| < tol:\n", " break\n", " # current guess too low!\n", " else:\n", " cL = c0 \n", " c0 = 0.5 * (c_H + c_L)\n", "\t\n", " elif k is falling:\n", " # check for convergence\n", " if |c(T) - c^*| < tol:\n", " break \n", " # current guess too high!\n", " else: \n", " cH = c0\n", " c0 = 0.5 * (cH + cL)\n", "\n", "The figure below gives a sense of how the \"forward shooting\" algorithm converges to an initial condition, $c_0^*$, that yields a good approximation of the stable manifold, $M_S$. \n", "\n", "
\n", " \"Forward\n", "
\n", "\n", "To solve the model using forward shooting we use the `solve_forward_shooting` method of the `Model` class. The basic syntax looks as follows:\n", "\n", " model.solve_forward_shooting(k0, h=1.0, tol=1e-6, mesg=True, integrator='dopri5', **kwargs)\n", " \n", "where the arguments are an initial condition for capital per effective worker, `k0`; a step size, `h`, to use when solving the IVP for given initial conditions using some finite-difference method defined by `integrator`; a convergence tolerance, `tol`, defining how close to the steady state is close enough; a boolean flag, `mesg`, specifying whether or not you wish to print progress messages; and, `**kwargs` a dictionary of optional keyword arguments for the IVP solver. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# check the docstring\n", "model.solve_forward_shooting?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set initial condition for capital per effective worker\n", "k0 = model.steady_state.values['k_star'] / 2\n", "\n", "# compute the lower branch of the stable manifold using forward shooting\n", "ms_lower = model.solve_forward_shooting(k0, h=1.0, tol=1.5e-4, mesg=True, integrator='dopri5')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the approximation of the stable manifold in phase space as follows..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=50, N=1000, arrows=False)\n", "\n", "# plot the trajectory \n", "model.plot_trajectory(ms_lower, color='r')\n", "\n", "# demarcate the initial condition\n", "c0 = ms_lower[0,2]\n", "plt.vlines(k0, 0, c0, linestyle='dashed', color='k')\n", "plt.xticks([k0], ['$k_0$'])\n", "plt.hlines(c0, 0, k0, linestyle='dashed', color='k')\n", "plt.yticks([c0], ['$c_0$'])\n", "\n", "# change the plot title\n", "plt.title('\"Forward shooting\" approximation of $M_S$', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise:\n", "\n", "Compute and plot the upper half of the $M_S$ starting from an initial condition for capital per effective worker that is 2 times larger than the steady state value." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set initial condition for capital per effective worker\n", "# INSERT YOUR CODE HERE!\n", "\n", "# compute the upper branch of the stable manifold using forward shooting\n", "# INSERT YOUR CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "# INSERT YOUR CODE HERE!\n", "\n", "# plot the upper branch of the stable manifold\n", "# INSERT YOUR CODE HERE!\n", "\n", "# plot the lower branch of the stable manifold\n", "model.plot_trajectory(ms_lower, color='r')\n", "\n", "# change the plot title\n", "plt.title('You have computed your first stable manifold!', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set initial condition for capital per effective worker\n", "k0 = 2 * model.steady_state.values['k_star']\n", "\n", "# compute the upper branch of the stable manifold using forward shooting\n", "ms_upper = model.solve_forward_shooting(k0, h=1.0, tol=1e-6, mesg=False, integrator='dopri5')\n", "\n", "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=50, N=1000, arrows=False)\n", "\n", "# plot the upper branch of the stable manifold\n", "model.plot_trajectory(ms_upper, color='r')\n", "\n", "# plot the lower branch of the stable manifold\n", "model.plot_trajectory(ms_lower, color='r')\n", "\n", "# change the plot title\n", "plt.title('\"Forward shooting\" approximation of $M_S$', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quick graphical comparison of the linear approximation to the full, non-linear saddle path shows that the linear approximation is tangent to the non-linear approximation at the steady-state equilibrium point." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=50, N=1000, arrows=False)\n", "\n", "# forward shooting approximation of stable manifold\n", "model.plot_trajectory(ms_upper_linear, color='purple')\n", "model.plot_trajectory(ms_lower_linear, color='purple', label='Linearization')\n", "\n", "# forward shooting approximation of stable manifold\n", "model.plot_trajectory(ms_upper, color='r')\n", "model.plot_trajectory(ms_lower, color='r', label='Forward shooting')\n", "\n", "# change the plot title\n", "plt.title('Linearization vs. \"forward-shooting\" approximations', fontsize=20, family='serif')\n", "plt.legend(loc=0, frameon=False, prop={'family':'serif'})\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving the model via reverse-shooting\n", "\n", "The difficulty in applying forward shooting methods to continuous-time, infinite-horizon problems is that the terminal state is extremely sensitive to the initial guess. However, this observation also implies that the initial state corresponding to any given terminal state should be relatively *insensitive* to the terminal state. Thus better strategy for solving continuous-time, infinite-horizon problems is to guess a value for the terminal state and integrate the system of differential equations backwards to find the correspondign initial condition. This method is known as \"reverse shooting.\" \n", "\n", "The phase diagram suggests that using forward shooting to compute the stable manifold, $M_S$, of the model might be difficult for at least some parameter values and initial conditions $k(0)$. Any initial deviation from $M_S$, however small, is magnified over time resulting in a path that increasingly departs from the exact solution. However, suppose that instead of computing the stable manifold, I wished instead to compute the unstable manifold, $M_U$. As the figure below suggests, deviations from $M_U$, however large, become smaller over time. \n", "\n", "
\n", "\n", "In fact, all I would need to do in order to compute a path that lies on the unstable manifold is to choose a point \"near\" the steady state and then integrate the system forward. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=130, N=1000, arrows=True)\n", "\n", "# compute the unstable manifold\n", "eps = 1e-5\n", "k_star, c_star = model.steady_state.values['k_star'], model.steady_state.values['c_star']\n", "\n", "mu_lower = model.integrate(t0=0, y0=[k_star + eps, c_star], h=0.1, T=300, integrator='dopri5')\n", "mu_upper = model.integrate(t0=0, y0=[k_star, c_star + eps], h=0.1, T=300, integrator='dopri5')\n", "\n", "# plot the unstable manifold\n", "model.plot_trajectory(mu_upper, color='r')\n", "model.plot_trajectory(mu_lower, color='r')\n", "\n", "# change the plot title\n", "plt.title('Computing $M_U$ is trivial!', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The \"reverse shooting\" approach transforms the system of equations in such a way so that the stable manifold of the original system becomes the unstable manifold of the transformed system and then solves for the unstable manifold by of the transformed system by shooting \"backwards\" from the steady state. \n", "\n", "Since time plays no direct role in the model, the household's choice of consumption per effective worker at time $t$ depends only on the value of capital per effective worker at time $t$. I can express this by defining a policy function, $c(k)$, such that $c(k)=c(k(t))=c(t)$. Furthermore, the original system of differential equations implies that the consumption policy function satisfies the following differential equation.\n", "\n", "\\begin{equation}\n", " \tc'(k) = \\frac{\\dot{c}}{\\dot{k}} = \\left(\\frac{c(k)}{\\theta}\\right)\\left(\\frac{f'(k) - \\delta - \\rho - \\theta g}{f(k) - (n+g+\\delta)k - c(k)}\\right) \\tag{4.1}\n", "\\end{equation}\n", "\n", "Since optimality requires the economy to converge to its steady state, the solution $c(k)$ must satisfy the boundary condition $c(k^*)=c^*$. The reverse shooting approach solves for the lower portion of the consumption policy function by choosing some initial step-size $\\epsilon > 0$, setting \n", "\n", "\\begin{align}\n", "k_0 =& k^*- \\epsilon \\notag \\\\\n", "c_0 =& c(k^* - \\epsilon) \\approx c^* - \\epsilon c'(k^*) \\notag\n", "\\end{align}\n", "\n", "and then integrating equation 4.1 *backward* using some IVP scheme. To compute the upper portion of the policy function using reverse shooting, simply integrate equation 4.1 *forward* from the initial condition\n", "\n", "\\begin{align}\n", "k_0 =& k^*+\\epsilon \\notag \\\\\n", "c_0 =& c(k^* + \\epsilon) \\approx c^* + \\epsilon c'(k^*)\\notag.\n", "\\end{align} \n", "\n", "The following is an example of the basic syntax used to approximate the saddle-path/stable manifold using reverse shooting.\n", "\n", " model.solve_reverse_shooting(k0=k0, h=1.0, eps=1e-5, integrator='dopri5')\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# as always check the docstring!\n", "model.solve_reverse_shooting?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# compute the upper branch of the stable manifold using reverse shooting\n", "k0 = 2 * k_star\n", "ms_upper_reverse = model.solve_reverse_shooting(k0=k0, h=1.0, eps=1e-5, integrator='dopri5')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# display the trajectory\n", "print ms_upper_reverse" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "plt.figure(figsize=(8,6))\n", "\n", "# plot the phase diagram\n", "model.plot_phase_diagram(gridmax=130, N=1000, arrows=False)\n", "\n", "# compute the stable manifold using reverse shooting\n", "k0 = k_star / 4\n", "ms_lower_reverse = model.solve_reverse_shooting(k0, h=1.0, eps=1e-5, integrator='dopri5')\n", "\n", "k0 = 4 * k_star\n", "ms_upper_reverse = model.solve_reverse_shooting(k0, h=1.0, eps=1e-5, integrator='dopri5')\n", "\n", "# plot the stable manifold\n", "plt.plot(ms_lower_reverse[:,0], ms_lower_reverse[:,1], color='r')\n", "plt.plot(ms_upper_reverse[:,0], ms_upper_reverse[:,1], color='r')\n", "\n", "# change the plot title\n", "plt.title('\"Reverse shooting\" approximation of $c(k)$', fontsize=20, family='serif')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Comparison: Speed vs. Accuracy\n", "\n", "You have now seen three approaches to approximating the saddle-path/stable manifold of the optimal growth model: linearization, forward shooting, and reverse shooting. In this section, we compare the methods in terms of their relative speed and accuracy.\n", "\n", "### Speed" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# calibrate a model with Cobb-Douglas production using US data \n", "ramsey.calibrate_ces(model, iso3_code='USA', bounds=[1950, 2011], sigma0=1.0, rho=0.04)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# define some initial condition for k\n", "k_star = model.steady_state.values['k_star']\n", "\n", "# define some bounds over which to approximate the consumption policy\n", "kmin, kmax = k_star / 4, 4 * k_star" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# time reverse shooting algorithm (best of 3)\n", "%timeit -n 1 -r 3 model.get_consumption_policy(kmin, kmax, method='forward_shooting', h=0.1, tol=1.5e-4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# time reverse shooting algorithm (best of 3)\n", "%timeit -n 1 -r 3 model.get_consumption_policy(kmin, kmax, method='reverse_shooting', h=0.1, eps=1.5e-4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# time linearization algorithm (best of 3)\n", "%timeit -n 1 -r 3 model.get_consumption_policy(method='linearization')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results illustrate a general tradeoff: linearization is typically much faster than reverse shooting, and reverse shooting is typically much faster than forward shooting. Linearization is fast because it is only computing a local approximation of the consumption policy function in a neighborhood of the steady state. Both the forward and reverse shooting methods, meanwhile, seek to globally approximate the consumption policy function over some interval. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Accuracy:\n", "\n", "To assess the relative accuracy of our numerical methods we would ideally like to compare their approximations against an analytical benchmark. Luckily for us it so happens that the optimal growth model with Cobb-Douglas production does have a few closed form solutions. \n", "\n", "The most general closed form solution requires that $\\theta=\\alpha$ and leads to a consumption policy function that is linear in $k$.\n", "\n", "\\begin{align}\n", "\tk(t)=& \\left[k_0^{1-\\alpha}e^{-(1-\\alpha)\\left(\\frac{\\delta + \\rho + \\alpha g}{\\alpha}\\right) t} + \\left(\\frac{\\alpha}{\\delta + \\rho + \\alpha g}\\right)\\left(1 - e^{-(1-\\alpha)\\left(\\frac{\\delta + \\rho + \\alpha g}{\\alpha}\\right) t}\\right)\\right]^{\\frac{1}{1-\\alpha}} \\tag{??}\\\\\n", "\tc(t) =& \\left(\\frac{(1 - \\alpha) \\delta + \\rho - \\alpha n}{\\alpha}\\right)k(t) \\tag{??}\n", "\\end{align}\n", "\n", "Here are some Python functions encoding this analytic solution..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def analytic_capital_1(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for the time path of capital (per person/\n", " effective person) in an optimal growth model with Cobb-Douglas production technology \n", " technology where the inverse elasticity of intertemporal substitution equals capital's \n", " share, alpha.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " \n", " params: (dict) Dictionary of model parameters.\n", " \n", " Returns:\n", " \n", " analytic_k: (array-like) (T,) array representing the time path of \n", " capital (per person/effective person).\n", "\n", " \"\"\"\n", " # extract parameter values\n", " g = params['g']\n", " n = params['n']\n", " alpha = params['alpha']\n", " delta = params['delta']\n", " rho = params['rho']\n", " \n", " # lambda governs the speed of convergence\n", " lmbda = ((delta + rho + alpha * g) / alpha) * (1 - alpha)\n", " \n", " # analytic solution for k at time t\n", " k_t = (((alpha / (delta + rho + alpha * g)) * (1 - np.exp(-lmbda * t)) + \n", " k0**(1 - alpha) * np.exp(-lmbda * t))**(1 / (1 - alpha)))\n", " \n", " return k_t\n", "\n", "def analytic_consumption_1(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for the time path of consumption (per person/\n", " effective person) in a optimal growth model with Cobb-Douglas production technology \n", " technology where the inverse elasticity of intertemporal substitution equals capital's \n", " share, alpha.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " \n", " params: (dict) Dictionary of model parameter values.\n", " \n", " Returns:\n", " \n", " analytic_c: (array-like) (T,) array representing the time path of \n", " consumption (per person/effective person).\n", "\n", " \"\"\"\n", " # extract parameter values\n", " g = params['g']\n", " n = params['n']\n", " alpha = params['alpha']\n", " delta = params['delta']\n", " rho = params['rho']\n", " \n", " # analytic solution consumption at time t\n", " c_t = ((((1 - alpha) * delta + rho - alpha * n) / alpha) * \n", " ramsey_analytic_capital(k0, t, params))\n", " \n", " return c_t\n", "\n", "def analytic_solution_1(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for an optimal growth model with Cobb-Douglas \n", " production technology technology where the inverse elasticity of intertemporal \n", " substitution equals capital's share, alpha.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " \n", " params: (dict) Dictionary of model parameters.\n", " \n", " Returns:\n", " \n", " analytic_traj: (array-like) (T,2) array representing the analytic \n", " solution trajectory.\n", "\n", " \"\"\"\n", " # analytic solution for capital at time t\n", " k_t = analytic_capital_1(k0, t, params) \n", " \n", " # analytic solution for consumption at time t\n", " c_t = analytic_consumption_1(k0, t, params)\n", "\n", " # combine into a (T, 2) array\n", " tup = (t[:,np.newaxis], k_t[:,np.newaxis], c_t[:,np.newaxis])\n", " analytic_traj = np.hstack(tup)\n", " \n", " return analytic_traj" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def analytic_consumption_policy_1(k):\n", " \"\"\"\n", " Computes the analytic solution for the consumption (per person/\n", " effective person) policy function in a Ramsey model with Cobb-Douglas \n", " production technology where the coefficient of relative risk aversion, \n", " theta, equals capital's share, alpha.\n", " \n", " Arguments:\n", " \n", " k: (array-like) Array of values for capital (per person/\n", " effective person) \n", " \n", " Returns:\n", " \n", " c_policy: (array-like) (T,) array representing the consumption \n", " (per person/effective person) policy function.\n", "\n", " \"\"\"\n", " # extract parameter values\n", " g = model.params['g']\n", " n = model.params['n']\n", " alpha = model.params['alpha']\n", " delta = model.params['delta']\n", " rho = model.params['rho']\n", " \n", " # analytic solution consumption at time t\n", " c_policy = (((1 - alpha) * delta + rho - alpha * n) / alpha) * k\n", " \n", " return c_policy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "analytic_params_1 = {'A0':1.0, 'rho':0.04, 'g':0.025, 'theta':0.33, 'L0':1.0, \n", " 'alpha':0.33, 'delta':0.1, 'sigma':1.0, 'n':0.025}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# update the model parameters\n", "model.update_model_parameters(analytic_params_1)\n", "\n", "# re-compute the steady state values\n", "model.steady_state.set_values()\n", "\n", "# extract steady state value\n", "k_star = model.steady_state.values['k_star']\n", "kmin, kmax = 0.5 * k_star, 2.0 * k_star" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute approximations of the consumption policy function using linearization, forward shooting, and reverse shooting..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# approximate consumption policy function using linearization\n", "linearized_policy = model.get_consumption_policy(method='linearization')\n", "\n", "# approximate consumption policy function using forward shooting\n", "forward_shooting_policy = model.get_consumption_policy(kmin, kmax, method='forward_shooting', deg=3)\n", "\n", "# approximate the consumption policy function using reverse shooting\n", "reverse_shooting_policy = model.get_consumption_policy(kmin, kmax, method='reverse_shooting', deg=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...and compare then to the analytic consumption policy function using the $L^2$ metric (i.e., sum of point-wise squared deviations)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# grid of points at which to compare methods\n", "kgrid = np.linspace(kmin, kmax, 1000)\n", "\n", "# linearization is accurate when consumption policy is linear\n", "model.compare_policies(analytic_consumption_policy_1, linearized_policy, kgrid, metric='L2')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# forward shooting looks pretty good\n", "model.compare_policies(analytic_consumption_policy_1, forward_shooting_policy, kgrid, metric='L2')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# reverse shooting also does extremely well\n", "model.compare_policies(analytic_consumption_policy_1, reverse_shooting_policy, kgrid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A linear consumption policy function is not a terribly useful benchmark for differentiating our numerical methods. In order to obtain a closed-form solution of the model that is sufficiently non-linear, I focus on the a special case of the model with a constant gross savings rate. Note that a constant gross savings rate implies that the representative household is saving a constant fraction of its income at each point in time. In order for it to be optimal for the representative household to have a constant gross savings rate, the discount rate must satisfy the following rather peculiar restriction.\n", "\n", "$$ \\rho = \\alpha\\theta(n + g + \\delta) - (\\delta + \\theta g) > 0$$\n", "\n", "Assuming that the above parameter restriction is satisfied, one can obtain the following closed form solution for the functions $k(t)$ and $c(t)$.\n", "\n", "\\begin{align}\n", "k(t)=& \\left[k_0^{1-\\alpha}e^{-\\lambda t} + \\left(\\frac{1}{\\theta(n+g+\\delta)}\\right)\\left(1 - e^{-\\lambda t}\\right)\\right]^{\\frac{1}{1-\\alpha}} \\\\\n", "c(t) =& \\left(\\frac{\\theta - 1}{\\theta}\\right)k(t)^{\\alpha}\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def analytic_capital_2(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for the time path of capital (per person/\n", " effective person) in a Ramsey model with Cobb-Douglas production technology \n", " with a constant gross savings rate.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " params: (dict) Dictionary of model parameters.\n", " \n", " Returns:\n", " \n", " analytic_k: (array-like) (T,) array representing the time path of \n", " capital (per person/effective person).\n", "\n", " \"\"\"\n", " # extract parameter values\n", " g = params['g']\n", " n = params['n']\n", " alpha = params['alpha']\n", " delta = params['delta']\n", " theta = params['theta']\n", " \n", " # lambda governs the speed of convergence\n", " lmbda = (1 - alpha) * (n + g + delta)\n", " \n", " # analytic solution for k at time t\n", " k_t = (((1 / (theta * (n + g + delta))) * (1 - np.exp(-lmbda * t)) + \n", " k0**(1 - alpha) * np.exp(-lmbda * t))**(1 / (1 - alpha)))\n", " \n", " return k_t\n", "\n", "def analytic_consumption_2(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for the time path of consumption (per person/\n", " effective person) in a Ramsey model with Cobb-Douglas production technology \n", " with a constant gross savings rate.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " params: (dict) Dictionary of model parameter values.\n", " \n", " Returns:\n", " \n", " analytic_c: (array-like) (T,) array representing the time path of \n", " consumption (per person/effective person).\n", "\n", " \"\"\"\n", " # extract parameter values\n", " alpha = params['alpha']\n", " theta = params['theta']\n", " \n", " # analytic solution consumption at time t\n", " c_t = ((theta - 1) / theta) * ramsey_analytic_capital(k0, t, params)**alpha\n", " \n", " return c_t\n", "\n", "def analytic_solution_2(k0, t, params):\n", " \"\"\"\n", " Computes the analytic solution for the Ramsey model with Cobb-Douglas\n", " production technology and a constant gross savings rate.\n", " \n", " Arguments:\n", " \n", " k0: (float) Initial value for capital (per person/effective person) \n", " t: (array-like) (T,) array of points at which the solution is \n", " desired.\n", " params: (dict) Dictionary of model parameters.\n", " \n", " Returns:\n", " \n", " analytic_traj: (array-like) (T,2) array representing the analytic \n", " solution trajectory.\n", "\n", " \"\"\"\n", " # analytic solution for capital at time t\n", " k_t = ramsey_analytic_capital(k0, t, params) \n", " \n", " # analytic solution for consumption at time t\n", " c_t = ramsey_analytic_consumption(k0, t, params)\n", "\n", " # combine into a (T, 2) array\n", " tup = (t[:,np.newaxis], k_t[:,np.newaxis], c_t[:,np.newaxis])\n", " analytic_traj = np.hstack(tup)\n", " \n", " return analytic_traj" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def analytic_consumption_policy_2(k):\n", " \"\"\"\n", " Computes the analytic solution for the consumption (per person/\n", " effective person) policy function in a Ramsey model with Cobb-Douglas \n", " production technology and a constant gross savings rate.\n", " \n", " Arguments:\n", " \n", " k: (array-like) Array of values for capital (per person/\n", " effective person) \n", " params: (dict) Dictionary of model parameter values.\n", " \n", " Returns:\n", " \n", " c_policy: (array-like) (T,) array representing the consumption \n", " (per person/effective person) policy function.\n", "\n", " \"\"\"\n", " # extract parameter values\n", " alpha = model.params['alpha']\n", " theta = model.params['theta']\n", " \n", " # analytic solution consumption at time t\n", " c_policy = ((theta - 1) / theta) * k**alpha\n", " \n", " return c_policy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def get_rho(params):\n", " \"\"\"\n", " Computes the required discount rate in order to insure that \n", " optimal consumption policy is consistent with the assumption of\n", " a constant gross savings rate.\n", " \n", " \"\"\"\n", " # extract parameter values\n", " g = params['g']\n", " n = params['n']\n", " alpha = params['alpha']\n", " delta = params['delta']\n", " theta = params['theta']\n", " \n", " rho = alpha * theta * (n + g + delta) - (delta + theta * g)\n", " params['rho'] = rho" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# new set of baseline params...note rho is missing!\n", "analytic_params_2 = {'A0':1.0, 'g':0.02, 'theta':2.5, 'L0':1.0, \n", " 'alpha':0.5, 'delta':0.04, 'sigma':1.0, 'n':0.025}\n", "\n", "# compute rho based on the other params\n", "get_rho(analytic_params_2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# checkout the required discount factor...peculiar restriction doesn't seem that crazy!\n", "analytic_params_2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# update the model parameters\n", "model.update_model_parameters(analytic_params_2)\n", "\n", "# re-compute the steady state values\n", "model.steady_state.set_values()\n", "\n", "# extract steady state value\n", "k_star = model.steady_state.values['k_star']\n", "\n", "kmin, kmax = 0.5 * k_star, 2.0 * k_star\n", "kgrid = np.linspace(kmin, kmax, 1000)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# approximate consumption policy function using linearization\n", "linearized_policy = model.get_consumption_policy(method='linearization')\n", "\n", "# approximate consumption policy function using forward shooting\n", "forward_shooting_policy = model.get_consumption_policy(kmin, kmax, method='forward_shooting', deg=3)\n", "\n", "# approximate the consumption policy function using reverse shooting\n", "reverse_shooting_policy = model.get_consumption_policy(kmin, kmax, method='reverse_shooting', deg=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Results? Linearization is pretty terrible *relative* to..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.compare_policies(analytic_consumption_policy_2, linearized_policy, kgrid, metric='L2')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...forward shooting which is roughly 11 orders of magnitude more accurate!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.compare_policies(analytic_consumption_policy_2, forward_shooting_policy, kgrid, metric='L2')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, reverse shooting is roughly 100 times more accurate than linearization and slightly more accurate that forward shooting." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "model.compare_policies(analytic_consumption_policy_2, reverse_shooting_policy, kgrid, metric='L2')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These results generalize: \n", "\n", "1. Unless the interval you are interested in is *very* small, then shooting methods will be signficantly more accurate (if somewhat slower) than linearization. \n", "2. Reverse shooting is typically faster and more accurate than forward shooting and thus is the generally preferable approach." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Task 5: Impulse response functions\n", "\n", "Suppose that we start the optimal growth model economy in steady-state with current parameter values and then at $t=2013$ shock the economy by doubling the discount rate. How should we expect the time-paths of capital, output, and consumption per effective worker to respond? How about the time paths of the return to capital and the real wage?\n", "\n", "We can assess the short-run comparative statics of the optimal growth model by plotting [impulse response functions (IRFs)](http://en.wikipedia.org/wiki/Impulse_response). The `plot_impulse_response` method of the `ramsey.Model` class can be used to plot three different types of impulse response functions (i.e., ``efficiency_units``, ``per_capita``, or ``levels``) following a multiplicative shock to any of the models exogenous parameters. \n", "\n", "The syntax for plotting the IRFs (computed using forward shooting) for capital per effective worker, consumption per effective worker, and the real return to capital in response to a doubling of the discount rate would look as follows.\n", "\n", " model.plot_impulse_response(variables=['k', 'c', 'r'], method='forward_shooting', kind='efficiency_units',\n", " param='rho', shock=1.0, T=50, log=False, figsize=(10,6))\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# check the docstring!\n", "model.plot_impulse_response?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# calibrate the model to the UK \n", "ramsey.calibrate_ces(model, iso3_code='GBR', bounds=[1950, 2011], \n", " sigma0=0.5, alpha0=0.5, theta0=2.5, rho=0.04)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# can also plot irfs for all variables at once!\n", "model.plot_impulse_response(variables='all', method='forward_shooting', kind='efficiency_units',\n", " param='rho', shock=2.0, T=50, figsize=(8,12), log=False)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the time of the change in $\\rho$, the value of $k$ (i.e., the stock of capital per effective worker) is predetermined and thus cannot change discontinuously. By contrast consumption per effective worker, $c$, is free to jump at the time of the shock. In the instant following the change in $\\rho$ the representative household increases its consumption. Why? The increase in $\\rho$ means that they are now discounting future utility from consumption more heavily which means that the value of current consumption has increased (relative to future consumption)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5\n", "\n", "This exercise asks you to consider the impact of a productivity growth slowdown on savings. Consider a Ramsey-Cass-Koopmans economy that is on its balanced growth path, and suppose that there is a permanent fall in the rate of technological growth, $g$.\n", "\n", "#### Part a)\n", "Making use of the `plot_phase_diagram` method show how, if at all, does the fall in $g$ affect the $\\dot{k}=0$ and $\\dot{c}=0$ curves?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# INSERT YOUR CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution\n", "%run -i exercise_5a.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Part b)\n", "\n", "What happens to consumption per effective worker, $c$, at the time of the change?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# INSERT YOUR CODE HERE!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# my solution \n", "%run -i exercise_5b.py" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "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.4.4" } }, "nbformat": 4, "nbformat_minor": 0 }