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

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

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

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

...is the Pratt-Arrow measure of absolute risk aversion.

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