Ramsey model: solution via shooting algorithm¶

In [1]:
%matplotlib inline

import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt

np.seterr(all='ignore')
Out[1]:
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
In [2]:
# Define needed functions

def ProductionFunction(k, params):
    α, ρ, σ, δ, n, g = params
    return k**α

def InterestRate(k, params):
    α, ρ, σ, δ, n, g = params
    
    return α*k**(α-1) - δ
    
def SteadyState(params):
    α, ρ, σ, δ, n, g = params
    
    r_star = (1+n)*(1+ρ)*(1+g)**σ - 1
    k_star = (α/(r_star+δ))**(1/(1-α))
    c_star = ProductionFunction(k_star, params) - (δ+n+g+n*g)*k_star
    
    return k_star, c_star

def ResourceConstraint(k, c, params):
    α, ρ, σ, δ, n, g = params
    
    return (ProductionFunction(k, params) + (1-δ)*k - c)/((1+n)*(1+g))

def EulerEquation(k, c, params):
    α, ρ, σ, δ, n, g = params
    
    k_next = ResourceConstraint(k, c, params)
    
    if k_next > 0:
        r_next = InterestRate(k_next, params)
        c_next = ((1+r_next)/(1+ρ)/(1+n))**(1/σ) * c / (1+g)
        return c_next
    else:
        return 0

def Constant_k(k, params):
    α, ρ, σ, δ, n, g = params
    
    return ProductionFunction(k, params) - (δ+n+g+n*g)*k
In [3]:
# Iterate path forward for given k_0 and c_0

def Path(c_0, k_0, params, T=100):
    
    T += 1
    
    k_t = np.zeros(T)
    c_t = np.zeros(T)
    
    k_t[0] = k_0
    c_t[0] = c_0
    
    for t in range(T-1):
        k_t[t+1] = ResourceConstraint(k_t[t], c_t[t], params)
        if k_t[t+1] > 0:
            c_t[t+1] = EulerEquation(k_t[t], c_t[t], params)
        else:
            k_t[t+1] = 0
            c_t[t+1] = 0
            
    return k_t, c_t

# Convergence criterion

def Path_crit(c_0, k_0, params, T=100):
    
    k_star, c_star = SteadyState(params)
    
    if np.sign( (k_0 - k_star) * (c_0 - Constant_k(k_0, params)) ) != 1:
        return np.inf
    
    else:
        
        k_t, c_t = Path(c_0, k_0, params, T)

        ss_diff = np.sqrt((k_t/k_star-1)**2 + (c_t/c_star-1)**2)
        
        return np.min(ss_diff)

Analytic solution case: $\sigma=1$ and $\delta=1$¶

In [4]:
# Analytic check: σ=1, δ=1 case

# Parameters
α = 0.7
ρ = 1/.95-1
σ = 1
δ = 1
n = 0
g = 0

params_analytic = α, ρ, σ, δ, n, g

# Steady state
print(SteadyState(params_analytic))
k_star, c_star = SteadyState(params_analytic)

(0.2566879516448986, 0.1293089681218662)
In [5]:
# Steady state
k_check = (α/((1+ρ)*(1+n)*(1+g)**σ-(1-δ)))**(1/(1-α))
c_check = Constant_k(k_check, params_analytic)
print(k_check, c_check)

0.2566879516448986 0.1293089681218662
In [6]:
# Convergence for k_0 < k_star

k_0_check = k_star / 1000
c_0_guess = Constant_k(k_0_check, params_analytic)-0.01

res = opt.minimize(Path_crit, c_0_guess, 
                   args=(k_0_check, params_analytic, 100), method='Powell')
print(res)

k_check, c_check = Path(res.x, k_0_check, params_analytic, T=25)

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 1.350134647115135e-07
       x: [ 1.027e-03]
     nit: 6
   direc: [[ 8.132e-16]]
    nfev: 217
In [7]:
# Convergence for k_0 < k_star

plt.plot(k_check, c_check, 'C2.-', lw=2)
plt.plot(k_star, c_star, 'ko')
plt.show()
In [8]:
# Convergence for k_0 > k_star

k_0_check_up = 1.5*k_star
c_0_guess_up = Constant_k(k_0_check_up, params_analytic)+0.01

res = opt.minimize(Path_crit, c_0_guess_up, 
                   args=(k_0_check_up, params_analytic, 100), method='Powell')
print(res)

