The Basic Real Business Cycle Model presented here simply extends the Stochastic Neo-Classical Growth Model to include an endogenous choice of labour supply (of how much to work). This is an important difference in terms of the Toolkit we know have a '$d$' variable: a choice variable that does not affect tomorrows state (in the Stochastic Neo-Classical Growth Model the only choice variable was '$aprime$', next periods state). The following code solves the model and then uses the results to reproduce the relevant findings of Aruoba, Fernandez-Villaverde, and Rubio-Ramirez (2006) looking at the sizes of numerical errors --- this serves the dual purpose of allowing you to see how measures of the numerical errors change with the choice of grids, and of providing an example of how to use the policy function created by the Value Function (Case 1) command.
The value function problem to be solved is, $$ V(k,z)=\sup_{k'} \{ \frac{(c^{\theta} (1-l)^{1-\theta})^{\tau}}{1-\tau}+\beta*E[V(k',z')|z] $$ subject to $$ c+i=exp(z) k^{\alpha} l^{1-\alpha} $$ $$ k'=(1-\delta)k+i $$ $$ z'=\rho z +\epsilon' \quad ,\epsilon \overset{iid}{\sim} N(0,\sigma_{\epsilon}^2) $$ where $k$ is physical capital, $i$ is investment, $c$ is consumption, $l$ is labour supply, $z$ is a productivity shock that follows an AR(1) process; $exp(z) k^{\alpha}$ is the production function, $\delta$ is the depreciation rate.
We will use the Tauchen Method to discretize the AR(1) shock.
The following code solves this model and then reproduces a number of relevant elements from Tables & Figures in Aruoba, Fernandez-Villaverde, & Rubio-Ramirez (2006) to the accuracy of numerical solutions; it contains many comments describing what is being done at each step. [Aruoba, Fernandez-Villaverde, & Rubio-Ramirez (2006) discuss solving this model with a wide variety of different numerical solution methods.]
In [1]:
using Distributions #The Tauchen Method uses the Distributions package to evaluate the normal cdf
using PyPlot #Needed for graphing
Next we include some general julia scripts that implement the Tauchen method, and perform value function iteration with a discrete state space (no interpolation at all).
If you want to see the contents of the scripts you can check them out on github.com/robertdkirkby/economics-with-julia
In [2]:
#Currently I just pull in all these scripts one-by-one.
#I intend to make some way to pull in all the relevant scripts in one go but am not yet sure how to do this in Julia
include("tauchenmethod.jl")
include("valuefniter_case1.jl")
include("valuefniter_case1_raw.jl")
include("createreturnfnmatrix_case1.jl")
Out[2]:
Define the parameters
In [3]:
UseAlternativeParams=1
#Discounting rate
beta = 0.9896
#Parameter values
alpha = 0.4
theta=0.357
rho = 0.95
tau=2
delta = 0.0196
sigmasq_epsilon=0.007^2
if UseAlternativeParams==1
tau=50
sigmasq_epsilon=0.035^2
end
Out[3]:
Define the grids (it is important that grids should be defined as column vectors)
In [4]:
n_d=10
n_a=100
n_z=7 #Note, for the figures to correctly use z=0 this must be an odd number (make it 39)
q=3 #Parameter for the Tauchen method
#Aruoba, Fernandez-Villaverde, & Rubio-Ramirez (2006) use 40 points for z and 25000 points for a (they use functional maximization so have no d grid).
#Step 1: Compute the steady state (just use these when picking grids)
varphi=(1/alpha*(1/beta-1+delta))^(1/(1-alpha))
Psi=theta/(1-theta)*(1-alpha)*varphi^(-alpha)
Omega=varphi^(1-alpha)-delta
K_ss=Psi/(Omega+varphi*Psi)
tic()
z_grid, pi_z=tauchenmethod(0,sigmasq_epsilon,rho,n_z,q) #transition matrix pi_z is z-by-zprime
time0=toc()
a_grid=linspace(0.01,5*K_ss,n_a)
d_grid=linspace(0,1,n_d)
println("Time to run Tauchen Method: $time0")
We now define the return function.
In [5]:
function returnfn(d,aprime,a,z)
#Note that you can parameters that have been previously defined inside the return function.
# d is decision variables that do not directly affect next state; here hours worked l
# aprime is next period k
# a is this period k
# z is z
F=-Inf;
c=exp(z)*(a^alpha)*(d^(1-alpha))-(aprime-(1-delta)*a);
if c[1]>0
F=(((c^theta)*((1-d).^(1-theta))).^(1-tau))/(1-tau);
end
return F
end
Out[5]:
In terms of using this Julia library to solve value function problems, we have now done all of the work (declaring grids for the decision variables, grids for the endogenous state, grids and transition matrix for the exogenous state, declaring the return function, and declaring the discount rate (beta)). [The parameters are needed as part of declaring these objects]
In [6]:
tic()
V0=zeros(n_a,n_z)
V,Policy =valuefniter_case1(V0, n_d,n_a,n_z,d_grid,a_grid,z_grid, pi_z, beta, returnfn)
toc()
Out[6]:
Done! The model is now solved.
Let's finish by graphing the value function to see what it looks like.
In [7]:
z0=ceil(n_z/2) #Just graph V for the 'middle' z value
plot(a_grid[10:n_a],V[10:n_a;z0])
# I start from the 10th point in a_grid as V approaches minus infinity near capital=0,
# so if the graph includes that section it just looks a bit silly.
Out[7]:
I plan to update this notebook later on to replicate a number of the figures and tables in Aruoba, Fernandez-Villaverde & Rubio-Ramirez (2006)