{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# PEA 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 with Julia\"\n",
"# Academic Press - Elsevier\n",
"# for comments, email at: petre(dot)caraiani(at)gmail(dot)com\n",
"\n",
"# set algorithm parameters\n",
"long = 200;\n",
"init = 50;\n",
"slong = init+long;\n",
"T = init+1:slong-1;\n",
"T1 = init+2:slong;\n",
"tol = 1e-6;\n",
"crit = 1;\n",
"\n",
"#set model parameters\n",
"gam = 1;\n",
"sigma = 1;\n",
"delta = 0.1;\n",
"beta = 0.95;\n",
"alpha = 0.3;\n",
"ab = 0;\n",
"rho = 0.9;\n",
"se = 0.01;\n",
"\n",
"#steady state\n",
"ksy =(alpha*beta)/(1-beta*(1-delta));\n",
"yss = ksy^(alpha/(1-alpha));\n",
"kss = yss^(1/alpha);\n",
"iss = delta*kss;\n",
"css = yss-iss;\n",
"csy = css/yss;\n",
"lss = css^(-sigma);"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"#generate own shocks\n",
"#using(DataFrames)\n",
"srand(1)\n",
"e = se*randn(slong,1);\n",
"#ee = DataFrame(e)\n",
"#writetable(\"output.csv\", ee);\n",
"#load shocks from Collard(2015) generated in Matlab to find the same solution\n",
"e=readdlm(\"e1_noise.csv\");"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"a = zeros(slong,1);\n",
"a[1]= ab+e[1];\n",
"for i=2:slong;\n",
"a[i]=rho*a[i-1]+(1-rho)*ab+e[i];\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#initial solution\n",
"ncont = 3; # of static equations\n",
"nbend = 1; # endogenous predetermined variables\n",
"nshoc = 1; # of shocks\n",
"nback = nbend+nshoc; # of state variables\n",
"nforw = 1; # of costate variables\n",
"nstat = nback+nforw; # of state and costate variables\n",
"Mcc = zeros(ncont,ncont);\n",
"Mcs = zeros(ncont,nstat);\n",
"Mss0 = zeros(nstat,nstat);\n",
"Mss1 = zeros(nstat,nstat);\n",
"Msc0 = zeros(nstat,ncont);\n",
"Msc1 = zeros(nstat,ncont);\n",
"Mse = zeros(nstat,nshoc);"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#setup the matrices\n",
"# Output\n",
"Mcc[1,1] = 1;\n",
"Mcs[1,1] = alpha;\n",
"Mcs[1,2] = 1;\n",
"\n",
"# investment\n",
"Mcc[2,1] = 1;\n",
"Mcc[2,2] = -iss/yss;\n",
"Mcc[2,3] = -css/yss;\n",
"\n",
"# consumption\n",
"Mcc[3,3] = -sigma;\n",
"Mcs[3,3] = 1;\n",
"\n",
"# capital\n",
"Mss0[1,1] = 1;\n",
"Mss1[1,1] = delta-1;\n",
"Msc1[1,2] = delta;\n",
"\n",
"# technology shock\n",
"Mss0[2,2] = 1;\n",
"Mss1[2,2] = -rho;\n",
"Mse[2,1] = 1;\n",
"\n",
"# Euler\n",
"Mss0[3,1] = (1-beta*(1-delta));\n",
"Mss0[3,3] = -1;\n",
"Mss1[3,3] = 1;\n",
"Msc0[3,1] = (1-beta*(1-delta));"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Solving the system \n",
"M0=inv(Mss0-Msc0*inv(Mcc)*Mcs);\n",
"M1=(Mss1-Msc1*inv(Mcc)*Mcs);\n",
"W=-M0*M1;\n",
"\n",
"# MU -> eigenvalues, P -> eigenvectors\n",
"v,w = eig(W)\n",
"r=[abs.(v) w']\n",
"vsize=size(v,1)\n",
"for i = 1:vsize\n",
" for j = i+1:vsize\n",
" if real(r[i,1])> real(r[j,1])\n",
" tmp = r[i,:];\n",
" r[i,:]=r[j,:];\n",
" r[j,:]=tmp;\n",
" elseif real(r[i,1])== real(r[j,1])\n",
" if imag(r[i,1])> imag(r[j,1])\n",
" tmp = r[i,:];\n",
" r[i,:]=r[j,:];\n",
" r[j,:]=tmp; \n",
" end; \n",
" end;\n",
" end; \n",
"end;\n",
"\n",
"lam= r[:,1];\n",
"P=r[:,2:4]';\n",
"\n",
"Q=inv(P);\n",
"lam=diagm(lam);"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"#Direct solution\n",
"Gamma=-inv(Q[nback+1:nstat,nback+1:nstat])*Q[nback+1:nstat,1:nback];\n",
"MSS=W[1:nback,1:nback]+W[1:nback,nback+1:nstat]*Gamma;\n",
"PI=inv(Mcc)*(Mcs[:,1:nback]+Mcs[:,nback+1:nstat]*Gamma);\n",
"MSE=[zeros(nbend,nshoc);eye(nshoc)];\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"S = zeros(nback,long+init);\n",
"S[:,1]= MSE*e[1];\n",
"for i = 2:long+init;\n",
" S[:,i]= MSS*S[:,i-1]+MSE*e[i];\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"lb = Gamma*S;\n",
"lb = lss*exp.(lb);\n",
"lbv= vec(lb);\n",
"k = log(kss)+S[1,:];\n",
"kv = vec(k);\n",
"ek = exp.(kv);\n",
"a = S[2,:];\n",
"av = vec(a);\n",
"ea = exp.(av);\n",
"T = init+1:init+long-1;\n",
"T1 = init+2:init+long;"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"XX = [ones(long-1,1) kv[T] av[T] kv[T].*kv[T] av[T].*av[T] kv[T].*av[T]];\n",
"yy = log.(beta*lb[T1].*(alpha*ea[T1].*ek[T1].^(alpha-1)+1-delta));\n",
"b0 = \\(XX,yy);"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solution[1.63964, -3.00701, 1.60307, 1.27846, 1.46655, -2.19256]\n",
"solution[1.76819, -3.27105, 1.53844, 1.41379, 1.47546, -2.12414]\n",
"solution[1.83915, -3.41625, 1.44735, 1.48787, 1.49376, -2.0291]\n",
"solution[1.88338, -3.50636, 1.37241, 1.53361, 1.51939, -1.9512]\n",
"solution[1.91351, -3.5675, 1.31829, 1.56453, 1.54463, -1.89512]\n",
"solution[1.93547, -3.61199, 1.28147, 1.58698, 1.56597, -1.85714]\n",
"solution[1.95228, -3.646, 1.2573, 1.60412, 1.58271, -1.83232]\n",
"solution[1.96553, -3.67281, 1.24174, 1.61765, 1.59534, -1.81647]\n",
"solution[1.97611, -3.69424, 1.2318, 1.62847, 1.60468, -1.80643]\n",
"solution[1.98458, -3.71142, 1.2254, 1.63715, 1.61152, -1.80003]\n",
"solution[1.99135, -3.72515, 1.22119, 1.6441, 1.61652, -1.79587]\n",
"solution[1.99673, -3.73607, 1.21832, 1.64963, 1.62019, -1.79307]\n",
"solution[2.00098, -3.74469, 1.2163, 1.654, 1.62291, -1.79111]\n",
"solution[2.00431, -3.75145, 1.21483, 1.65743, 1.62494, -1.78968]\n",
"solution[2.00691, -3.75674, 1.21371, 1.66011, 1.62646, -1.78861]\n",
"solution[2.00894, -3.76086, 1.21286, 1.66219, 1.62762, -1.78779]\n",
"solution[2.01052, -3.76406, 1.2122, 1.66381, 1.6285, -1.78715]\n",
"solution[2.01174, -3.76654, 1.21169, 1.66507, 1.62918, -1.78666]\n",
"solution[2.01269, -3.76847, 1.2113, 1.66605, 1.6297, -1.78628]\n",
"solution[2.01343, -3.76996, 1.21099, 1.6668, 1.63011, -1.78598]\n",
"solution[2.01401, -3.77113, 1.21076, 1.66739, 1.63042, -1.78575]\n",
"solution[2.01445, -3.77204, 1.21058, 1.66785, 1.63067, -1.78558]\n",
"solution[2.0148, -3.77274, 1.21044, 1.66821, 1.63086, -1.78545]\n",
"solution[2.01508, -3.7733, 1.21034, 1.66849, 1.631, -1.78535]\n",
"solution[2.01529, -3.77373, 1.21026, 1.66871, 1.63112, -1.78527]\n",
"solution[2.01546, -3.77407, 1.2102, 1.66888, 1.63121, -1.78522]\n",
"solution[2.01559, -3.77434, 1.21015, 1.66901, 1.63128, -1.78517]\n",
"solution[2.0157, -3.77455, 1.21012, 1.66912, 1.63133, -1.78514]\n",
"solution[2.01578, -3.77471, 1.21009, 1.6692, 1.63137, -1.78512]\n",
"solution[2.01584, -3.77484, 1.21007, 1.66927, 1.6314, -1.7851]\n",
"solution[2.01589, -3.77494, 1.21006, 1.66932, 1.63143, -1.78508]\n",
"solution[2.01593, -3.77502, 1.21005, 1.66936, 1.63145, -1.78507]\n",
"solution[2.01596, -3.77508, 1.21004, 1.66939, 1.63147, -1.78506]\n",
"solution[2.01599, -3.77513, 1.21003, 1.66942, 1.63148, -1.78506]\n",
"solution[2.01601, -3.77517, 1.21003, 1.66944, 1.63149, -1.78505]\n",
"solution[2.01602, -3.77521, 1.21002, 1.66945, 1.6315, -1.78505]\n",
"solution[2.01603, -3.77523, 1.21002, 1.66946, 1.6315, -1.78505]\n",
"solution[2.01604, -3.77525, 1.21002, 1.66947, 1.63151, -1.78504]\n",
"solution[2.01605, -3.77526, 1.21002, 1.66948, 1.63151, -1.78504]\n",
"solution[2.01606, -3.77528, 1.21001, 1.66949, 1.63151, -1.78504]\n",
"solution[2.01606, -3.77529, 1.21001, 1.66949, 1.63152, -1.78504]\n",
"solution[2.01607, -3.77529, 1.21001, 1.6695, 1.63152, -1.78504]\n",
"solution[2.01607, -3.7753, 1.21001, 1.6695, 1.63152, -1.78504]\n",
"solution[2.01607, -3.77531, 1.21001, 1.6695, 1.63152, -1.78504]\n",
"solution[2.01607, -3.77531, 1.21001, 1.6695, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77531, 1.21001, 1.66951, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77531, 1.21001, 1.66951, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77532, 1.21001, 1.66951, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77532, 1.21001, 1.66951, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77532, 1.21001, 1.66951, 1.63152, -1.78504]\n",
"solution[2.01608, -3.77532, 1.21001, 1.66951, 1.63152, -1.78504]\n"
]
}
],
"source": [
"#Main Loop\n",
"iter=1;\n",
"crit=1;\n",
"while crit>tol\n",
"k = zeros(slong+1,1);\n",
"lb = zeros(slong,1);\n",
"X = zeros(slong,length(b0));\n",
"k[1]= kss;\n",
"for i= 1:slong;\n",
" X[i,:]= [1 log.(k[i]) a[i] log(k[i])*log(k[i]) a[i]*a[i] log(k[i])*a[i]];\n",
" lb[i] = exp.(X[i,:]'*b0);\n",
" k[i+1] = exp.(a[i])*k[i]^alpha+(1-delta)*k[i]-lb[i]^(-1/sigma);\n",
"end\n",
"y = beta*lb[T1].*(alpha*exp.(a[T1]).*k[T1].^(alpha-1)+1-delta);\n",
"bt = X[T,:]\\log.(y);\n",
"b = copy(gam*bt+(1-gam)*b0);\n",
"crit = maximum(abs.(b-b0));\n",
"b0 = copy(b);\n",
"println(\"solution\",b0)\n",
"iter=iter+1;\n",
"end;\n"
]
},
{
"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": [
"#Pkg.add(\"Plots\")\n",
"#Pkg.add(\"PlotlyJS\")\n",
"#using PlotlyJS\n",
"using Plots\n",
"plotly() # Choose the Plotly.jl backend for web interactivity\n",
"#plot(k',c',linewidth=1,label=\"Consumption\",title=\"Consumption vs capital stock\")\n",
"scatter(k[T],k[T1],label=\"Capital Stock\",title=\"k(t+1) vs. k(t)\")\n",
"\n",
"#Plots.plotlyjs() # specify backend\n",
"\n",
"#scatter!(k[T],k[T1], color=:blue, label=\"k\")"
]
}
],
"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
}