{ "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", "