Consider the equilibrium conditions for a basic RBC model without labor:
\begin{align} C_t^{-\sigma} & = \beta E_t \left[C_{t+1}^{-\sigma}(\alpha A_{t+1} K_{t+1}^{\alpha-1} + 1 - \delta)\right]\\ Y_t & = A_t K_t^{\alpha}\\ I_t & = K_{t+1} - (1-\delta)K_t\\ Y_t & = C_t + I_t\\ \log A_t & = \rho_a \log A_{t-1} + \epsilon_t \end{align}In the nonstochastic steady state, we have:
\begin{align} K & = \left(\frac{\alpha A}{1/\beta+\delta-1}\right)^{\frac{1}{1-\alpha}}\\ Y & = AK^{\alpha}\\ I & = \delta K\\ C & = Y - I \end{align}Given values for the parameters $\beta$, $\sigma$, $\alpha$, $\delta$, and $A$, steady state values of capital, output, investment, and consumption are easily computed.
# 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
linearsolve
¶To initialize the model, we need to first set the model's parameters. We do this by creating a Pandas Series variable called parameters
:
# 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
Next, we need to define a function that returns the equilibrium conditions of the model. The function will take as inputs two vectors: one vector of "current" variables and another of "forward-looking" or one-period-ahead variables. The function will return an array that represents the equilibirum conditions of the model. We'll enter each equation with all variables moved to one side of the equals sign. For example, here's how we'll enter the produciton fucntion:
production_function = technology_current*capital_current**alpha - output_curent
Here the variable production_function
stores the production function equation set equal to zero. We can enter the equations in almost any way we want. For example, we could also have entered the production function this way:
production_function = 1 - output_curent/technology_current/capital_current**alpha
One more thing to consider: the natural log in the equation describing the evolution of total factor productivity will create problems for the solution routine later on. So rewrite the equation as:
\begin{align} A_{t+1} & = A_{t}^{\rho_a}e^{\epsilon_{t+1}}\\ \end{align}So the complete system of equations that we enter into the program looks like:
\begin{align} C_t^{-\sigma} & = \beta E_t \left[C_{t+1}^{-\sigma}(\alpha Y_{t+1} /K_{t+1}+ 1 - \delta)\right]\\ Y_t & = A_t K_t^{\alpha}\\ I_t & = K_{t+1} - (1-\delta)K_t\\ Y_t & = C_t + I_t\\ A_{t+1} & = A_{t}^{\rho_a}e^{\epsilon_{t+1}} \end{align}Now let's define the function that returns the equilibrium conditions:
# 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
euler_eqn = p.beta*fwd.c**-p.sigma*(p.alpha*fwd.y/fwd.k+1-p.delta) - cur.c**-p.sigma
# Production function
production_fuction = cur.a*cur.k**p.alpha - cur.y
# Capital evolution
capital_evolution = fwd.k - (1-p.delta)*cur.k - cur.i
# Goods market clearing
market_clearing = cur.c + cur.i - cur.y
# Exogenous technology
technology_proc = cur.a**p.rhoa- fwd.a
# Stack equilibrium conditions into a numpy array
return np.array([
euler_eqn,
production_fuction,
capital_evolution,
market_clearing,
technology_proc
])
Notice that inside the function we have to define the variables of the model form the elements of the input vectors variables_forward
and variables_current
.
To initialize the model, we need to specify the total number of state variables in the model, the number of state variables with exogenous shocks, the names of the endogenous variables, and the parameters of the model.
It is essential that the variable names are ordered in the following way: First the names of the endogenous variables with the state variables with exogenous shocks, then the state variables without shocks, and finally the control variables. Ordering within the groups doesn't matter.
# Initialize the model
rbc = ls.model(equations = equations,
n_states=2,
n_exo_states=1,
var_names=['a','k','c','y','i'],
parameters=parameters)
Next, we need to compute the nonstochastic steady state of the model. The .compute_ss()
method can be used to compute the steady state numerically. The method's default is to use scipy's fsolve()
function, but other scipy root-finding functions can be used: root
, broyden1
, and broyden2
. The optional argument options
lets the user pass keywords directly to the optimization function. Check out the documentation for Scipy's nonlinear solvers here: http://docs.scipy.org/doc/scipy/reference/optimize.html
# Compute the steady state numerically
guess = [1,1,1,1,1]
rbc.compute_ss(guess)
print(rbc.ss)
a 1.000000 k 34.398226 c 2.589794 y 3.449750 i 0.859956 dtype: float64
Note that the steady state is returned as a Pandas Series. Alternatively, you could compute the steady state directly and then sent the rbc.ss
attribute:
# Steady state solution
p = parameters
K = (p.alpha*p.A/(1/p.beta+p.delta-1))**(1/(1-p.alpha))
C = p.A*K**p.alpha - p.delta*K
Y = p.A*K**p.alpha
I = Y - C
rbc.set_ss([p.A,K,C,Y,I])
print(rbc.ss)
a 1.000000 k 34.398226 c 2.589794 y 3.449750 i 0.859956 dtype: float64
Now we use the .log_linear()
method to find the log-linear appxoximation to the model's equilibrium conditions. That is, we'll transform the nonlinear model into a linear model in which all variables are expressed as log-deviations from the steady state. Specifically, we'll compute the matrices $A$ and $B$ that satisfy:
where the vector $x_{t}$ denotes the log deviation of the endogenous variables from their steady state values.
# Find the log-linear approximation around the non-stochastic steady state
rbc.log_linear_approximation()
print('The matrix A:\n\n',np.around(rbc.a,4),'\n\n')
print('The matrix B:\n\n',np.around(rbc.b,4))
The matrix A: [[ 0.00000e+00 -8.30000e-03 -3.59900e-01 8.30000e-03 0.00000e+00] [ 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00] [ 0.00000e+00 3.43982e+01 0.00000e+00 0.00000e+00 0.00000e+00] [ 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00] [-1.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00]] The matrix B: [[-0. -0. -0.3599 -0. -0. ] [-3.4497 -1.2074 -0. 3.4497 -0. ] [-0. 33.5383 -0. -0. 0.86 ] [-0. -0. -2.5898 3.4497 -0.86 ] [-0.9 -0. -0. -0. -0. ]]
Finally, we need to obtain the solution to the log-linearized model. The solution is a pair of matrices $F$ and $P$ that specify:
We use the .klein()
method to find the solution.
# Solve the model
rbc.solve_klein(rbc.a,rbc.b)
# Display the output
print('The matrix F:\n\n',np.around(rbc.f,4),'\n\n')
print('The matrix P:\n\n',np.around(rbc.p,4))
The matrix F: [[ 0.2297 0.513 ] [ 1. 0.35 ] [ 3.3197 -0.1408]] The matrix P: [[0.9 0. ] [0.083 0.9715]]
One the model is solved, use the .impulse()
method to compute impulse responses to exogenous shocks to the state. The method creates the .irs
attribute which is a dictionary with keys equal to the names of the exogenous shocks and the values are Pandas DataFrames with the computed impulse respones. You can supply your own values for the shocks, but the default is 0.01 for each exogenous shock.
# Compute impulse responses and plot
rbc.impulse(T=41,t0=1,shocks=None,percent=True)
print('Impulse responses to a 0.01 unit shock to A:\n\n',rbc.irs['e_a'].head())
Impulse responses to a 0.01 unit shock to A: e_a a k c y i 0 0.0 0.000 0.000000 0.000000 0.000000 0.000000 1 1.0 1.000 0.000000 0.229718 1.000000 3.319739 2 0.0 0.900 0.082993 0.249318 0.929048 2.976083 3 0.0 0.810 0.155321 0.265744 0.864362 2.667127 4 0.0 0.729 0.218116 0.279348 0.805341 2.389390
Plotting is easy.
rbc.irs['e_a'][['a','k','c','y','i']].plot(lw='5',alpha=0.5,grid=True).legend(loc='upper right',ncol=2)
rbc.irs['e_a'][['e_a','a']].plot(lw='5',alpha=0.5,grid=True).legend(loc='upper right',ncol=2)
<matplotlib.legend.Legend at 0x7fe197b9db80>
Creating a stochastic simulation of the model is straightforward with the .stoch_sim()
method. In the following example, I create a 151 period (including t=0) simulation by first simulating the model for 251 periods and then dropping the first 100 values. The standard deviation of the shock to $A_t$ is set to 0.00763. The seed for the numpy random number generator is set to 0.
rbc.stoch_sim(T=121,drop_first=100,cov_mat=np.array([0.00763**2]),seed=0,percent=True)
rbc.simulated[['k','c','y','i']].plot(lw='5',alpha=0.5,grid=True).legend(loc='upper right',ncol=4)
rbc.simulated[['a']].plot(lw='5',alpha=0.5,grid=True).legend(ncol=4)
rbc.simulated['e_a'].plot(lw='5',alpha=0.5,grid=True).legend(ncol=4)
<matplotlib.legend.Legend at 0x7fe197fb7790>