Replicate CIA model from Chapter 4 of Monetary Theory and Policy, 2nd edition by Carl Walsh.
# Import numpy, pandas, linearsolve, scipy.optimize, matplotlib.pyplot
import numpy as np
import pandas as pd
import linearsolve as ls
from scipy.optimize import root,fsolve,broyden1,broyden2
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
alpha = 0.36
beta = 0.989
delta = 0.019
eta = 1
psi = 1.34
sigma = 2
A = 1
rhoa = 0.95
gamma = 0.8
phi=0.5
r_ss = 1/beta
yk_ss= 1/alpha*(r_ss-1+delta)
ck_ss = yk_ss-delta
def func(n):
'''Funciton to compute steady state labor'''
return (1-alpha)/psi*beta*yk_ss**((sigma-alpha)/(1-alpha))*ck_ss**(-sigma) - (1-n)**-eta*n**sigma
n_ss = root(func,0.3)['x'][0]
nk_ss = (yk_ss)**(1/(1-alpha))
k_ss = n_ss/nk_ss
y_ss = yk_ss*k_ss
c_ss = ck_ss*k_ss
m_ss = c_ss
a_ss = 1
u_ss = 1
pi_ss = 1
lam_ss = beta*c_ss**-sigma
mu_ss = (1-beta)*c_ss**-sigma
# Store steady state values in a list
ss = [a_ss,u_ss,m_ss,k_ss,pi_ss,r_ss,n_ss,c_ss,lam_ss,mu_ss,y_ss]
# Load parameter values into a Pandas Series
parameters = pd.Series({
'alpha':alpha,
'beta':beta,
'delta':delta,
'eta':eta,
'psi':psi,
'sigma':sigma,
'rhoa':rhoa,
'gamma':gamma,
'phi':phi,
'n_ss':n_ss,
'yk_ss':yk_ss,
'ck_ss':ck_ss
})
Solve and simulate the log-linearized model.
# Define function to compute equilibrium conditions
def equations(variables_forward,variables_current,parameters):
# Parameters
p = parameters
# Variables
fwd = variables_forward
cur = variables_current
# Household Euler equation
foc1 = p.alpha*cur.k+(1-p.alpha)*cur.n + fwd.a - cur.y
foc2 = p.ck_ss*fwd.m + fwd.k - (1-p.delta)*cur.k - p.yk_ss*cur.y
foc3 = p.alpha*p.yk_ss*(fwd.y - fwd.k) - cur.r
foc4 = fwd.lam + cur.r - cur.lam
foc5 = (1+p.eta*p.n_ss/(1-p.n_ss))*cur.n - cur.y - cur.lam
foc6 = cur.r + fwd.pi - cur.rn
foc7 = -p.sigma*fwd.maux-fwd.pi - cur.lam
foc8 = cur.m-cur.pi+cur.u - fwd.m
foc9 = cur.maux - fwd.m
foc10= p.gamma*cur.u+p.phi*fwd.a - fwd.u
foc11= p.rhoa*cur.a - fwd.a
# Stack equilibrium conditions into a numpy array
return np.array([
foc1,
foc2,
foc3,
foc4,
foc5,
foc6,
foc7,
foc8,
foc9,
foc10,
foc11,
])
# Initialize the model
model = ls.model(equations = equations,
n_states=4,
n_exo_states=3,
var_names=['a', 'u', 'm', 'k', 'lam', 'pi', 'rn', 'r', 'n', 'y','maux'],
parameters = parameters)
# Compute the steady state numerically
guess = 0*np.array([1,1,10,10,1,1,0.5,2,1,1,1])
model.compute_ss(guess,method='fsolve')
# Construct figure and axes
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
# Iterate over different degrees of persistence for money growth shock
for gamma in [0.5,0.8]:
model.parameters['gamma'] = gamma
# Solve the model
model.approximate_and_solve(log_linear=False)
# Compute impulse responses and plot
model.impulse(T=17,t0=1,shocks=None,percent=True)
# Plot
y = model.irs['e_u']['y']
n = model.irs['e_u']['n']
rn = model.irs['e_u']['rn']
pi = model.irs['e_u']['pi']
tme=rn.index
ax1.plot(tme,y,lw=5,alpha=0.5,label='y ($\gamma='+str(gamma)+'$)')
ax1.plot(tme,n,'--',lw=5,alpha=0.5,label='n ($\gamma='+str(gamma)+'$)')
ax1.grid(True)
ax1.legend(loc='lower right')
ax2.plot(tme,rn,lw=5,alpha=0.5,label='Rn ($\gamma='+str(gamma)+'$)')
ax2.plot(tme,pi,'--',lw=5,alpha=0.5,label='$\pi$ ($\gamma='+str(gamma)+'$)')
ax2.grid(True)
ax2.legend()
Approximate, solve, and simulate the log-linearized model.
# Define function to compute equilibrium conditions
def equations(variables_forward,variables_current,parameters):
# Parameters
p = parameters
# Variables
fwd = variables_forward
cur = variables_current
# Household Euler equation
foc_1 = cur.a**p.rhoa - fwd.a
foc_2 = cur.u**p.gamma*cur.a**p.phi - fwd.u
foc_3 = cur.lam+cur.mu - cur.c**-p.sigma
foc_4 = cur.lam*(1-p.alpha)*cur.y/cur.n - p.psi*(1-cur.n)**-p.eta
foc_5 = p.beta*(fwd.lam*cur.Rn)/fwd.pi - cur.lam
foc_6 = p.beta*(fwd.mu+fwd.lam)/fwd.pi - cur.lam
foc_7 = p.beta*(fwd.lam*(p.alpha*fwd.y/fwd.k+1-p.delta)) - cur.lam
foc_8 = cur.a*cur.k**alpha*cur.n**(1-p.alpha) - cur.y
foc_9 = cur.c+fwd.k-(1-p.delta)*cur.k - cur.y
foc_10 = cur.m/cur.pi*cur.u - fwd.m
foc_11 = fwd.m - cur.c
# Stack equilibrium conditions into a numpy array
return np.array([
foc_1,
foc_2,
foc_3,
foc_4,
foc_5,
foc_6,
foc_7,
foc_8,
foc_9,
foc_10,
foc_11
])
# Initialize the model
varNames=['a','u','m','k','pi','Rn','n','c','lam','mu','y']
parameters = parameters[['alpha','beta','delta','eta','psi','sigma','rhoa','gamma','phi']]
model = ls.model(equations = equations,
n_states=4,
n_exo_states=3,
var_names=varNames,
parameters = parameters)
# Set the steady state using exact values calculated above
model.set_ss(ss)
# Construct figure and axes
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
# Iterate over different degrees of persistence for money growth shock
for gamma in [0.5,0.8]:
model.parameters['gamma'] = gamma
# Find the log-linear approximation around the non-stochastic steady state
model.approximate_and_solve()
# Compute impulse responses and plot
model.impulse(T=17,t0=1,shocks=None,percent=True)
# Plot
y = model.irs['e_u']['y']
n = model.irs['e_u']['n']
rn = model.irs['e_u']['Rn']
pi = model.irs['e_u']['pi']
tme=rn.index
ax1.plot(tme,y,lw=5,alpha=0.5,label='y ($\gamma='+str(gamma)+'$)')
ax1.plot(tme,n,'--',lw=5,alpha=0.5,label='n ($\gamma='+str(gamma)+'$)')
ax1.grid(True)
ax1.legend(loc='lower right')
ax2.plot(tme,rn,lw=5,alpha=0.5,label='Rn ($\gamma='+str(gamma)+'$)')
ax2.plot(tme,pi,'--',lw=5,alpha=0.5,label='$\pi$ ($\gamma='+str(gamma)+'$)')
ax2.grid(True)
ax2.legend()