{
"cells": [
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Aiyagari Model adapted from George Hall code, http://people.brandeis.edu/~ghall/econ303/\n",
"# tested in Julia 0.6\n",
"# this code is part of chapter 6, \"Heterogeneous Agents Models\" 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 parameter values\n",
"sigma = 1.50; # risk aversion \n",
"beta = 0.98; # subjective discount factor \n",
"prob = [ .8 .2; .5 .5]; # prob(i,j) = probability (s(t+1)=sj | s(t) = si)\n",
"delta = 0.97; # 1 - depreciation\n",
"A = 1.00; # production technology\n",
"alpha = 0.25; # capital's share of income\n",
"theta = 0.05; # non-rental income if unemployed is theta*wage\n",
"Kstart = 10.0; # initial value for aggregate capital stock\n",
"g = 0.20; # relaxation parameter\n",
"\n",
"# form capital grid\n",
"maxkap = 20; # maximum value of capital grid \n",
"inckap = 0.025; # size of capital grid increments\n",
"nkap = trunc(Int,maxkap/inckap+1); # number of grid points: make it integer - Julia indexes must be integer\n",
"\n",
"#global variables\n",
"decis = zeros(nkap,2);\n",
"lambda = zeros(nkap,2);\n",
"probk = zeros(nkap,1);\n",
"\n",
"# calculate aggregate labor supply\n",
"D = [0.0 0.0; 0.0 0.0];\n",
"ed,ev = eig(prob);\n",
"\n",
"# make a matrix from the eigenvalues\n",
"edm = diagm(ed)\n",
"(emax,inmax) = findmax(edm)\n",
"\n",
"D[inmax,inmax] = emax;\n",
"pinf = ev*D*inv(ev);\n",
"pempl = pinf[inmax,inmax];\n",
"N = 1.0*pempl + theta*(1-pempl);\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ITERATING ON K\n",
"\n",
"Iter metric meanK Kold\n",
"1 0.5499633921667572 4.500366078332428 10.0\n",
"2 0.4544142973494472 4.855752699010765 8.900073215666486\n",
"3 0.3404741462790561 5.336361597447647 8.091209112335342\n",
"4 0.21417551289286063 5.9253049236885325 7.540239609357803\n",
"5 0.10094935807960213 6.488675647864649 7.21725267222395\n",
"6 0.033215148107074696 6.836655109672292 7.07153726735209\n",
"7 0.009285198741766663 6.959336392381948 7.024560835816131\n",
"8 0.0006063274285822773 7.015767221563982 7.011515947129295\n"
]
}
],
"source": [
"liter = 1;\n",
"maxiter = 50;\n",
"toler = 0.001;\n",
"metric = 10.0;\n",
"K = Kstart;\n",
"Kold= K;\n",
"wage=1.0;\n",
"rent=1.0;\n",
"println(\"ITERATING ON K\");\n",
"println(\"\"); \n",
"println(\"Iter metric meanK Kold\");\n",
"\n",
"\n",
"# loop to find fixed point for agregate capital stock\n",
"while (metric[1] > toler) & (liter <= maxiter);\n",
"\n",
"\n",
"# calculate rental rate of capital and wage\n",
" #\n",
" wage = (1-alpha) * A * K^(alpha) * N^(-alpha);\n",
" rent = (alpha) * A * K^(alpha-1) * N^(1-alpha);\n",
" # \n",
" # tabulate the utility function such that for zero or negative\n",
" # consumption utility remains a large negative number so that\n",
" # such values will never be chosen as utility maximizing \n",
" #\n",
" util1=-10000*ones(nkap,nkap); # utility when employed \n",
" util2=-10000*ones(nkap,nkap); # utility when unemployed \n",
" for i=1:nkap;\n",
" kap=(i-1)*inckap;\n",
" for j=1:nkap; \n",
" kapp = (j-1)*inckap;\n",
" cons1 = wage + (rent + delta)*kap - kapp; \n",
" if cons1 > .0;\n",
" util1[j,i]=(cons1)^(1-sigma)/(1-sigma);\n",
" end;\n",
" cons2 = theta*wage + (rent + delta)*kap - kapp;\n",
" if cons2 > .0;\n",
" util2[j,i]=(cons2)^(1-sigma)/(1-sigma);\n",
" end;\n",
" end;\n",
" end;\n",
" #\n",
" # initialize some variables\n",
" #\n",
" v = zeros(nkap,2);\n",
" tdecis1 = zeros(nkap,2);\n",
" tdecis2 = zeros(nkap,2);\n",
"\n",
" test = 10;\n",
" rs,cs = size(util1)\n",
" r1=zeros(cs,cs);\n",
" r2=zeros(cs,cs);\n",
"\n",
" #\n",
" # iterate on Bellman's equation and get the decision \n",
" # rules and the value function at the optimum\n",
"\n",
"\n",
" while test != 0;\n",
" for i=1:cs;\n",
" r1[:,i]=util1[:,i]+beta*(prob[1,1]*v[:,1]+ prob[1,2]*v[:,2]);\n",
" r2[:,i]=util2[:,i]+beta*(prob[2,1]*v[:,1]+ prob[2,2]*v[:,2]);\n",
" end;\n",
" (tv1,inds1)=findmax(r1,1);\n",
" tdecis1 =map(x->ind2sub(r1, x)[1], inds1) #to find the relative position of the max in a column\n",
" (tv2,inds2)=findmax(r2,1);\n",
" tdecis2 =map(x->ind2sub(r2, x)[1], inds2)\n",
"\n",
" tdecis=[tdecis1' tdecis2'];\n",
" tv=[tv1' tv2'];\n",
"\n",
" test=maximum((tdecis-decis));\n",
" copy!(v, tv);\n",
" copy!(decis, tdecis);\n",
" end; \n",
"\n",
" decis=(decis-1)*inckap;\n",
"\n",
"# form transition matrix\n",
"# trans is the transition matrix from state at t (row)\n",
"# to the state at t+1 (column) \n",
"# The eigenvector associated with the unit eigenvalue\n",
"# of trans' is the stationary distribution. \n",
" \n",
" g2=spzeros(cs,cs);\n",
" g1=spzeros(cs,cs);\n",
" for i=1:cs\n",
" g1[i,tdecis1[i]]=1;\n",
" g2[i,tdecis2[i]]=1;\n",
" end\n",
" trans=[ prob[1,1]*g1 prob[1,2]*g1; prob[2,1]*g2 prob[2,2]*g2];\n",
" trans=trans';\n",
" probst = (1/(2*nkap))*ones(2*nkap,1);\n",
" test=1;\n",
" while test > 10.0^(-8);\n",
" probst1 = trans*probst;\n",
" test = maximum(abs.(probst1-probst));\n",
" copy!(probst, probst1);\n",
" end;\n",
" \n",
"# vectorize the decision rule to be conformable with probst\n",
"# calculate new aggregate capital stock meanK \n",
"\n",
"kk=vec(decis);\n",
"meanK=probst'*kk;\n",
"\n",
"# calculate measure over (k,s) pairs\n",
"# lambda has same dimensions as decis\n",
"lambda=reshape(probst, cs,2)\n",
"\n",
"\n",
"# calculate stationary distribution of k\n",
"d1,v1=eig(prob');\n",
"d1m = diagm(d1)\n",
"dmax,imax=findmax(diag(d1m))\n",
"\n",
"probst1=v1[:,imax];\n",
"ss=sum(probst1);\n",
"probst1=probst1./ss;\n",
"probk=sum(lambda',1)'\n",
"# form metric and update K\n",
"Kold= K;\n",
"Knew= g*meanK[1] + (1-g)*Kold;\n",
"metric = abs.((Kold-meanK)./Kold);\n",
"K = Knew;\n",
"println(liter,\" \", metric[1],\" \",meanK[1],\" \",Kold);\n",
"liter = liter+1;\n",
"end;"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PARAMETER VALUES\n",
"\n",
"sigma beta delta A alpha theta\n",
"1.5 0.98 0.97 1.0 0.25 0.05\n",
"\n",
"EQUILIBRIUM RESULTS \n",
"\n",
"K N wage rent\n",
"7.0115 0.7286 1.321 0.0458\n"
]
}
],
"source": [
"#print results\n",
"println(\"PARAMETER VALUES\");\n",
"println(\"\");\n",
"println(\"sigma beta delta A alpha theta\"); \n",
"println(sigma,\" \",beta, \" \",delta,\" \", A, \" \",alpha,\" \",theta);\n",
"println(\"\"); \n",
"println(\"EQUILIBRIUM RESULTS \");\n",
"println();\n",
"println(\"K N wage rent\");\n",
"println(round(Kold,4),\" \", round(N,4), \" \",round(wage,4),\" \", round(rent,4) );"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"markov (generic function with 1 method)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# function markov(T,n,s0,V);\n",
"# chain generates a simulation from a Markov chain of dimension\n",
"# the size of T\n",
"# T is transition matrix\n",
"# n is number of periods to simulate\n",
"# s0 is initial state\n",
"# V is the quantity corresponding to each state\n",
"# state is a matrix recording the number of the realized state at time t\n",
"\n",
"function markov(T,n,s0,V);\n",
" \n",
"r,c = size(T);\n",
"v1,v2 = size(V);\n",
"\n",
"#rand('uniform');\n",
"#using Distributions\n",
"X=rand(Uniform(0,1), n-1,1)\n",
"\n",
"state=zeros(2,99);\n",
"chain=[];\n",
"s=zeros(r,1);\n",
"s[s0]=1\n",
"cum=T*triu(ones(size(T)));\n",
"ppi =[];\n",
"\n",
"state[:,1]=s\n",
"\n",
"for k=1:length(X);\n",
"#k=1\n",
" state[:,k]=s;\n",
" ppi =[0 s'*cum];\n",
" ss1= convert(Array{Float64},((X[k].<=ppi[2:r+1])))\n",
" ss2= convert(Array{Float64},(X[k].>ppi[1:r]))\n",
" s=(ss1.*ss2) \n",
" s=reshape(s,2,1)\n",
"end \n",
"\n",
"chain=V*state\n",
"return state,chain \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SIMULATING LIFE HISTORY\n"
]
}
],
"source": [
"# simulate life histories of the agent\n",
"#rand('uniform');\n",
"using Distributions\n",
"\n",
"println(\"SIMULATING LIFE HISTORY\");\n",
"k = Kold; # initial level of capital \n",
"n = 100; # number of periods to simulate\n",
"s0 = 1; # initial state \n",
"hist = zeros(n-1,2);\n",
"cons = zeros(n-1,1);\n",
"invest = zeros(n-1,1);\n",
"\n",
"grid = collect(0:inckap:maxkap);\n",
"r,c = size(prob);\n",
"T = prob;\n",
"V = collect(1:r)';\n",
"\n",
"state, chain= markov(prob,n,s0,V)\n",
"\n",
"chain=convert(Array{Int64,2},chain);\n",
"state=convert(Array{Int64,2},state);"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"\n",
"for i = 1:n-1;\n",
" hist[i,:] = [ k chain[i] ];\n",
" I1 = trunc(Int,k/inckap) ;\n",
" I2 = trunc(Int,k/inckap) + 1;\n",
" if I1 == 0;\n",
" I1=1;\n",
" println(\"N.B. I1 = 0\");\n",
" end;\n",
" if I2 > nkap;\n",
" I2 = nkap;\n",
" println(\"N.B. I2 > nkap\");\n",
" end;\n",
" weight = (grid[I2,1] - k)/inckap; \n",
" kprime = weight*(decis[I1,chain[i]]) + (1-weight)*(decis[I2,chain[i]]);\n",
" if chain[i] == 1;\n",
" cons[i] = wage + (rent + delta)*k - kprime;\n",
" elseif chain[i] == 2;\n",
" cons[i] = wage*theta + (rent + delta)*k - kprime;\n",
" else;\n",
" println(\"something is wrong with chain\");\n",
" chain\n",
" end;\n",
" k = kprime;\n",
" invest[i] = kprime;\n",
"end;\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#plots\n",
"#Income distribution\n",
"income = [ (rent*grid + wage) (rent*grid + wage*theta) ] ; \n",
"#sort income\n",
"pinc = sortrows(hcat(vec(income), 1:length(vec(income))), by = x -> x[1])[:,1]\n",
"index= sortrows(hcat(vec(income), 1:length(vec(income))), by = x -> x[1])[:,2]\n",
"index=trunc.(Int,index)\n",
"plambda = vec(lambda);"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"plotly() \n",
"plot(pinc,plambda[index], linewidth=1,title=\"Income Distribution\", ylabel=\"% of agents\", label=\"income level\")"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" \n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#plots\n",
"#Capital distribution\n",
"#using Plots\n",
"plotly() \n",
"plot(grid,probk, linewidth=1,title=\"Capital Distribution\", ylabel=\"% of agents\", label=\"capital\")\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.6.0",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}