Here I demonstrate how how relatively straightforward it is to appoximate, solve, and simulate a DSGE model using linearsolve
. In the example that follows, I describe the procedure more carefully.
In the block of code that immediately follows, I input the model, solve for the steady state, compute the log-linear approximation of the equilibirum conditions, and compute some impulse responses following a shock to technology $A_t$.
# Import numpy, pandas, linearsolve, matplotlib.pyplot
import numpy as np
import pandas as pd
import linearsolve as ls
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
# Input model parameters
parameters = pd.Series(dtype=float)
parameters['alpha'] = .35
parameters['beta'] = 0.99
parameters['delta'] = 0.025
parameters['rhoa'] = .9
parameters['sigma'] = 1.5
parameters['A'] = 1
# Funtion that evaluates the equilibrium conditions
def equations(variables_forward,variables_current,parameters):
# Parameters
p = parameters
# Variables
fwd = variables_forward
cur = variables_current
# Household Euler equation
euler_eqn = p.beta*fwd.c**-p.sigma*(p.alpha*fwd.a*fwd.k**(p.alpha-1)+1-p.delta) - cur.c**-p.sigma
# Goods market clearing
market_clearing = cur.c + fwd.k - (1-p.delta)*cur.k - cur.a*cur.k**p.alpha
# Exogenous technology
technology_proc = p.rhoa*np.log(cur.a) - np.log(fwd.a)
# Stack equilibrium conditions into a numpy array
return np.array([
euler_eqn,
market_clearing,
technology_proc
])
# Initialize the model
model = ls.model(equations = equations,
n_states=2,
n_exo_states = 1,
var_names=['a','k','c'],
parameters = parameters)
# Compute the steady state numerically
guess = [1,1,1]
model.compute_ss(guess)
# Find the log-linear approximation around the non-stochastic steady state and solve
model.approximate_and_solve()
# Compute impulse responses and plot
model.impulse(T=41,t0=5,shocks=None)
fig = plt.figure(figsize=(12,4))
ax1 =fig.add_subplot(1,2,1)
model.irs['e_a'][['a','k','c']].plot(lw='5',alpha=0.5,grid=True,ax=ax1).legend(loc='upper right',ncol=3)
ax2 =fig.add_subplot(1,2,2)
model.irs['e_a'][['e_a','a']].plot(lw='5',alpha=0.5,grid=True,ax=ax2).legend(loc='upper right',ncol=2)
<matplotlib.legend.Legend at 0x7feb88af4280>