Example based on Diaz-Gimenez (2001) - Linear Quadratic Approximations: An Introduction (Chapter 2 in 'Computational Methods for the Study of Dynamic Economies', edited by Marimon & Scott)
This model is also used by Aldrich, Fernandez-Villaverde, Gallant, & Rubio-Ramirez (2011) - "Tapping the supercomputer under your desk: Solving dynamic equilibrium models with graphics processors". But they do use slightly different parameters to those used here.
The Stochastic Neo-Classical Growth Model was first developed in Brock & Mirman (1972), although our treatment of the model here is based on Díaz-Giménez (2001). If you are unfamiliar with the model a full description can be found in Section 2.2 of Díaz-Giménez (2001); he also briefly discusses the relation to the social planners problem, how the problem looks in both sequential and recursive formulation, and how we know that value function iteration will give us the correct solution. In what follows I assume you are familiar with these issues and simply give the value function problem. Our concentration here is on how to solve this problem using the toolkit.
The value function problem to be solved is, $$ V(k,z)=\sup_{k'} \{ \frac{c^{1-\gamma}}{1-\gamma}+\beta*E[V(k',z')|z] \} $$ subject to $$ c+i=exp(z) k^{\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, $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 first thing we need to do is to call some packages that are needed for some statistical operations (eg. the Tauchen Method uses some functions relating to the normal distribution) and for graphs.
In [2]:
using Distributions #The Tauchen Method uses the Distributions package to evaluate the normal cdf
using PyPlot #Needed for graphing
# Some options that can be passed to various commands. We just let them use the defaults in this example. This section
# tells you what those (unused) defaults are.
# The following options are used by a number of commands
# Tolerance=10^(-9.0) #This is the degree of 'accuracy' that various commands will aim to acheive
# Parallel=0 #If =0, no parallelization, =1 parallel on cpu, =2 parallel on gpu (where available)
# Verbose=0 #If =1, gives feedback/output on the (internal) progress of the codes
# The following options are only used by the value fn iteration codes
# Howards=80 #Number of iterations used by Howards Improvement Algorithm (part of value fn iteration). Default is 80.
# PolIndOrVal=1 #If =1, PolicyIndexes will contain the indexes of the optimal policy;
#If =2, PolicyIndexes will contain the values of the optimal policy
# LowMemory=0 #If =0 it solves the value function iteration in a way that is faster, but very memory intensive.
#If =1, slower but less memory.
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 [9]:
# there is certainly a better way of including all the scripts as a bunch, but I am not yet sure how to do this
include("tauchenmethod.jl")
include("valuefniter_case1.jl")
include("createreturnfnmatrix_case1.jl")
include("valuefniter_case1_nod_raw.jl")
Out[9]:
Now we declare various paramaters and grids.
In [9]:
# The sizes of the grids
n_z=4;
n_k=512;
# Discounting rate
beta = 0.96;
# Give the parameter values
alpha = 0.33;
gamma=1; #gamma=1 is log-utility
rho = 0.95;
delta = 0.10;
sigmasq_epsilon=0.09;
In [10]:
#If you set this parameter to 0 then the parameters will all be set to those used by Aldrich, Fernandez-Villaverde, Gallant, & Rubio-Ramirez (2011)
Javier=1;
if Javier==0
n_z=4;
beta=0.984;
gamma=2;
alpha=0.35;
delta=0.01;
rho=0.95;
sigma_epsilon=0.005;
sigmasq_epsilon=sigma_epsilon^2;
end
In [11]:
# Step 1: Compute the steady state
K_ss=((alpha*beta)/(1-beta*(1-delta)))^(1/(1-alpha));
X_ss= delta*K_ss;
#These are not really needed; we just use them to determine the grid on
#capital. I mainly calculate them to stay true to original article.
Out[11]:
In [12]:
## Create grids (grids are defined as a column vectors)
k_grid=linspace(0,20*K_ss,n_k); # Grids should always be declared as column vectors
# We use the Tauchen Method to define the grid for the exogenous technology shock, z.
q=3; # A parameter needed for the Tauchen Method
tic()
z_grid, pi_z =tauchenmethod(0,sigmasq_epsilon,rho,n_z,q,0,0); #[states, transmatrix]=TauchenMethod_Param(mew,sigmasq,rho,znum,q), transmatix is (z,zprime)
time=toc()
Out[12]:
We now define the return function.
In [16]:
#Now create the return function
#StochasticNeoClassicalGrowthModel_ReturnFn
function returnfn(kprime,k,z)#,gamma,alpha,delta)
F=-Inf;
c=exp(z)*k^alpha-(kprime-(1-delta)*k);
if c[1]>0 #Must use c[1] rather than c as otherwise Julia makes c a 1x1 array, and is then unable to do a 'lessthan' operation on the 1x1 array.
if gamma==1
F=log(c);
else
F=(c^(1-gamma))/(1-gamma);
end
end
return F
end
Out[16]:
In terms of using this Julia library to solve value function problems, we have now done all of the work (declaring 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 some other models you may also want to declare some further decision variables, other than just the endogenous state]
In [17]:
# Solve
# Do the value function iteration. Returns both the value function itself, and the optimal policy function.
d_grid=0; #no d variable
n_d=0;
V0=ones(n_k,n_z);
tic()
V, Policy = valuefniter_case1(V0, n_d,n_k,n_z,d_grid,k_grid,z_grid, pi_z, beta, returnfn);
toc()
Out[17]:
Done! The model is now solved.
Let's finish by graphing the value function to see what it looks like.
In [15]:
#I don't yet know how to do 3d plots in Julia, so let's just do a 2d plot for each of the z values.
plot(k_grid,V[:,1],k_grid,V[:,2],k_grid,V[:,3],k_grid,V[:,4])
Out[15]:
Finished!
A couple of comments on how this compares to doing the same thing in Matlab.
Writing the code for Julia was just as easy as writing in for Matlab. The run times on my laptop are slightly faster for Julia than Matlab, despite this being an almost direct copy of code I orginially wrote for Matlab, and which is optimized for speed in Matlab (in that everything is vectorized. I wrote many different variants of the valuefniter_case1 code in Matlab before I settled on the one implemented here as being the fastest in Matlab).
The main difficulty with using Julia was that debugging and the code editor are not as mature and easy to use as in Matlab. For example, using iJulia notebooks the error messages tell me the relevant line of code, but the iJulia cells do not tell me line numbers, so I had to manually count them. Also, you can't mouse-over end's to see which if-for-while they relate to.
One issue that you need to know is that Julia treats 1x1 array's very differently to scalars. This often caused me issues with the return function. For example, c in returnfn is considered a 1x1 array. Doing 'lessthan' on a 1x1 array generates an error. Thus the code says 'c[1]>0' instead of just 'c>0'. You have to be careful using '^' in the return function for the same reason; better to use '.^'.
There are also some other minor issues with Julia. For example I had to edit commands where I normally use length(a) as an input [eg. zeros(length(a),1) or for i=1:length(a)] as Julia would not accept length(a) as a argument. I am not sure if these reflect design choices or are points that will be ironed out as Julia matures.
Also, if you come from coding Matlab, just be aware that in Julia if you set b=a, and then modify b (eg. b[2]=3), you will also modify a. You need to set b=copy(a) to avoid this.