\n",
" \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",
" \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
}