<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A Cash-in-Advance Model
Replicate CIA model from Chapter 4 of Monetary Theory and Policy, 2nd edition by Carl Walsh.
Step1:  Set parameters
Step2:  Compute exact steady state
Step3:  Linear model
Solve and simulate the log-linearized model.
Step4:  Nonlinear model
Approximate, solve, and simulate the log-linearized model.
<ASSISTANT_TASK:>
Python Code:
# 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
Explanation: A Cash-in-Advance Model
Replicate CIA model from Chapter 4 of Monetary Theory and Policy, 2nd edition by Carl Walsh.
End of explanation
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
Explanation: Set parameters
End of explanation
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
})
Explanation: Compute exact steady state
End of explanation
# 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()
Explanation: Linear model
Solve and simulate the log-linearized model.
End of explanation
# 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()
Explanation: Nonlinear model
Approximate, solve, and simulate the log-linearized model.
End of explanation
<END_TASK>