In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy
import scipy.interpolate as interp
import scipy.optimize as opt
%matplotlib inline
In [2]:
# parameters
rho = 0.05 # discount rate
alpha = 0.3 # capital elasticity
delta = 0.1 # depreciation
theta = 1.0 # crra param
In [3]:
# algorithm
S0 = 4096 # max number of steps in perturbation simulation
dt0 = 0.01 # default time step size
K = 512 # size of capital grid
ptol = 1.0e-4 # initial displacement for perturbation simulation
dt_max = 1.0 # max time step in variable delt schemes
In [4]:
# steady state
k_ss = (alpha/(rho+delta))**(1.0/(1.0-alpha))
i_ss = delta*k_ss
y_ss = k_ss**alpha
c_ss = y_ss - i_ss
In [5]:
def simulate_ramsey(k0, c0, k_lo=0, k_hi=None, gk_targ=0.001, tau=0, dt=None, direc=-1, S=S0):
    # initial values
    k, c = k0, c0
    kpath, cpath = [k], [c]

    for t in range(S-1):
        # growth rates
        gk = k**(alpha-1) - delta - c/k
        gc = (1/theta)*((1-tau)*(alpha*k**(alpha-1)-delta)-rho)

        # adaptive step size
        dt1 = direc*np.minimum(dt_max, gk_targ/np.abs(gk)) if dt is None else dt
        k *= 1 + dt1*gk
        c *= 1 + dt1*gc

        # store
        kpath.append(k)
        cpath.append(c)

        # terminal condition
        if k_lo is not None and k <= k_lo:
            break
        if k_hi is not None and k >= k_hi:
            break

    return np.array(kpath), np.array(cpath)
In [6]:
def forward_ramsey(k0, cf, dt=dt0, S=S0):
    # initial values
    k, c = k0, cf(k0)
    kpath, cpath = [k], [c]

    for t in range(S-1):
        # growth rates
        gk = k**(alpha-1) - delta - c/k

        # adaptive step size
        k *= 1 + dt*gk
        c = cf(k)

        # store
        kpath.append(k)
        cpath.append(c)

    return np.array(kpath), np.array(cpath)

Part (b)

In [7]:
# capital grid around steady state
k_lo, k_hi = 0.1*k_ss, 2.0*k_ss
k_grid = np.linspace(k_lo, k_hi, K)

# generate upper and lower paths
kvec1, cvec1 = simulate_ramsey(k_ss-ptol, c_ss-ptol, k_lo, k_hi)
kvec2, cvec2 = simulate_ramsey(k_ss+ptol, c_ss+ptol, k_lo, k_hi)

# combined path
kvec = np.concatenate([kvec1[::-1], kvec2])
cvec = np.concatenate([cvec1[::-1], cvec2])

# interpolate
cpol = interp.interp1d(kvec, cvec)
c_grid = cpol(k_grid)
In [8]:
# plot results
fig, ax = plt.subplots()
ax.plot(k_grid, c_grid);
ax.scatter([k_ss], [c_ss]);

Part (c)

In [9]:
T = dt0*S0

k0 = 0.9*k_ss

tfwd = np.linspace(0, T, S0)
kfwd, cfwd = forward_ramsey(k0, cpol)

kconv = tfwd[np.argmax(np.abs(kfwd-k_ss) < 0.01*k_ss)]
cconv = tfwd[np.argmax(np.abs(cfwd-c_ss) < 0.01*c_ss)]
print(f'k converge (1%) = {kconv:.1f} years')
print(f'c converge (1%) = {cconv:.1f} years')

fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
ax0.plot(tfwd, kfwd)
ax0.scatter(0, k_ss)
ax0.hlines(k_ss, 0, T, linestyle='--', linewidth=1)
ax0.set_xlabel('Time (years)')
ax0.set_title('Capital (k)')
ax1.plot(tfwd, cfwd)
ax1.scatter(0, c_ss)
ax1.hlines(c_ss, 0, T, linestyle='--', linewidth=1)
ax1.set_xlabel('Time (years)')
ax1.set_title('Consumption (c)');
k converge (1%) = 12.7 years
c converge (1%) = 9.7 years

Part (d)

In [10]:
tau = 0.2
In [11]:
# steady state
k_ss1 = (alpha/(delta+rho/(1-tau)))**(1.0/(1.0-alpha))
i_ss1 = delta*k_ss1
y_ss1 = k_ss1**alpha
c_ss1 = y_ss1 - i_ss1
In [12]:
# capital grid around steady state
k_lo, k_hi = 0.1*k_ss1, 2.0*k_ss1

# generate upper and lower paths
kvec1, cvec1 = simulate_ramsey(k_ss1-ptol, c_ss1-ptol, k_lo, k_hi, tau=tau)
kvec2, cvec2 = simulate_ramsey(k_ss1+ptol, c_ss1+ptol, k_lo, k_hi, tau=tau)

# combined path
kvec = np.concatenate([kvec1[::-1], kvec2])
cvec = np.concatenate([cvec1[::-1], cvec2])

# interpolate
cpol1 = interp.interp1d(kvec, cvec)
In [13]:
T = dt0*S0

k0 = k_ss

tfwd = np.linspace(0, T, S0)
kfwd, cfwd = forward_ramsey(k0, cpol1)

fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
ax0.plot(tfwd, kfwd)
ax0.scatter(0, k_ss)
ax0.hlines(k_ss1, 0, T, linestyle='--', linewidth=1)
ax0.set_xlabel('Time (years)')
ax0.set_title('Capital (k)')
ax1.plot(tfwd, cfwd)
ax1.scatter(0, c_ss)
ax1.hlines(c_ss1, 0, T, linestyle='--', linewidth=1)
ax1.set_xlabel('Time (years)')
ax1.set_title('Consumption (c)');
In [14]:
T = 5
dt = 0.01
S = int(T/dt)
S1 = S + S0
T1 = dt*S1

def cons_err(c0):
    kfwd, cfwd = simulate_ramsey(k_ss, c0, dt=dt, S=S, direc=1)
    k1, c1 = kfwd[-1], cfwd[-1]
    cp = cpol1(kfwd[-1])
    return c1 - cp

c_opt = opt.bisect(cons_err, 0.8*c_ss, 1.2*c_ss)

kfwd1, cfwd1 = simulate_ramsey(k_ss, c_opt, dt=dt, S=S, direc=1)
k1, c1 = kfwd1[-1], cfwd1[-1]
kfwd2, cfwd2 = forward_ramsey(k1, cpol1)

tfwd = np.linspace(0, T1, S1)
kfwd = np.concatenate([kfwd1, kfwd2])
cfwd = np.concatenate([cfwd1, cfwd2])

fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
ax0.plot(tfwd, kfwd)
ax0.scatter([0, 5], [k_ss, k1])
ax0.hlines(k_ss1, 0, T1, linestyle='--', linewidth=1)
ax0.set_xlabel('Time (years)')
ax0.set_title('Capital (k)')
ax1.plot(tfwd, cfwd)
ax1.scatter([0, 5], [c_ss, c1])
ax1.hlines(c_ss1, 0, T1, linestyle='--', linewidth=1)
ax1.set_xlabel('Time (years)')
ax1.set_title('Consumption (c)');
In [ ]: