{ "cells": [ { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Huggett 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", "theta = 0.05; # non-interest income if unemployed\n", "wage = 1.00; # non-interest income if employed\n", "Rstart = 1.021; # initial gross interest rate\n", "F = -2.0; # borrowing constraint parameter \n", "g = 0.60; # relaxation parameter\n", "\n", "#initialize variables\n", "A = 1.0;\n", "Aold = 1.0;\n", "Anew = 1.0;\n", "meanA = 1.0;\n", "\n", "#form asset grid\n", "maxast = 8; # maximum value of asset grid \n", "minast = -5; # minimum value of asset grid\n", "incast = 0.5; # size of asset grid increments\n", "nasset = trunc(Int,((maxast-minast)/incast+1)); # number of grid points\n", "assetp = 1.0;\n", "\n", "#global variables\n", "decis = zeros(nasset,2);\n", "tdecis = zeros(nasset,2);\n", "lambda = zeros(nasset,2);\n", "\n", "# loop to find R such that sum(lambda*A) = 0 \n", "liter = 1;\n", "maxiter = 50;\n", "toler = 0.0001;\n", "step = 0.05;\n", "R = Rstart;\n", "flag = 1;\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ITERATING ON R\n", "\n", "Iter R A newstep\n", "1 1.021 7.381503742789281 -0.05\n", "2 0.9709999999999999 -0.42868278522192527 0.025\n", "3 0.9959999999999999 0.4769046747885715 -0.0125\n", "4 0.9834999999999999 0.01692921100922419 -0.0125\n", "5 0.971 -0.42868278522192527 0.00625\n", "6 0.97725 -0.4286827763313216 0.00625\n", "7 0.9834999999999999 0.01692921100922419 -0.003125\n", "8 0.9803749999999999 0.01692921967959421 -0.003125\n", "9 0.9772499999999998 -0.4286827763313216 0.0015625\n", "10 0.9788124999999999 -0.4286827763313216 0.0015625\n", "11 0.9803749999999999 0.01692921967959421 -0.00078125\n", "12 0.9795937499999999 0.01692921967959421 -0.00078125\n", "13 0.9788125 -0.4286827763313216 0.000390625\n", "14 0.979203125 0.016929215870976413 -0.0001953125\n", "15 0.9790078124999999 0.016929259815237283 -0.0001953125\n", "16 0.9788124999999999 -0.4286827763313216 9.765625e-5\n" ] } ], "source": [ "println(\"ITERATING ON R\");\n", "println(\"\");\n", "println(\"Iter R A newstep\");\n", "while (flag != 0) && (liter <= maxiter);\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(nasset,nasset); # utility when employed \n", " util2=-10000*ones(nasset,nasset); # utility when unemployed \n", "\n", "for i=1:nasset\n", " asset=(i-1)*incast + minast;\n", " for j=1:nasset\n", " assetp = (j-1)*incast + minast;\n", " cons = wage + R*asset - assetp;\n", " if assetp >= F && cons > 0;\n", " util1[j,i]=(cons)^(1-sigma)/(1-sigma);\n", " end;\n", " end\n", " for j=1:nasset\n", " assetp = (j-1)*incast + minast;\n", " cons = theta*wage + R*asset - assetp;\n", " if assetp>= F && cons > 0;\n", " util2[j,i]=(cons)^(1-sigma)/(1-sigma);\n", " end;\n", " end;\n", " end;\n", "\n", "# initialize some variables\n", "\n", " v = zeros(nasset,2);\n", " tdecis1 = zeros(nasset,2);\n", " tdecis2 = zeros(nasset,2); \n", "\n", " \n", " test1 = 10;\n", " test2 = 10;\n", " rs,cs = size(util1);\n", " r1=zeros(nasset,nasset);\n", " r2=zeros(nasset,nasset);\n", "\n", "# iterate on Bellman's equation and get the decision \n", "# rules and the value function at the optimum \n", " \n", " while (test1 != 0) || (test2 > .1);\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", "\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", " test1=maximum((tdecis-decis));\n", " test2=maximum(abs.(tv-v));\n", " copy!(v, tv);\n", " copy!(decis, tdecis); \n", " end;\n", " decis=(decis-1)*incast + minast;\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", " \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*nasset))*ones(2*nasset,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 asset meanA \n", "aa=vec(decis);\n", "meanA=(probst'*aa)[1];\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", "lambda=reshape(probst, cs,2)\n", "probk=sum(lambda',1); # stationary distribution of capital - sum by each column\n", "probk=probk'\n", "\n", " \n", "if liter == 1;\n", " A=copy(meanA);; \n", "if meanA > 0.0;\n", " step=copy(-step);\n", "end;\n", "end;\n", "\n", "Aold = copy(A);\n", "Anew = copy(meanA);\n", "\n", "if sign(Aold) != sign(Anew)\n", " step = copy(-.5*step);\n", "end;\n", "println(liter,\" \",R,\" \",meanA,\" \",step);\n", "if abs.(step) >= toler;\n", " R=copy(R+step);\n", "else;\n", " flag = 0;\n", "end;\n", " A=copy(Anew);\n", " liter = liter+1;\n", " \n", "end;" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "#display solution\n", "# calculate consumption and expected utility\n", "\n", "grid = collect(minast:incast:maxast) ; \n", "congood = wage*(ones(nasset,1)) + R*grid - grid[tdecis[:,1]];\n", "conbad = theta*wage*(ones(nasset,1)) + R*grid - grid[tdecis[:,2]];\n", "consum = [congood conbad ];\n", "cons2 = [congood.^2 conbad.^2];\n", "meancon = sum(diag(lambda'*consum));\n", "meancon2 = sum(diag(lambda'*cons2));\n", "varcon = ( meancon2 - meancon^2 );\n", "UTILITY = ((complex(consum)).^(1-sigma))./(1-sigma);\n", "UCEU2 = sum(diag(lambda'*UTILITY));\n" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PARAMETER VALUES\n", "\n", "sigma beta F theta\n", "1.5 0.98 -2.0 0.05\n", "\n", "EQUILIBRIUM RESULTS\n", "\n", "R A UCEU meancon varcon\n", "0.9788 -0.4287 -2.5003 0.7377 0.0634\n" ] } ], "source": [ "#print out results\n", "println(\"PARAMETER VALUES\");\n", "println(\"\");\n", "println(\"sigma beta F theta\"); \n", "println(sigma,\" \",beta,\" \",F,\" \", theta);\n", "println(\"\"); \n", "println(\"EQUILIBRIUM RESULTS\");\n", "println(\"\");\n", "println(\"R A UCEU meancon varcon\");\n", "println(round(R,4),\" \",round(meanA,4),\" \",real(round(UCEU2,4)),\" \",round(meancon,4),\" \",round(varcon,4));" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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 }