k_check_up, c_check_up = Path(res.x, k_0_check_up, params_analytic, T=25)

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 3.1401956880202074e-08
       x: [ 1.717e-01]
     nit: 6
   direc: [[-8.493e-18]]
    nfev: 191
In [9]:
plt.plot(k_check_up, c_check_up, 'C2.-', lw=2)
plt.plot(k_star, c_star, 'ko')
plt.show()
In [10]:
# Number of periods
T = 25+1

k_t = np.zeros(T)
c_t = np.zeros(T)

# Solving the model forward: analytic solution
k_t[0] = k_0_check
c_t[0] = (1-α/(1+ρ)) * k_t[0]**α

for t in range(T-1):
    k_t[t+1] = (α/(1+ρ) * k_t[t]**α)/(1+n)/(1+g)
    c_t[t+1] = (1-α/(1+ρ)) * k_t[t+1]**α

# Store results for future use
k_simple, c_simple = k_t, c_t
In [11]:
# Number of periods
T = 25+1

k_t = np.zeros(T)
c_t = np.zeros(T)

# Solving the model forward: analytic solution
k_t[0] = k_0_check_up
c_t[0] = (1-α/(1+ρ)) * k_t[0]**α

for t in range(T-1):
    k_t[t+1] = (α/(1+ρ) * k_t[t]**α)/(1+n)/(1+g)
    c_t[t+1] = (1-α/(1+ρ)) * k_t[t+1]**α

# Store results for future use
k_simple_up, c_simple_up = k_t, c_t
In [12]:
# Convergence plot

plt.hlines(k_star, 0, T-1, color='C3', lw=1, linestyle='--')
plt.hlines(c_star, 0, T-1, color='C0', lw=1, linestyle='--')

plt.plot(k_simple, color='C3', lw=2, label='$k$')
plt.plot(c_simple, color='C0', lw=2, label='$c$')

plt.title('Convergence to steady state over time')

plt.xlabel('Period')
plt.xlim(0, T-1)

plt.yticks([c_star, k_star], ['$c^{*}$', '$k^{*}$'])
plt.ylim(0, 0.3)

plt.legend(loc='lower right', frameon=False)

plt.show()
In [13]:
# Phase diagram

kk = np.linspace(0, 1.5*k_star, 1000)
cc = np.linspace(0, 1.8*c_star, 1000)

k_GR = (α/(δ+n+g+n*g))**(1/(1-α))

plt.hlines(c_star, 0, k_star, 'C0', lw=1, linestyle='--')
plt.vlines(k_GR, 0, Constant_k(k_GR, params_analytic), 'C6', lw=1, linestyle='--')

plt.plot(k_check, c_check, 'C2', lw=3)
plt.plot(k_simple, c_simple, 'k--', lw=1)

plt.plot(kk, Constant_k(kk, params_analytic), 'C0', lw=2, label='$k_{t+1}=k_{t}$')
plt.plot(kk**0 * k_star, cc, 'C3', lw=2, label='$c_{t+1}=c_{t}$')

plt.plot(k_check_up, c_check_up, 'C2', lw=3, label='Transition path (numerical)')
plt.plot(k_simple_up, c_simple_up, 'k--', lw=1, label='Transition path (analytical)')

plt.plot(k_star, c_star, 'ko', label='Steady state')

plt.title('Phase diagram of the Ramsey model')
plt.xlabel('Capital per worker $k$')
plt.ylabel('Consumption per worker $c$')
plt.legend(loc='upper left', frameon=False)

plt.xlim(0, kk[-1])
plt.ylim(0, cc[-1])

plt.xticks([0, k_star, k_GR, kk[-1]], [0, '$k^{*}$', '$k^{*}_{GR}$', '$k$'])
plt.yticks([c_star, cc[-1]], ['$c^{*}$', '$c$'])

plt.show()

General case¶

In [14]:
# General case

# Parameters (for g>0 variables become per effective labor)
α = 0.33
ρ = 0.04
σ = 2
δ = 0.1
n = 0.01
g = 0.02

params = α, ρ, σ, δ, n, g

# Steady state
print(SteadyState(params))
k_star, c_star = SteadyState(params)

(2.2297027166708214, 1.0126239481814872)
In [15]:
# Convergence from below: Find the function minimum, starting from an initial guess

k_0 = k_star / 200

c_0_guess = Constant_k(k_0, params) - 0.01

result = opt.minimize(Path_crit, c_0_guess, args=(k_0, params, 100), method='Powell')
print(result)

