In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import linearsolve as ls
np.set_printoptions(suppress=True)
from scipy.optimize import root,fsolve,broyden1,broyden2
%matplotlib inline
In [2]:
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
In [3]:
r_ss = 1/beta
yk_ss= 1/alpha*(r_ss-1+delta)
ck_ss = yk_ss-delta
def func(n):
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
In [4]:
ss = [a_ss,u_ss,m_ss,k_ss,pi_ss,r_ss,n_ss,c_ss,lam_ss,mu_ss,y_ss]
In [5]:
def equilibrium_equations(vars_fwd,vars_cur,parameters):
# Parameters
p = parameters
# Define the variables: predetermined or state variables come first
a_fwd, u_fwd, m_fwd, k_fwd, lam_fwd, pi_fwd, rn_fwd, r_fwd, n_fwd, y_fwd, maux_fwd = vars_fwd
a_cur, u_cur, m_cur, k_cur, lam_cur, pi_cur, rn_cur, r_cur, n_cur, y_cur, maux_cur = vars_cur
# Household Euler equation
foc1 = [y_cur,p.alpha*k_cur+(1-p.alpha)*n_cur + a_fwd]
foc2 = [p.yk_ss*y_cur,p.ck_ss*m_fwd + k_fwd-(1-p.delta)*k_cur]
foc3 = [r_cur,p.alpha*p.yk_ss*(y_fwd-k_fwd)]
foc4 = [lam_cur,lam_fwd+r_cur]
foc5 = [y_cur+lam_cur,(1+p.eta*p.n_ss/(1-p.n_ss))*n_cur]
foc6 = [rn_cur, r_cur+pi_fwd]
foc7 = [lam_cur,-p.sigma*maux_fwd-pi_fwd]
foc8 = [m_fwd,m_cur-pi_cur+u_cur]
foc9 = [m_fwd,maux_cur]
foc10= [u_fwd,p.gamma*u_cur+p.phi*a_fwd]
foc11= [a_fwd,p.rhoa*a_cur]
# Stack equilibrium conditions into a numpy array
return np.array([
foc1,
foc2,
foc3,
foc4,
foc5,
foc6,
foc7,
foc8,
foc9,
foc10,
foc11,
])
# Initialize the model
parameter_vals = [alpha,beta,delta,eta,psi,sigma,rhoa,gamma,phi,n_ss,yk_ss,ck_ss]
parameter_names = ['alpha','beta','delta','eta','psi','sigma','rhoa','gamma','phi','n_ss','yk_ss','ck_ss']
parameters = pd.Series(parameter_vals,index = parameter_names)
model = ls.model(equations = equilibrium_equations,
nstates=4,
varNames=['a', 'u', 'm', 'k', 'lam', 'pi', 'rn', 'r', 'n', 'y',' maux'],
shockNames=['ea','eu','em','ek'],
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')
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
for gamma in [0.5,0.8]:
model.parameters['gamma'] = gamma
# Find the log-linear approximation around the non-stochastic steady state
model.log_linear(islinear=True)
# Solve the model
model.solve_klein(model.a,model.b)
# Compute impulse responses and plot
model.impulse(T=17,t0=1,shock=None,percent=True)
y = model.irs['eu']['y']
n = model.irs['eu']['n']
rn = model.irs['eu']['rn']
pi = model.irs['eu']['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()
In [6]:
def equilibrium_equations(vars_fwd,vars_cur,parameters):
# Define parameters
p = parameters
# Define the variables: predetermined or state variables come first
a_fwd, u_fwd, m_fwd, k_fwd, pi_fwd, Rn_fwd, n_fwd, c_fwd, lam_fwd, mu_fwd, y_fwd = vars_fwd
a_cur, u_cur, m_cur, k_cur, pi_cur, Rn_cur, n_cur, c_cur, lam_cur, mu_cur, y_cur = vars_cur
# Household Euler equation
foc_1 = [a_fwd,a_cur**p.rhoa]
foc_2 = [u_fwd,u_cur**p.gamma*a_cur**p.phi]
foc_3 = [c_cur**-p.sigma, lam_cur+mu_cur]
foc_4 = [p.psi*(1-n_cur)**-p.eta,lam_cur*(1-p.alpha)*y_cur/n_cur]
foc_5 = [lam_cur, p.beta*(lam_fwd*Rn_cur)/pi_fwd]
foc_6 = [lam_cur, p.beta*(mu_fwd+lam_fwd)/pi_fwd]
foc_7 = [lam_cur, p.beta*(lam_fwd*(p.alpha*y_fwd/k_fwd+1-p.delta))]
foc_8 = [y_cur,a_cur*k_cur**alpha*n_cur**(1-p.alpha)]
foc_9 = [y_cur,c_cur+k_fwd-(1-p.delta)*k_cur]
foc_10 = [m_fwd,m_cur/pi_cur*u_cur]
foc_11 = [c_cur,m_fwd]
# 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']
shockNames=['ea','eu','em','ek']
parameter_vals = [alpha,beta,delta,eta,psi,sigma,rhoa,gamma,phi]
parameter_names = ['alpha','beta','delta','eta','psi','sigma','rhoa','gamma','phi']
parameters = pd.Series(parameter_vals,index = parameter_names)
model = ls.model(equations = equilibrium_equations,
nstates=4,
varNames=varNames,
shockNames=shockNames,
parameters = parameters)
guess = ss
# model.compute_ss(guess,method='fsolve')
model.set_ss(ss)
# model.set_ss(ss)
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)
for gamma in [0.5,0.8]:
model.parameters['gamma'] = gamma
# Find the log-linear approximation around the non-stochastic steady state
model.log_linear()
# Solve the model
model.solve_klein(model.a,model.b)
# Compute impulse responses and plot
model.impulse(T=17,t0=1,shock=None,percent=True)
y = model.irs['eu']['y']
n = model.irs['eu']['n']
rn = model.irs['eu']['Rn']
pi = model.irs['eu']['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()