In [1]:




In [2]:




In [3]:

%matplotlib inline




In [4]:

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn

import pycollocation



# Textbook example: Ramsey-Cass-Koopmans model

## Household behavior

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).

### Solution to the household problem

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}

#### Derivation of the Euler equation

From first-order condition for $c$ we have that...

$$e^{-\rho t}\frac{\partial u}{\partial c}=\lambda(t)$$

...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.

## Firm behavior

Inputs to production are capital, $K$, labor, $L$, and technology, $A$. Representative firm has a production function, $F$, that is homogenous of degree 1...

$$Y(t) = F\bigg(A(t), K(t), L(t)\bigg)$$

...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.$$

## Market clearing equilibrium

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:

1. Dynamic adjustment processes necessary to actually clear markets occur at time scales that are very small relevant to the time scale of the model.
2. ??

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}}$$

## De-trending the model

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}

## Boundary conditions

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}^*.$$

## Full model

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}

## Example: HARA risk aversion and Cobb-Douglas production

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}



## Solving the model with pyCollocation

### Defining a pycollocation.TwoPointBVP instance



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,
)



### Finding a good initial guess for $\tilde{k}(t)$

Theory tells us that, starting from some initial condition $\tilde{k}_0$, the solution to the Solow model converges monotonically toward its long run equilibrium value $\tilde{k}^*$. Our initial guess for the solution should preserve this property...



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



### Solving the model



In [274]:

pycollocation.solvers.Solver?



### Polynomial basis functions



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()






### B-spline basis functions



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()






# Generic Ramsey-Cass-Koopmans model

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]:

True




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()






### Phase space plots

#### 2D phase space



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]:

<matplotlib.text.Text at 0x109c263c8>



#### 3D phase space



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()




/Users/drpugh/anaconda/envs/pycollocation-dev/lib/python3.5/site-packages/matplotlib/collections.py:590: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
if self._edgecolors == str('face'):




In [ ]: