In [1]:
%load_ext autoreload
In [2]:
%autoreload 2
In [3]:
%matplotlib inline
In [4]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn
import pycollocation
Suppose that there exist a representative household with the following lifetime utility...
$$U \equiv \int_{t=0}^{\infty} e^{-\rho t}u(c(t))N(t)dt$$...where the flow utility function $u(C(t))$ is assumed to be a concave function of per capita consumption, $c(t)$. Note that $N(t)$, the size of the household, is assumed to grow at a constant and exogenous rate $n$:
$$ N(t) = N(0)e^{nt}.$$At each instant in time, the representative household faces the following budget constraint...
$$\dot{K}(t) = r(t)K(t) + W(t)lN(t) - P(t)c(t)N(t)$$...where $r(t)$ is the real interest rate, $K(t)$ is the total amount of household capital; $W(t)$ is the real wage paid for labor, $l$ is the household's labor endowment; $P(t)$ is the price level of the consumption good (which we are free to normalize to 1 for convenience).
Form the Hamiltonian...
$$ H(t, K, c, \lambda) \equiv e^{-\rho t}u(c(t))N(t) + \lambda(t)\bigg[r(t)K(t) + W(t)lN(t) - c(t)N(t)\bigg] $$...differentiate with respect to control variables $c$ and the state variable $K$...
\begin{align} \frac{\partial H}{\partial c} \equiv& e^{-\rho t}\frac{\partial u}{\partial c}N(t) - \lambda(t)N(t) \\ \frac{\partial H}{\partial K} \equiv& r(t)\lambda(t) \end{align}...the state and co-state equations are...
\begin{align} \dot{K}(t) = \frac{\partial H}{\partial \lambda} =& r(t)K(t) + W(t)lN(t) - c(t)N(t) \\ \dot{\lambda} = -\frac{\partial H}{\partial K} =& -r(t)\lambda(t)\\ \end{align}...differentiate with respect to time...
$$ -\rho e^{-\rho t}\frac{\partial u}{\partial c} + e^{-\rho t}\frac{\partial^2 u}{\partial c^2}\dot{c}(t)=\dot{\lambda}(t)$$...using the co-state equation we can write...
$$ -\rho e^{-\rho t}\frac{\partial u}{\partial c} + e^{-\rho t}\frac{\partial^2 u}{\partial c^2}\dot{c}(t)=-r(t)\lambda(t)$$...using the first-order condition again we can write...
$$ -\rho e^{-\rho t}\frac{\partial u}{\partial c} + e^{-\rho t}\frac{\partial^2 u}{\partial c^2}\dot{c}(t)=-r(t)e^{-\rho t}\frac{\partial u}{\partial c}$$...finally, after a bit of algebra we find that the consumption behavior of the representative household is described by the following Euler equation...
$$ \frac{\dot{c}(t)}{c(t)} = \frac{1}{RRA(c(t))}\bigg[r(t) - \rho\bigg] $$...where...
$$ RRA(c(t)) = -\frac{c\frac{\partial^2 u}{\partial c^2}}{\frac{\partial u}{\partial c}}$$...is the Pratt-Arrow measure of relative risk aversion. Consumption Euler equation says that consumption growth is proportional to the gap between the real interest rate $r(t)$ and the discount rate $\rho$; and inversely proportional to risk preferences.
Inputs to production are capital, $K$, labor, $L$, and technology, $A$. Representative firm has a production function, $F$, that is homogenous of degree 1...
\begin{equation} Y(t) = F\bigg(A(t), K(t), L(t)\bigg) \end{equation}...note that technology, $A$, is assumed to grow at a constant, exogenous rate $g$
$$ A(t) = A(0)e^{gt}. $$Firms choose demands for $K$ and $L$ in order to maximise profits...
\begin{align} \Pi(t) =& P(t)Y(t) - (r(t) + \delta)K(t) - W(t)L(t) \\ =& F\bigg(A(t), K(t), L(t)\bigg) - (r(t) + \delta)K(t) - W(t)L(t) \end{align}Standard assumptions of perfect competition in factor markets as well as the production of output goods implies that inputs to production are paid their marginal products. Thus the real interest rate $r(t)$ and the wage, $W(t)$ are...
\begin{align} r(t) =& \frac{\partial F}{\partial K} - \delta \\ W(t) =& \frac{\partial F}{\partial L}. \end{align}Note that our homogeneity assumption on $F$ combined with Euler's other theorem are enough to insure that the representative firm earns zero profit...
$$ \Pi(t) = F\bigg(A(t), K(t), L(t)\bigg) - \bigg[K\frac{\partial F}{\partial K} + L\frac{\partial F}{\partial L}\bigg] = 0. $$Impose market clearing equilibrium assumption for capital and labor markets (market for consumption goods clears by Walras Law). In particular for the labor market...
$$ \text{labor demand} \equiv L(t) = lN(t) \equiv \text{labor supply}. $$Note that the explicit assumption of market clearing equilibrium really imposes two implicit assumptions:
Assumption that factor markets clear allows us to combine household and firm behavior to get a system of ordinary differential equations...
\begin{align} \dot{c}(t) =& \frac{1}{ARA(c(t))}\bigg[\frac{\partial F}{\partial K} - \delta - \rho\bigg] \\ \dot{K}(t) =& F\bigg(A(t), K(t), lN(t)\bigg) - \delta K(t) - c(t)N(t) \\ \end{align}...where...
$$ ARA(c(t)) = -\frac{\frac{\partial^2 u}{\partial c^2}}{\frac{\partial u}{\partial c}}$$Recall that technology and household size (i.e., population!) are growing at exogenous rates. Want to de-trend the model so that we can analyze fixed point equilibrium of the dynamical system.
Homogeneity assumption on $F$ tells us that...
$$\frac{1}{A(t)N(t)}F\bigg(A(t), K(t), lN(t)\bigg) = F\bigg(1, \frac{K(t)}{A(t)N(t)}, l\bigg) = f(k(t)) $$...and...
$$ \frac{\partial F}{\partial K} = \frac{\partial F}{\partial K}\bigg(1, \frac{K(t)}{A(t)N(t)}, l\bigg) = \frac{\partial f}{\partial k} = f'(k(t)) $$...where...
$$ \tilde{k}(t) = \frac{K(t)}{A(t)N(t)} $$...is now defined to be capital per unit effective labor supply. Similarly we can de-trend per capita consumption $c(t)$ using the following change of variables...
\begin{align} c(t) =& A(t)\tilde{c}(t) \\ \dot{c} =& \dot{A}(t)\tilde{c}(t) + A(t)\dot{\tilde{c}}(t) \end{align}...where...
$$ \tilde{c}(t) \equiv \frac{c(t)}{A(t)}$$is consumption per unit of effective labor supply.
Using these results, we can re-write the above system of differential equations as follows...
\begin{align} \dot{\tilde{k}}(t) =& f(\tilde{k}(t)) - (g + n + \delta)\tilde{k}(t) - \tilde{c}(t)\\ \dot{\tilde{c}}(t) =& \frac{1}{A(0)e^{gt}ARA\bigg(A(0)e^{gt}\tilde{c}(t)\bigg)}\bigg[f'(k(t)) - \delta - \rho\bigg] - g\tilde{c}(t) \end{align}...where $k$ and $c$ are now measured per unit of effective labor supply. Note that the equation of motion for consumption per unit effective labor supply is no longer time invariant. However, if marginal utility of consumption for the representative household is homogenous of degree $k$, then the equation of motion for $\tilde{c}$ becomes time-invariant...
\begin{align} \dot{\tilde{c}}(t) =& \frac{1}{ARA\bigg(\tilde{c}(t)\bigg)}\bigg[f'(k(t)) - \delta - \rho\bigg] - g\tilde{c}(t) \end{align}To complete the model we need to specifcy some boundary conditions. The initial conditions for technology, $A(0)$, and household size (i.e., population), $N(0)$, and capital, $K(0)$, are assumed given. Therefore...
$$ \tilde{k}(0) = \frac{K(0)}{A(0)N(0)}. $$We will impose the following terminal condition on consumption per unit effective labor supply...
$$ \lim_{t \rightarrow \infty} \tilde{c}(t) = \tilde{c}^*. $$The full model is completely specified by the following system of ordinary conditions and boundary conditions...
\begin{align} \dot{\tilde{k}}(t) =& f(\tilde{k}(t)) - (g + n + \delta)\tilde{k}(t) - \tilde{c}(t),\ \tilde{k}(0) = \tilde{k}_0 \\ \dot{\tilde{c}}(t) =& \frac{1}{A(0)e^{gt}ARA\bigg(A(0)e^{gt}\tilde{c}(t)\bigg)}\bigg[f'(\tilde{k}(t)) - \delta - \rho\bigg] - g\tilde{c}(t),\ \lim_{t \rightarrow \infty} \tilde{c}(t) = c^* \end{align}Suppose representative household preferences are consistent with Hyperbolic Relative Risk Aversion (HARA)...
$$ ARA(c) = \frac{1}{ac + b} $$...in this case the equation of motion for consumption per unit effective labor reduces to...
$$\dot{\tilde{c}}(t) = \bigg[a\tilde{c}(t) + \frac{b}{A(0)e^{gt}}\bigg]\bigg[f'(\tilde{k}(t)) - \delta - \rho\bigg] - g\tilde{c}(t) $$...note that when $b \ne 0$ this is a non-autonomous differential equation and the representative household has time-varying risk preferences. On the other hand, if $b = 0$ then the representative household has standard CRRA risk preferences (which do not depend directly on time).
In [5]:
def hara(t, c, a, b, **params):
"""
Hyperbolic Absolute Risk Aversion (HARA).
Notes
-----
For Constant Absolute Risk Aversion (CARA), set a=0; for
Constant Relative Risk Aversion (CRRA), set b=0.
"""
return 1 / (a * c + b)
def cobb_douglas_output(k_tilde, alpha, l, **params):
return k_tilde**alpha * l**(1 - alpha)
def cobb_douglas_mpk(k_tilde, alpha, l, **params):
return alpha * k_tilde**(alpha - 1) * l**(1 - alpha)
def c_tilde_dot(t, k_tilde, c_tilde, A0, delta, g, rho, **params):
A = A0 * np.exp(g * t)
r = cobb_douglas_mpk(k_tilde, **params) - delta
return ((r - rho) / (A * hara(t, A * c_tilde, **params))) - g * c_tilde
def k_tilde_dot(t, k_tilde, c_tilde, delta, g, n, **params):
return cobb_douglas_output(k_tilde, **params) - c_tilde - (g + n + delta) * k_tilde
def standard_ramsey_model(t, k_tilde, c_tilde, A0, delta, g, n, rho, **params):
out = [k_tilde_dot(t, k_tilde, c_tilde, delta, g, n, **params),
c_tilde_dot(t, k_tilde, c_tilde, A0, delta, g, rho, **params)]
return out
def initial_condition(t, k_tilde, c_tilde, A0, K0, N0, **params):
return [k_tilde - (K0 / (A0 * N0))]
def terminal_condition(t, k_tilde, c_tilde, **params):
return [c_tilde - equilibrium_consumption(**params)]
def equilibrium_capital(a, alpha, b, delta, g, l, n, rho, **params):
return ((a * alpha * l**(1 - alpha)) / (a * (delta + rho) + g))**(1 / (1 - alpha))
def equilibrium_consumption(a, alpha, b, delta, g, l, n, rho, **params):
kss = equilibrium_capital(a, alpha, b, delta, g, l, n, rho)
return cobb_douglas_output(kss, alpha, l) - (g + n + delta) * kss
To complete the model we need to define some parameter values.
In [6]:
# set b=0 for CRRA...
params = {'a': 1.0, 'b': 0.0, 'g': 0.02, 'n': 0.02, 'alpha': 0.15,
'delta': 0.04, 'l': 1.0, 'K0': 1.0, 'A0': 1.0, 'N0': 1.0,
'rho': 0.02}
In [14]:
pycollocation.problems.TwoPointBVP?
In [7]:
standard_ramsey_bvp = pycollocation.problems.TwoPointBVP(bcs_lower=initial_condition,
bcs_upper=terminal_condition,
number_bcs_lower=1,
number_odes=2,
params=params,
rhs=standard_ramsey_model,
)
In [8]:
def initial_mesh(t, T, num, problem):
# compute equilibrium values
cstar = equilibrium_consumption(**problem.params)
kstar = equilibrium_capital(**problem.params)
ystar = cobb_douglas_output(kstar, **problem.params)
# create the mesh for capital
ts = np.linspace(t, T, num)
k0 = problem.params['K0'] / (problem.params['A0'] * problem.params['N0'])
ks = kstar - (kstar - k0) * np.exp(-ts)
# create the mesh for consumption
s = 1 - (cstar / ystar)
y0 = cobb_douglas_output(k0, **problem.params)
c0 = (1 - s) * y0
cs = cstar - (cstar - c0) * np.exp(-ts)
return ts, ks, cs
In [274]:
pycollocation.solvers.Solver?
In [9]:
polynomial_basis = pycollocation.basis_functions.PolynomialBasis()
solver = pycollocation.solvers.Solver(polynomial_basis)
boundary_points = (0, 200)
ts, ks, cs = initial_mesh(*boundary_points, num=1000, problem=standard_ramsey_bvp)
basis_kwargs = {'kind': 'Chebyshev', 'domain': boundary_points, 'degree': 25}
k_poly = polynomial_basis.fit(ts, ks, **basis_kwargs)
c_poly = polynomial_basis.fit(ts, cs, **basis_kwargs)
initial_coefs = np.hstack([k_poly.coef, c_poly.coef])
nodes = polynomial_basis.roots(**basis_kwargs)
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs,
nodes, standard_ramsey_bvp)
In [10]:
ts, _, _ = initial_mesh(*boundary_points, 1000, standard_ramsey_bvp)
k_soln, c_soln = solution.evaluate_solution(ts)
plt.plot(ts, k_soln)
plt.plot(ts, c_soln)
plt.show()
In [11]:
k_resids, c_resids = solution.evaluate_residual(ts)
plt.plot(ts, k_resids)
plt.plot(ts, c_resids)
plt.show()
In [12]:
k_normalized_resids, c_normalized_resids = solution.normalize_residuals(ts)
plt.plot(ts, np.abs(k_normalized_resids))
plt.plot(ts, np.abs(c_normalized_resids))
plt.yscale('log')
plt.show()
In [13]:
bspline_basis = pycollocation.basis_functions.BSplineBasis()
solver = pycollocation.solvers.Solver(bspline_basis)
boundary_points = (0, 200)
ts, ks, cs = initial_mesh(*boundary_points, num=250, problem=standard_ramsey_bvp)
tck, u = bspline_basis.fit([ks, cs], u=ts, k=5, s=0)
knots, coefs, k = tck
initial_coefs = np.hstack(coefs)
basis_kwargs = {'knots': knots, 'degree': k, 'ext': 2}
nodes = np.linspace(*boundary_points, num=249)
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs,
nodes, standard_ramsey_bvp)
In [14]:
ts, _, _ = initial_mesh(*boundary_points, num=1000, problem=standard_ramsey_bvp)
k_soln, c_soln = solution.evaluate_solution(ts)
plt.plot(ts, k_soln)
plt.plot(ts, c_soln)
plt.show()
In [15]:
k_resids, c_resids = solution.evaluate_residual(ts)
plt.plot(ts, k_resids)
plt.plot(ts, c_resids)
plt.show()
In [16]:
k_normalized_resids, c_normalized_resids = solution.normalize_residuals(ts)
plt.plot(ts, np.abs(k_normalized_resids))
plt.plot(ts, np.abs(c_normalized_resids))
plt.yscale('log')
plt.show()
Can we refactor the above code so that we can solve a Ramsey-Cass-Koopmans model for arbitrary intensive production functions and risk preferences? Yes!
In [17]:
from pycollocation.tests import models
Example usage...
In [18]:
def ces_output(k_tilde, alpha, l, sigma, **params):
gamma = (sigma - 1) / sigma
if gamma == 0:
y = k_tilde**alpha * l**(1 - alpha)
else:
y = (alpha * k_tilde**gamma + (1 - alpha) * l**gamma)**(1 / gamma)
return y
def ces_mpk(k_tilde, alpha, l, sigma, **params):
y = ces_output(k_tilde, alpha, l, sigma)
gamma = (sigma - 1) / sigma
if gamma == 0:
mpk = alpha * (y / k_tilde)
else:
mpk = alpha * k_tilde**(gamma - 1) * (y / (alpha * k_tilde**gamma + (1 - alpha) * l**gamma))
return mpk
def ces_equilibrium_capital(a, alpha, b, delta, g, l, n, rho, sigma, **params):
"""Steady state value for capital stock (per unit effective labor)."""
gamma = (sigma - 1) / sigma
if gamma == 1:
kss = ((a * alpha * l**(1 - alpha)) / (a * (delta + rho) + g))**(1 / (1 - alpha))
else:
kss = l * ((1 / (1 - alpha)) * (((a * (delta + rho) + g) / (a * alpha))**(gamma / (1 - gamma)) - alpha))**(-1 / gamma)
return kss
In [19]:
ces_params = {'a': 0.5, 'b': 1.0, 'g': 0.02, 'n': 0.02, 'alpha': 0.15,
'delta': 0.04, 'l': 5.0, 'K0': 2.0, 'A0': 1.0, 'N0': 1.0,
'rho': 0.02, 'sigma': 2.0}
In [20]:
generic_ramsey_bvp = models.RamseyCassKoopmansModel(hara,
ces_output,
ces_equilibrium_capital,
ces_mpk,
ces_params)
In [37]:
polynomial_basis = pycollocation.basis_functions.PolynomialBasis()
solver = pycollocation.solvers.Solver(polynomial_basis)
boundary_points = (0, 500)
ts, ks, cs = initial_mesh(*boundary_points, num=1000, problem=generic_ramsey_bvp)
basis_kwargs = {'kind': 'Chebyshev', 'domain': boundary_points, 'degree': 30}
k_poly = polynomial_basis.fit(ts, ks, **basis_kwargs)
c_poly = polynomial_basis.fit(ts, cs, **basis_kwargs)
initial_coefs = np.hstack([k_poly.coef, c_poly.coef])
nodes = polynomial_basis.roots(**basis_kwargs)
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs,
nodes, generic_ramsey_bvp)
In [38]:
solution.result.success
Out[38]:
In [39]:
k_soln, c_soln = solution.evaluate_solution(ts)
plt.plot(ts, k_soln)
plt.plot(ts, c_soln)
plt.show()
In [40]:
k_normalized_resids, c_normalized_resids = solution.normalize_residuals(ts)
plt.plot(ts, np.abs(k_normalized_resids))
plt.plot(ts, np.abs(c_normalized_resids))
plt.yscale('log')
plt.show()
In [45]:
css = generic_ramsey_bvp.equilibrium_consumption(**generic_ramsey_bvp.params)
kss = ces_equilibrium_capital(**generic_ramsey_bvp.params)
plt.plot(k_soln / kss, c_soln / css)
plt.xlabel(r'$\frac{\tilde{k}}{\tilde{k}^*}$', fontsize=20)
plt.ylabel(r'$\frac{\tilde{c}}{\tilde{c}^*}$', fontsize=20, rotation='horizontal')
plt.title("Phase space for the Ramsey-Cass-Koopmans model")
Out[45]:
In [42]:
from mpl_toolkits.mplot3d import Axes3D
In [44]:
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(k_soln / kss, c_soln / css, ts, label='Ramsey model trajectory')
plt.xlabel(r'$\frac{\tilde{k}}{\tilde{k}^*}$', fontsize=20)
plt.ylabel(r'$\frac{\tilde{c}}{\tilde{c}^*}$', fontsize=20)
ax.set_zlabel('$t$')
ax.legend()
plt.show()
In [ ]: