{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Projection implementation of stochastic growth model adapted from Fabrice Collard's Matlab code, http://fabcol.free.fr/\n",
"# tested in Julia 0.6.1.1\n",
"# this code is part of chapter 5, \"Advanced Numerical Techniques\" from the book: \"Introduction to Quantitative Macroeconomics using Julia\"\n",
"# Academic Press - Elsevier\n",
"# for comments, email at: petre(dot)caraiani(at)gmail(dot)com\n",
"\n",
"# The Optimal Growth Model\n",
"# Collocation method (Markov Chain case) \n",
"\n",
"global kmin, ksup, XX, kt;\n",
"global nstate, nbk, ncoef, XX, XT, PI;\n",
"\n",
"# Parameters \n",
"nbk = 4; # Degree of polynomials (capital)\n",
"nodes = nbk + 1; # Nodes\n",
"nstate= 2; \n",
"ncoef = nbk +1 ; # of coefficients\n",
" \n",
"# Structural Parameters \n",
"delta = 0.1;\n",
"beta = 0.95;\n",
"alpha = 0.3;\n",
"sigma = 1.5;\n",
"ysk =(1-beta*(1-delta))/(alpha*beta);\n",
"ksy = 1/ysk;\n",
"ys = ksy^(alpha/(1-alpha));\n",
"ks = ys^(1/alpha);\n",
"is = delta*ks;\n",
"cs = ys-is;\n",
"ab = 0;\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"function hernodes(nstate)\n",
"TOL = sqrt(2.2204e-16);\n",
"MAXIT = 30;\n",
"PIM4 = pi^(-1/4);\n",
"\n",
"n = nstate;\n",
"m = (n+1)/2;\n",
"m = convert(Int64, floor(m));\n",
"x = zeros(n,1);\n",
"w = zeros(n,1);\n",
"\n",
"z = 0;\n",
"z1 = 0;\n",
"pp = 0;\n",
"for i= 1:m\n",
"\n",
"# Initialize the first four roots\n",
"if i == 1; \n",
" z= sqrt(2*n+1) - 1.85575 * (2*n+1)^(-1/6);\n",
"elseif i== 2;\n",
" z= z - 1.14 * (n^.426)/z;\n",
"elseif i== 3;\n",
" z= 1.86 * z - .86 * x[1];\n",
"elseif i== 4;\n",
" z= 1.91 * z - .91 * x[2];\n",
"else;\n",
" z= 2 * z - x[i-2];\n",
"end;\n",
"\n",
"for its= 1:MAXIT\n",
" p1 = PIM4;\n",
" p2 = 0;\n",
" for j = 1:n\n",
" p3 = p2;\n",
" p2 = p1;\n",
" p1 = z*sqrt(2/j)*p2 - sqrt((j-1)/j)*p3;\n",
" end;\n",
"\n",
" pp = p2 * sqrt(2*n);\n",
" z1 = z;\n",
" z = z1 - p1/pp;\n",
"if abs(z-z1) < TOL; break; end;\n",
"end;\n",
"\n",
"x[i] = z;\n",
"x[n+1-i]= -z;\n",
"w[i]= 2/(pp*pp);\n",
"w[n+1-i]= w[i];\n",
"end; \n",
" \n",
"return x,w\n",
" end;"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"function transprob(ym,wm,m,r,s)\n",
"# this function computes the transition matrix\n",
"# Variables:\n",
"# ym is the vector of quadrature points\n",
"# wm is the vector of weights\n",
"# m is the mean of the process\n",
"# r is the rho of the process\n",
"# s is the conditional std.dev. of the process\n",
" \n",
"n,n0=size(ym) ; # get the number of quadrature points n\n",
"xx = ones(n,1); \n",
"xx=round.(Int64,vec(xx)); #make it integer: indices value in Julia must be integer\n",
"x=ym[:,xx] ; # get x, nxn matrix, whose ji element is consumption\n",
" # growth in state j - so x is the value of\n",
" # consumption growth at time t, if the state is\n",
" # j at time t and i at time t+1 (notice it's\n",
" # constant across i)\n",
"\n",
"y=x' ; # also an nxn matrix whose ji element is\n",
" # consumption growth in state i - so y is the\n",
" # value of consumption growth at time t+1, if\n",
" # the state is j at time t and i at time t+1\n",
" # (notice it's constant across j)\n",
"\n",
"w=wm[:,xx]' ; # converts the weight vector to a nxn matrix\n",
" # whose ji element is w(i) for all j\n",
" \n",
"f=(y-m*(1-r)-x*r) ;\n",
"c1=exp.(-(f.*f)./(2.0*s*s))./(sqrt(2.0*pi)*s) ;\n",
"\n",
"\n",
"f=(y-m*(1-r)-m*r) ;\n",
"c2=exp.(-(f.*f)./(2.0*s*s))./(sqrt(2.0*pi)*s) ;\n",
"\n",
"p=(c1.*w)./c2 ; # builds the transition matrix with elemets\n",
" # p(j,i)=f[y(i)|x(j)]w(i)/f[y(i)|mu]\n",
"\n",
"sm=sum(p',1)' ; # creates column vector with elements\n",
" # s(j)=sum(i) p(j,i)\n",
"\n",
"p=p./sm[:,xx];\n",
"\n",
"return p;\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Markov Chain technological process\n",
"rho = 0.8;\n",
"se = 0.2;\n",
"ma = 0;\n",
"agrid,wmat=hernodes(nstate);\n",
"agrid=agrid*sqrt(2)*se; \n",
"PI=transprob(agrid,wmat,0,rho,se)\n",
"at=agrid+ma;"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#define functions rcheb transfo itransfo\n",
"function rcheb(nn);\n",
"mod=nn-floor(nn/2)*2;\n",
"n1=floor(nn/2);\n",
"k=collect(1:n1).';\n",
"r1=cos.((2*k-1)*pi/(2*nn));\n",
"r1=[r1 -r1];\n",
"if mod==1;\n",
" r1=[r1 0];\n",
"end;\n",
"r1=vec(r1);\n",
"rr=real(sort!(r1));\n",
"return rr\n",
"end;\n",
"\n",
"function transfo(x,xmin,xmax);\n",
" z=(2*(x-xmin)/(xmax-xmin))-1; \n",
"return z\n",
"end;\n",
"\n",
"function itransfo(x,xmin,xmax);\n",
"z=0.5*(x+1)*(xmax-xmin)+xmin;\n",
"return z\n",
"end;\n",
"\n",
"function cheb(xx,nn);\n",
"cc=real(cos.(kron(acos.(complex(xx)),nn))) \n",
"return cc\n",
" end;"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# grid for the capital stock\n",
"kmin = log(1.2);\n",
"ksup = log(6);\n",
"rk = rcheb(nodes); #roots\n",
"kt = exp.(itransfo(rk,kmin,ksup)) #grid\n",
"vnbk = collect(0:nbk).';\n",
"XX = cheb(rk,vnbk);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"#initial conditions\n",
"a0=repmat([-0.2 0.65 0.04 0 0]',nstate,1);\n",
"a0=vec(a0);\n",
"param=[alpha beta delta sigma]';"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# function makepoly\n",
"\n",
"function makepoly(XA,XW);\n",
"nba = size(XA,2);\n",
"nba1 = size(XA,1);\n",
"nbw = size(XW,2);\n",
"nmax = max(nba,nbw);\n",
"XX = Float64[] ; \n",
" \n",
"for i=1:nbw\n",
" for j=1:nba \n",
" XX=[XX; kron(XW[:,i],XA[:,j])];\n",
" end\n",
"end\n",
"return XX\n",
" end;\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"f = function(theta);\n",
"RHS=Float64[];\n",
"LHS=Float64[];\n",
"resid=Float64[]; \n",
"ct = 0.0; \n",
"c1 = 0.0; \n",
"lt = length(theta);\n",
"lt = round.(Int64,lt/nstate);\n",
"theta = reshape(theta,lt,nstate);\n",
"\n",
"for i = 1:nstate\n",
" \n",
" ct = exp.(XX*theta[:,i]);\n",
" k1 = exp.(at[i])*kt.^alpha+(1-delta)*kt-ct;\n",
" rk1 = transfo(log.(k1),kmin,ksup);\n",
" vnbk = collect(0:nbk).'; \n",
" xk1 = cheb(complex(rk1),vnbk);\n",
" aux = 0;\n",
"\n",
" for j=1:nstate;\n",
" c1 = exp.(xk1*theta[:,j]);\n",
" aux = aux+PI[i,j]*beta*(alpha*exp.(at[j])*k1.^(alpha-1)+1-delta).*c1.^(-sigma);\n",
" end;\n",
" \n",
" RHS = [RHS;-sigma*log.(ct)];\n",
" LHS = [LHS;log.(aux)];\n",
"end;\n",
"resid = LHS-RHS;\n",
"return vec(resid);\n",
"end;\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Display final rule:[0.161955 0.0163026; 0.343494 0.383178; 0.00949691 0.0084339; 4.27686e-5 -7.84181e-5; -1.06076e-5 -1.22279e-6]\n"
]
}
],
"source": [
"#Pkg.add(\"NLsolve\");\n",
"using NLsolve;\n",
"sol = nlsolve(not_in_place(f),a0);\n",
"fsol = reshape(sol.zero,ncoef,nstate);\n",
"println(\"Display final rule:\",fsol)\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"#final\n",
"lt = length(sol.zero);\n",
"nb = 1000;\n",
"kt = collect(kmin:(ksup-kmin)/(nb-1):ksup);\n",
"rk = transfo(kt,kmin,ksup);\n",
"rk = vec(rk);\n",
"XX = cheb(rk,vnbk);\n",
"kt = exp.(kt);\n",
"ct=[];k1=[];ii=[];\n",
"for i=1:nstate\n",
" ct = cat(i,ct, exp.(XX*fsol[:,i]));\n",
" k1 = cat(i, k1, exp.(at[i])*kt.^alpha+(1-delta)*kt-ct[:,i]);\n",
" ii = cat(i,ii, exp.(at[i])*kt.^alpha-ct[:,i]);\n",
" end; "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plotly() \n",
"plot(kt,ct, linewidth=1,title=\"Consumption vs Capital Stock\", label=\"Consumption\")\n"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Julia 0.6.1",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}