c_0 = result.x

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 4.406202981709587e-06
       x: [ 1.725e-01]
     nit: 6
   direc: [[ 5.521e-04]]
    nfev: 261
In [16]:
# Convergence from above: Find the function minimum, starting from an initial guess

k_0_up = 2 * k_star

c_0_guess_up = Constant_k(k_0_up, params) + 0.01

result = opt.minimize(Path_crit, c_0_guess_up, args=(k_0_up, params, 100), method='Powell')
print(result)

c_0_up = result.x

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 1.7536952922686197e-06
       x: [ 1.367e+00]
     nit: 6
   direc: [[-2.514e-03]]
    nfev: 270
In [17]:
# Phase diagram (realistic params)

k_t, c_t = Path(c_0, k_0, params, 100)
k_t_up, c_t_up = Path(c_0_up, 2*k_star, params, 100)

kk = np.linspace(0, 1.5*k_star, 1000)
cc = np.linspace(0, 1.5*c_star, 1000)

plt.hlines(c_star, 0, k_star, 'C0', lw=1, linestyle='--')

plt.plot(k_t, c_t, 'C2', lw=2)
plt.plot(kk, Constant_k(kk, params), 'C0', lw=2, label='$\hat{k}_{t+1}=\hat{k}_{t}$')
plt.plot(kk**0 * k_star, cc, 'C3', lw=2, label='$\hat{c}_{t+1}=\hat{c}_{t}$')
plt.plot(k_t_up, c_t_up, 'C2', lw=2, label='Transition path')
plt.plot(k_star, c_star, 'ko', label='Balanced Growth Path')

plt.title('Phase diagram of the Ramsey model')
plt.xlabel('Capital per effective labor $\hat{k}$')
plt.ylabel('Consumption per effective labor $\hat{c}$')
plt.legend(loc='upper left', frameon=False)

plt.xlim(0, kk[-1])
plt.ylim(0, cc[-1])

plt.xticks([0, k_star, kk[-1]], [0, '$\hat{k}^{*}$', '$\hat{k}$'])
plt.yticks([c_star, cc[-1]], ['$\hat{c}^{*}$', '$\hat{c}$'])

plt.show()
In [18]:
# Streamplot of alternative paths

kk_g, cc_g = np.meshgrid(kk, cc)

def Δk(k, c):
    return k**α-(δ+n+g+n*g)*k-c

def Δc(k, c):
    return (((1+α*k**(α-1)-δ)/(1+ρ)/(1+n))**(1/σ)/(1+g)-1)*c

plt.streamplot(kk_g, cc_g, Δk(kk_g, cc_g), Δc(kk_g, cc_g), color='C7', density=0.7)

plt.plot(k_t, c_t, 'C2', lw=2)
plt.plot(kk, Constant_k(kk, params), 'C0', lw=2)
plt.plot(kk**0 * k_star, cc, 'C3', lw=2)
plt.plot(k_t_up, c_t_up, 'C2', lw=2)
plt.plot(k_star, c_star, 'ko')

plt.title('Streamplot of alternative paths')

plt.xlim(0, kk[-1])
plt.ylim(0, cc[-1])

plt.xticks([0, k_star, kk[-1]], [0, '$\hat{k}^{*}$', '$\hat{k}$'])
plt.yticks([c_star, cc[-1]], ['$\hat{c}^{*}$', '$\hat{c}$'])

plt.show()
In [19]:
s_t = 1 - c_t/ProductionFunction(k_t, params)
s_star = 1-c_star/ProductionFunction(k_star, params)

T = 100

plt.plot(k_t, lw=2)
plt.hlines(k_star, 0, T, lw=1, linestyle='--')
plt.title('Capital per effective labor $\hat{k}$')
plt.xlabel('Period')
plt.xlim(0, T)

plt.yticks([k_star], ['$\hat{k}^{*}$'])

plt.show()

plt.plot(c_t, lw=2)
plt.hlines(c_star, 0, T, lw=1, linestyle='--')
plt.title('Consumption per effective labor $\hat{c}$')
plt.xlabel('Period')
plt.xlim(0, T)

plt.yticks([c_star], ['$\hat{c}^{*}$'])

plt.show()

plt.plot(100*s_t, lw=2)
plt.hlines(100*s_star, 0, T, lw=1, linestyle='--')
plt.title('Saving rate $s$') # (%)')
plt.xlabel('Period')
plt.xlim(0, T)

plt.yticks([100*s_star], ['$s^{*}$'])

plt.show()
In [ ]: