I highly recommend that you take a glance through Frank Ramsey's original paper, entitled A Mathematical Theory of Saving. In particular the introduction where Ramsey states that the common practice amongst economists and decision theorists of discounting the future was "ethicially indefensible and arises mearly from the weakness of the imagination." Ramsey leaves little doubt that his paper was normative and not positive in character. Finally, there is not a single citation at the end of Ramsey's paper: he simply signed his name! Apparently, for this model at least, Ramsey owed no intellectual debts.
%load_ext autoreload
%autoreload 2
# These libaries will be doing the bulk of the heavy lifting!
import numpy as np
import pandas as pd
import sympy as sp
from scipy import linalg
# these models contain all of the code
import ramsey
import pwt
# matplotlib is for graphics
import matplotlib as mpl
import matplotlib.pyplot as plt
from sympy.interactive import printing
printing.init_printing(use_latex=True)
# define symbolic variables
k = sp.symbols('k')
# define symbolic parameters
alpha, sigma = sp.symbols('alpha,sigma')
# define symbolic expressions
rho = (sigma - 1) / sigma
f = (alpha * k**rho + (1 - alpha))**(1 / rho)
The code in the cell below will compute the derivative of $f$ with respect to $\sigma$ and simplify it as much as possible. You job is to simplify it further! HINT: You may wish to make use of the following relationship:
$$ \alpha_k(k) = \frac{kf'(k)}{f(k)} = \frac{\alpha k^{\sigma-1}{\sigma}}{\alpha k^{\sigma-1}{\sigma} + (1 + \alpha)} \tag{1}$$# always check the docstring!
sp.diff?
# differentiate f with respect to sigma and simplify
df_dsigma = sp.diff(f, sigma).simplify()
# can you simplify this result further?
df_dsigma
The code in the cell below computes...
$$\lim_{\sigma \rightarrow 1}\ \frac{df}{d\sigma}$$What is the sign of the resulting limit? Is is positive? Negative? Or indeterminate? What does this imply about the sign of a small shock to $\sigma$ when the production function is Cobb-Douglas?
# always read the docstring
sp.limit?
# limit as sigma -> 1
sp.limit(e=df_dsigma, z=sigma, z0=1)
Determining the sign of $\frac{df}{d \sigma}$
# create a numpy function from the symbolic expression for df_dsigma
df_dsigma_num = sp.lambdify(args=(k, alpha, sigma), expr=df_dsigma, modules='numpy')
plt.figure(figsize=(8,6))
kgrid = np.linspace(0, 10, 1000)
# specify some numeric values for alpha and sigma (you may wish to change these!)
tmp_alpha = 0.5
tmp_sigma = 0.5
plt.plot(kgrid, df_dsigma_num(kgrid, tmp_alpha, tmp_sigma),
label=r'$\alpha=%g, \sigma=%g$' %(tmp_alpha, tmp_sigma))
plt.axhline(0, ls='dashed', color='k')
# add axes label, title, legend, etc
plt.xlabel('Capital per effective worker, $k$', fontsize=15, family='serif')
plt.ylabel(r"$\frac{df}{d\sigma}$", fontsize=25, family='serif', rotation='horizontal')
plt.title('The effect of a change in $\sigma$ on $f(k)$', fontsize=20, family='serif')
plt.legend(loc=0, frameon=False)
plt.show()
kgrid = np.linspace(0, 10, 1000)
plt.plot(kgrid, sp.lambdify((k, alpha, sigma), alpha_k * sp.log(k) - sp.log(f), modules='numpy')(kgrid, 0.5, 1.5))
plt.axhline(0, ls='dashed', color='k')
plt.ylim(-0.5, 0.1)
plt.show()
plt.plot(kgrid, sp.lambdify((k, alpha, sigma), f, modules='numpy')(kgrid, 0.5, 0.5))
plt.plot(kgrid, kgrid, 'k--')
plt.show()
# compute marginal product
mpk = f.diff(k)
mpk
sp.limit(mpk, sigma, 1)
mpk.diff(sigma)
dmpk_dsigma = mpk.diff(sigma)
dmpk_dsigma_num = sp.lambdify((k, alpha, sigma), dmpk_dsigma, modules='numpy')
plt.figure(figsize=(8,6))
kgrid = np.linspace(0, 10, 1000)
plt.plot(kgrid, dmpk_dsigma_num(kgrid, 0.9, 0.9))
plt.axhline(0, ls='dashed', color='k')
#plt.ylim(-0.5, 0.1)
plt.xlabel('Capital per effective worker, $k$', fontsize=15, family='serif')
plt.ylabel(r"$\frac{df'}{d\sigma}$", fontsize=25, family='serif', rotation='horizontal')
plt.title('The effect of a change in $\sigma$ on MPK with $\sigma<1$', fontsize=20, family='serif')
plt.show()
cs = plt.contour(X, Y, numeric_tmp(X, Y, 15, 1.0, 0.5, 0.5, 0.5), 30)
plt.clabel(cs)
plt.xlabel('x', fontsize=15, family='serif')
plt.ylabel('y', fontsize=15, family='serif', rotation='horizontal')
plt.show()
Copy the code for the following Python functions from the lab-2.ipynb
:
ces_output
marginal_product_capital
equation_motion_capital
consumption_euler_equation
crra_utility
ramsey_jacobian
analytic_k_star
analytic_c_star
Define a Python dictionary called baseline_params
containing the model parameters with the following values:
Create an instance of the ramsey.Model
class called model
using the functions and parameters defined above and then add the steady state functions to the model
object using the code from the lab. Finally compute the steady state values. If you have done everything correctly, then you should get the following...
{c_star: 3.07831598008, k_star: 18.1077410593}
def ces_output(t, k, params):
"""
Constant elasticity of substitution (CES) production function.
Arguments:
t: (array) Time.
k: (array) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.
Returns:
y: (array-like) Output (per person/ effective person)
"""
# extract params
alpha = params['alpha']
sigma = params['sigma']
gamma = (sigma - 1) / sigma
# nest Cobb-Douglas as special case
if gamma == 0:
y = k**alpha
else:
y = (alpha * k**gamma + (1 - alpha))**(1 / gamma)
return y
def marginal_product_capital(t, k, params):
"""
Marginal product of capital for CES production function.
Arguments:
t: (array-like) Time.
k: (array-like) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.
Returns:
mpk: (array-like) Derivative of output with respect to capital, k.
"""
# extract params
alpha = params['alpha']
sigma = params['sigma']
gamma = (sigma - 1) / sigma
# nest Cobb-Douglas as special case
if gamma == 0:
mpk = alpha * k**(alpha - 1)
else:
mpk = alpha * k**(gamma - 1) * (alpha * k**gamma + (1 - alpha))**((1 / gamma) - 1)
return mpk
def equation_motion_capital(t, vec, params):
"""
Equation of motion for capital (per person/ effective person).
Arguments:
t: (array-like) Time.
vec: (array-like) Vector of endogenous variables [k, c] where k is
capital (per person/effective person) and c is consumption
(per person/effective person).
params: (dict) Dictionary of parameter values.
Returns:
k_dot: (array-like) Rate of change in the stock of captial
(per person/effective person).
"""
# extract params
n = params['n']
g = params['g']
delta = params['delta']
# unpack the vector of endogenous variables
k, c = vec
# formula for the equation of motion for capital
k_dot = ces_output(t, k, params) - (n + g + delta) * k - c
return k_dot
def consumption_euler_equation(t, vec, params):
"""
Equation of motion for consumption (per person/ effective person).
Arguments:
t: (array-like) Time.
vec: (array-like) Vector of endogenous variables [k, c] where k is
capital (per person/effective person) and c is consumption
(per person/effective person).
params: (dict) Dictionary of parameter values.
Returns:
c_dot: (array-like) Rate of change in consumption (per person/effective
person).
"""
# extract params
g = params['g']
delta = params['delta']
rho = params['rho']
theta = params['theta']
# unpack the vector of endogenous variables
k, c = vec
# marginal product of capital
mpk = marginal_product_capital(t, k, params)
# formula for the equation of motion for consumption
c_dot = ((mpk - delta - rho - theta * g) / theta) * c
return c_dot
def crra_utility(c, params):
"""
The optimal growth model assumes that the representative individual has CRRA utility.
utility. Individual utility is a function of consumption per effective person, c. The
parameter, 0 < theta, plays a dual role:
1) theta is the agent's coefficient of relative risk aversion
2) 1 / theta is the agent's elasticity of inter-temporal substitution
N.B.: If theta = 1, then CRRA utility is equivalent to log utility!
Arguments:
c: (array) Consumption per effective person
params: (dict) Dictionary of model parameters.
Returns:
utility: (array) An array of values respresenting utility from consumption per
effective worker.
"""
# extract the params
theta = params['theta']
if theta == 1:
utility = np.log(c)
else:
utility = (c**(1 - theta) - 1) / (1 - theta)
return utility
def ramsey_jacobian(t, vec, params):
"""
The Jacobian for the optimal growth model is computed using symbolic
differentiation and then evaluated numerically.
Arguments:
t: (array-like) Time.
vec: (array-like) Vector of endogenous variables [k, c] where k is
capital (per person/effective person) and c is consumption
(per person/effective person).
params: (dict) Dictionary of parameter values.
Returns:
jac: (array-like) Jacobian matrix of partial derivatives.
"""
# declare two symbolic variables
k, c = sp.symbols('k,c')
# represent the system of equations as a SymPy matrix
k_dot = lambda k, c: equation_motion_capital(t, [k, c], params)
c_dot = lambda k, c: consumption_euler_equation(t, [k, c], params)
F = sp.Matrix([k_dot(k, c), c_dot(k, c)])
# now computing the Jacobian is trivial!
symbolic_jacobian = F.jacobian([k, c])
# lambdify the function in order to evaluate it numerically
numeric_jacobian = sp.lambdify(args=[k, c], expr=symbolic_jacobian, modules='numpy')
# output num_jac evaluated at vec
evaluated_jacobian = np.array(numeric_jacobian(vec[0], vec[1]))
return evaluated_jacobian
def analytic_k_star(params):
"""Steady-state level of capital (per person/effective person)."""
# extract params
n = params['n']
g = params['g']
alpha = params['alpha']
delta = params['delta']
rho = params['rho']
theta = params['theta']
sigma = params['sigma']
gamma = (sigma - 1) / sigma
# nest Cobb-Douglas as special case
if gamma == 0:
k_star = (alpha / (delta + rho + theta * g))**(1 / (1 - alpha))
else:
k_star = (1 - alpha)**(1 / gamma) * ((alpha / (delta + rho + theta * g))**(gamma / (gamma - 1)) - alpha)**(-1 / gamma)
return k_star
def analytic_c_star(params):
"""Steady-state level of consumption (per person/effective person)."""
# extract params
n = params['n']
g = params['g']
delta = params['delta']
# compute k_star
k_star = analytic_k_star(params)
# compute c_star
c_star = ces_output(0, k_star, params) - (n + g + delta) * k_star
return c_star
# start by defining a dictionary of parameter values
baseline_params = {'alpha':0.5, 'n':0.01, 'g':0.015, 'rho':0.04, 'theta':2.5,
'delta':0.04, 'sigma':1.0, 'A0':1.0, 'L0':1.0}
# we create an instance of the optimgal growth model as follows...
model = ramsey.Model(output=ces_output,
mpk=marginal_product_capital,
k_dot=equation_motion_capital,
c_dot=consumption_euler_equation,
utility=crra_utility,
jacobian=ramsey_jacobian,
params=baseline_params)
# create a dictionary containing the steady state expressions
steady_state_funcs = {'k_star':analytic_k_star, 'c_star':analytic_c_star}
# pass it as an arg to the set_functions method
model.steady_state.set_functions(func_dict=steady_state_funcs)
# compute the steady state values!
model.steady_state.set_values()
# display the values as follows
model.steady_state.values
{c_star: 3.07831598008, k_star: 18.1077410593}
Plot a phase diagram for the model using the plot_phase_diagram
method.
plt.figure(figsize=(8,6))
model.plot_phase_diagram(gridmax=100, arrows=True)
plt.show()
Plot the phase diagram following a 5% decrease in the value of $\sigma$. Clearly explain why each of the $\dot{c}=0$ and the $\dot{k}=0$ shifts.
plt.figure(figsize=(8,6))
model.plot_phase_diagram(gridmax=100, param='sigma', arrows=True, shock=0.95, reset=True)
plt.show()
model.solve_linearization(
plt.figure(figsize=(8,6))
k_star = model.steady_state.values['k_star']
kgrid = np.linspace(0, 2.0 * k_star, 1000)
# specify some numeric values for alpha and sigma
tmp_alpha = model.params['alpha']
tmp_sigma = model.params['sigma']
plt.plot(kgrid, df_dsigma_num(kgrid, tmp_alpha, tmp_sigma),
label=r'$\alpha=%g, \sigma=%g$' %(tmp_alpha, tmp_sigma))
plt.axhline(0, ls='dashed', color='k')
# add axes label, title, legend, etc
plt.xlabel('Capital per effective worker, $k$', fontsize=15, family='serif')
plt.ylabel(r"$\frac{df}{d\sigma}$", fontsize=25, family='serif', rotation='horizontal')
plt.title('The effect of a change in $\sigma$ on $f(k)$', fontsize=20, family='serif')
plt.legend(loc=0, frameon=False)
plt.show()
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) <ipython-input-424-38d65052a851> in <module>() 8 tmp_sigma = model.params['sigma'] 9 ---> 10 plt.plot(kgrid, df_dsigma_num(kgrid, tmp_alpha, tmp_sigma), 11 label=r'$\alpha=%g, \sigma=%g$' %(tmp_alpha, tmp_sigma)) 12 /Users/clarissasweet/Library/Enthought/Canopy_32bit/User/lib/python2.7/site-packages/numpy/__init__.pyc in <lambda>(k, alpha, sigma) ZeroDivisionError: float division by zero
Figure(640x480)
plt.figure(figsize=(8,6))
# specify some numeric values for alpha and sigma (you may wish to change these!)
tmp_alpha = model.params['alpha']
tmp_sigma = model.params['sigma']
plt.plot(kgrid, dmpk_dsigma_num(kgrid, tmp_alpha, tmp_sigma),
label=r'$\alpha=%g, \sigma=%g$' %(tmp_alpha, tmp_sigma))
plt.axhline(0, ls='dashed', color='k')
plt.xlabel('Capital per effective worker, $k$', fontsize=15, family='serif')
plt.ylabel(r"$\frac{df'}{d\sigma}$", fontsize=25, family='serif', rotation='horizontal')
plt.title(r'The effect of a change in $\sigma$ on MPK', fontsize=20, family='serif')
plt.legend(loc=0, frameon=False)
plt.show()
dmpk_dsigma_num(1, tmp_alpha, tmp_sigma)
def savings_rate(k):
c = consumption_policy_func(k)
y = model.output(0, k, model.params)
s = 1 - (c / y)
return s
kgrid = np.linspace(kmin, kmax, 1000)
plt.plot(kgrid, model.get_capitals_share(kgrid))
plt.show()
savings_rate(k_star)
0.35802484338211749
k_star
19.30052917806588
kgrid = np.linspace(kmin, kmax, 1000)
plt.plot(kgrid, savings_rate(kgrid))
plt.show()
kgrid = np.linspace(kmin, kmax, 1000)
plt.plot(kgrid, consumption_policy_func(kgrid))
plt.show()
kgrid = np.linspace(kmin, kmax, 1000)
plt.plot(savings_rate(kgrid), model.get_capitals_share(kgrid))
plt.show()
model.get_capitals_share(kgrid)