Skip to content

Instantly share code, notes, and snippets.

@nilshg

nilshg/ijulia

Last active August 29, 2015 14:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save nilshg/0f4489c2a8cb99f3142f to your computer and use it in GitHub Desktop.
Save nilshg/0f4489c2a8cb99f3142f to your computer and use it in GitHub Desktop.
Download ZIP
Julia Income and Wealth
Raw
ijulia
{
"metadata": {
"language": "Julia",
"name": ""
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"using PyPlot\n",
"using DataFrames\n",
"using Distributions\n",
"using Optim"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stderr",
"text": [
"INFO: Loading help data...\n"
]
}
],
"prompt_number": 1
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Parameters"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"tol = 1e-9 # Tolerance Level for Iterations\n",
"T = 40 # Number of working life periods\n",
"r = 0.015 # Real interest rate\n",
"R = 1 + r # Handy for notation\n",
"dbeta = 0.966 # Discount factor (Note: 1/beta = 1.035)\n",
"gamma = 0.8 # Strenght of habits (gamma=0 --> no habits, gamma=1 --> only relative consumption matters)\n",
"lambda = 0.8 # Persistence of habits (lambda=0 --> no persistence)\n",
"sigma = 2 # Coefficient of relative risk aversion\n",
"agents = 100 # Number of agents sharing an alpha/beta combination\n",
"bs = 1000 # Different alpha/beta combinations in the population;\n",
"\n",
"mu_a = 2 # Mean of the income intercept\n",
"mu_b = 0.005 # Mean of income growth rates\n",
"var_a = 0.005 # Cross-sectional variance of alpha\n",
"var_b = 0.00037 # Cross-sectional variance of beta\n",
"var_eta = 0.029 # Variance of persistent shock\n",
"var_eps = 0.047 # Variance of transitory shock\n",
"ar1rho = 0.82 # Persistence of AR(1) component of income\n",
"FPU = 0.65 # Lambda;\n",
"\n",
"init_beta = 0.04\n",
"init_z = 0\n",
"init_varbeta = 0.01\n",
"init_varz = 0.01\n",
"\n",
"Sgridpoints = 15\n",
"wgridpoints = 12\n",
"hgridpoints = 11\n",
"ygridpoints = 8\n",
"\n",
"# Utility function\n",
"function u(c, h, gamma = 0.8, lambda = 0.8, sigma = 2)\n",
" (c*h^(-gamma))^(1-sigma)/(1-sigma)\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Income Distribution"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"guvenen_distribution = true\n",
"\n",
"if guvenen_distribution # Use distribution from Guvenen's paper\n",
" Yit = readcsv(\"C:/Users/tew207/Dropbox/QMUL/PhD/Code/Guvenen FORTRAN Code/LaborReal.csv\")\n",
" alfabeta = readcsv(\"C:/Users/tew207/Dropbox/QMUL/PhD/Code/Guvenen FORTRAN Code/alfabeta.csv\")\n",
" alpha = alfabeta[:, 1]\n",
" alpha = reshape(repmat(alpha, 1, agents)', agents*bs, 1)\n",
" beta = alfabeta[:, 2]\n",
" beta = reshape(repmat(beta, 1, agents)', agents*bs, 1)\n",
"else # Draw a new distribution\n",
" # Draw some alphas and betas\n",
" alpha = zeros(bs)\n",
" beta1 = zeros(bs)\n",
" for i = 1:bs\n",
" alpha[i] = mu_a + sqrt(var_a)*randn()\n",
" beta1[i] = mu_b + sqrt(var_b)*randn()\n",
" end\n",
"\n",
" # Sort \n",
" sort!(beta1);\n",
"\n",
" beta2 = zeros(length(beta1))\n",
" for i in [1:length(beta2)]\n",
" if i < length(beta2)/2\n",
" beta2[i] = beta1[i] + sqrt(var_b)*randn()\n",
" else \n",
" beta2[i] = beta1[i] + sqrt(var_b)*randn()+0.03\n",
" end\n",
" end\n",
" # Set break point at which growth rates switches from beta to beta2\n",
" BR = 30;\n",
" \n",
" # Create one beta matrix that holds each agents beta for each year\n",
" beta = zeros(agents*bs, T)\n",
"\n",
" for t = 1:T\n",
" for i = 1:bs\n",
" if t < BR\n",
" beta[(i-1)*agents+1:agents*i,t] = beta1[i]\n",
" else\n",
" beta[(i-1)*agents+1:agents*i,t] = beta2[i]\n",
" end\n",
" end\n",
" end\n",
"\n",
" # Reshape alpha into 100,000x1 vector with each of the 1,000 unique values repeated 100 times\n",
" alpha = reshape(repmat(alpha,1,100)', agents*bs, 1)\n",
" \n",
" # Draw the income distribution:\n",
" Yit = zeros(bs*agents, T)\n",
" z = zeros(bs*agents, T)\n",
" \n",
" for t = 1:T\n",
" for i = 1:bs*agents \n",
" Yit[i, t] = exp(alpha[i] + beta[i, t]*t + z[t] + sqrt(var_eps)*randn()) \n",
" if t < T\n",
" z[i, t+1] = rho*z[i, t] + sqrt(var_eta)*randn()\n",
" end\n",
" end\n",
" end\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stderr",
"text": [
"Warning: imported binding for beta overwritten in module Main\n"
]
}
],
"prompt_number": 3
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Belief Calculation"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"H = [ones(1,T); linspace(1,T,T)'; ones(1,T)] \n",
"F = [1 0 0; 0 1 0; 0 0 ar1rho] \n",
"Q = [0 0 0; 0 0 0; 0 0 var_eta] \n",
"S_f_i = zeros(3, agents*bs,T) \n",
"S_0_i = repmat([mean(alpha); mean(beta); init_z], 1, agents*bs) \n",
"P_f = zeros(3, 3, T) \n",
"P_0 = [0.005 -0.0002 0; -0.0002 0.0001 0; 0 0 0.0885] # Directly out of Guvenen's paper \n",
"\n",
"# Add some prior knowledge on beta\n",
"for i = 1:agents*bs\n",
" S_0_i[2, i] = 0.65*beta[i] + 0.35*S_0_i[2, i];\n",
"end\n",
" \n",
" \n",
"# Forecast from initial beliefs\n",
"S_f_i[:,:,1] = S_0_i; \n",
"P_f[:, :, 1] = P_0;\n",
"\n",
"# Evolution of Var-Cov-Matrix\n",
"for t = 1:T-1\n",
" Ht = H[:, t]\n",
" Pf = P_f[:, :, t]\n",
" P_f[:, :, t+1] = F*( Pf - Pf * Ht .* (Ht' * Pf * Ht + var_eps).^(-1.0) * Ht' * Pf )*F' + Q\n",
"end\n",
"\n",
"# Calculate Standard Deviation (needed later on)\n",
"stdy = zeros(1,T)\n",
"for t = 1:T\n",
" std = sqrt( H[:, t]' * P_f[:, 3*t-2:3*t] * H[:, t] )\n",
" stdy[t] = std[1]\n",
"end \n",
"\n",
"# Calculate Beliefs\n",
"tic()\n",
"for t = 1:T-1\n",
" Pc = P_f[:, :, t]\n",
" Hc = H[:, t]\n",
" for i = 1:agents*bs\n",
" Sfc = S_f_i[:, i, t]\n",
" S_f_i[:, i, t+1] = F*( Sfc + Pc * Hc .* (Hc' * Pc * Hc + R).^(-1.0) .* (log(Yit[i, t]) - Hc' * Sfc) )\n",
" end\n",
" if t % 5 == 0\n",
" @printf \"In Belief calculation, period is %d\\n\" t\n",
" end\n",
"end\n",
"toc()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"In Belief calculation, period is "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"5\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"In Belief calculation, period is 10\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"In Belief calculation, period is 15\n",
"In Belief calculation, period is "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"20\n",
"In Belief calculation, period is "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"25\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"In Belief calculation, period is 30\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"In Belief calculation, period is 35\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 28.206507302 seconds\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"28.206507302"
]
}
],
"prompt_number": 4
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Grid Construction"
]
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Belief Grid"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"maxalph = zeros(1, T);\n",
"minalph = zeros(1, T);\n",
"maxbeta = zeros(1, T);\n",
"minbeta = zeros(1, T);\n",
"maxzet = zeros(1, T);\n",
"minzet = zeros(1,T);\n",
"\n",
"for t = 1:T\n",
" maxalph[1, t] = maximum(S_f_i[1, :, t]) \n",
" maxbeta[1, t] = maximum(S_f_i[2, :, t]) \n",
" maxzet[1, t] = maximum(S_f_i[3, :, t]) \n",
" minalph[t] = minimum(S_f_i[1, :, t]) \n",
" minbeta[t] = minimum(S_f_i[2, :, t]) \n",
" minzet[t] = minimum(S_f_i[3, :, t]) \n",
"end;\n",
"\n",
"ALPHA = zeros(Sgridpoints-1, T)\n",
"BETA = zeros(Sgridpoints-1, T)\n",
"Z = zeros(Sgridpoints-1, T)\n",
"\n",
"for t = 1:T\n",
" ALPHA[:, t] = linspace(minalph[t], maxalph[t], Sgridpoints-1)';\n",
" BETA[:, t] = linspace(minbeta[t], maxbeta[t], Sgridpoints-1)';\n",
" Z[:, t] = linspace(minzet[t], maxzet[t], Sgridpoints-1)';\n",
"end\n",
"\n",
"# Matrix with minimum belief triplet in all columns\n",
"MINBELIEF = zeros(3, agents*bs, T)\n",
"for t = 1:T\n",
" MINBELIEF[1, :, t] = minalph[t]\n",
" MINBELIEF[2, :, t] = minbeta[t]\n",
" MINBELIEF[3, :, t] = minzet[t]\n",
"end\n",
"\n",
"# Matrix that measures how far the an agent's belief is away from the minimum\n",
"DIST = S_f_i - MINBELIEF\n",
"\n",
"# Matrix that measures how large one bin is\n",
"DIFF = zeros(3, agents*bs, T)\n",
"for t = 1:T\n",
" DIFF[1, :, t] = ALPHA[2, t] - ALPHA[1, t]\n",
" DIFF[2, :, t] = BETA[2, t] - BETA[1, t]\n",
" DIFF[3, :, t] = Z[2, t] - Z[1, t]\n",
"end\n",
"\n",
"BELIEFGRID = ceil(DIST./DIFF)\n",
"BELIEFGRID = max(BELIEFGRID, 1);\n",
"\n",
"Qt = zeros(T, 1) \n",
"S_f = zeros(3, 1200, 40)\n",
"for t = 1:T\n",
" a = unique(BELIEFGRID[:, :, t]', 1)\n",
" S_f[:, 1:size(a, 1), t] = a'\n",
" Qt[t] = size(a, 1) \n",
"end\n",
"\n",
"for t = 1:T\n",
" for i = 1:Qt[t]\n",
" S_f[1, i, t] = ALPHA[S_f[1, i, t], t]\n",
" S_f[2, i, t] = BETA[S_f[2, i, t], t]\n",
" S_f[3, i, t] = Z[S_f[3, i, t], t]\n",
" end\n",
"end\n",
"\n",
"Ybelief = zeros(ygridpoints, 1200, T)\n",
"for t = 1:T\n",
" for j = 1:Qt[t]\n",
" ymean = H[:, t]'*S_f[:, j, t]\n",
" y_low = exp(ymean - 3*stdy[t])\n",
" y_high = exp(ymean + 3*stdy[t])\n",
" Ybelief[:, j, t] = linspace(y_low[1], y_high[1], ygridpoints)\n",
" end\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Habit Grid"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"hgrid = zeros(hgridpoints, T)\n",
"for t = 1:T\n",
" hgrid[:, t] = linspace(0.1, 2*maximum(Ybelief[ygridpoints, :, t]), hgridpoints)\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "heading",
"level": 4,
"metadata": {},
"source": [
"Wealth Grid"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"power = 2\n",
"wmin = zeros(1, T)\n",
"wmax = zeros(1, T)\n",
"wgrid = zeros(wgridpoints, T)\n",
"wgridexp = zeros(wgridpoints, T)\n",
"\n",
"for t = 1:T\n",
" wmin[t] = -2*minimum(Ybelief[1, 1:Qt[t], t])\n",
" wmax[t] = 4*maximum(Ybelief[ygridpoints, :, t])\n",
" \n",
" wdistexp = (wmax[t] - wmin[t])^(1/power)\n",
" winc = wdistexp/(wgridpoints-1)\n",
" for i = 1: wgridpoints\n",
" wgridexp[i, t] = (i-1)*winc\n",
" end\n",
" wgrid[:, t] = wgridexp[:, t].^power + wmin[t]\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Compute Policy Functions"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"wprime = zeros(wgridpoints, hgridpoints, ygridpoints, 1200, T)\n",
"c = zeros(wgridpoints, hgridpoints, ygridpoints, 1200, T)\n",
"V = zeros(wgridpoints, hgridpoints, ygridpoints, 1200, T);\n",
"\n",
"for w = 1:wgridpoints\n",
" for h = 1:hgridpoints\n",
" for y = 1:ygridpoints\n",
" for q = 1:Qt[T]\n",
" if wgrid[w, T] < 0\n",
" c[w, h, y, q, T] = -(1/(5*wgrid[w, T]))\n",
" else\n",
" c[w, h, y, q, T] = wgrid[w, T]\n",
" end\n",
" ht = hgrid[h, T]\n",
" V[w, h, y, q, T] = u(c[w, h, y, q, T], ht)\n",
" end\n",
" end\n",
" end\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function regs(w, h, y, a, b, z, t)\n",
" eps01 = abs(w)+1\n",
" eps11 = 0.1\n",
" ylow = y/(exp(a+b*t+z)+0.4)\n",
" \n",
" X = [ 1; w; h; y; a; b; z; \n",
" w^2; w^3; h^2; h^3; y^2; y^3; a^2; b^2; b^3; z^2; z^3; \n",
" sqrt(w+eps01); sqrt(y+eps11); sqrt(w+eps01)*y; sqrt(w+eps01)*a; sqrt(w+eps01)*b; sqrt(w+eps01)*z; \n",
" sqrt(y+eps11)*w; sqrt(y+eps11)*a; sqrt(y+eps11)*b; sqrt(y+eps11)*z; sqrt(w+eps01)*sqrt(y+eps11); \n",
" sqrt(h); sqrt(w+eps01)*h; sqrt(y)*h; sqrt(w+eps01)*sqrt(h); sqrt(y)*sqrt(h); \n",
" w*y; w*a; w*b; w*z; (w^2)*y; (w^2)*a; (w^2)*b; (w^2)*z; w*(y^2); w*(a^2); w*(b^2); w*(z^2) ; \n",
" w*h; y*h; a*h; b*h; z*h; w^2*h; y^2*h; h*a^2; h*b^2; h*z^2; \n",
" y*a; y*b; y*z; (y^2)*a; (y^2)*b; (y^2)*z; y*(a^2); y*(b^2); y*(z^2); \n",
" a*b; a*z; (a^2)*b; (a^2)*z; a*(b^2); a*(z^2); b*z; (b^2)*z; b*(z^2) ; \n",
" w*y*a; w*y*b; w*y*z; w*a*b; w*a*z; w*b*z; y*a*b; y*a*z; y*b*z; a*b*z; \n",
" w*y*h; w*h*a; w*h*b; w*h*z; y*a*h; y*b*h; y*h*z; \n",
" w*y*a*b; w*y*a*z; w*y*b*z; w*a*b*z; y*a*b*z; \n",
" w*(a^3); w*(b^3); w*(z^3); (w^2)*(y^2); (w^2)*(a^2); (w^2)*(b^2); (w^2)*(z^2); y*(b^3); y*(z^3); \n",
" (y^2)*(a^2); (y^2)*(b^2); (y^2)*(z^2); w*(h^3); w^2*h^2; y*(h^3); y^2*h^2; \n",
" sqrt(w+eps01)*a*b; sqrt(w+eps01)*y*a; sqrt(w+eps01)*y*b; sqrt(w+eps01)*y*z; sqrt(w+eps01)*a*z; \n",
" sqrt(w+eps01)*b*z; sqrt(w+eps01)*h*y; sqrt(w+eps01)*h*b; \n",
" sqrt(y+eps11)*w*a; sqrt(y+eps11)*w*b; sqrt(y+eps11)*w*z; sqrt(y+eps11)*a*b; \n",
" sqrt(y+eps11)*a*z; sqrt(y+eps11)*b*z; sqrt(y+eps11)*w*a*z; sqrt(y+eps11)*w*b*z; sqrt(y+eps11)*a*b*z; \n",
" sqrt(w+eps01)*(a^2); sqrt(w+eps01)*(y^2); sqrt(w+eps01)*(b^2); sqrt(w+eps01)*(z^2); sqrt(w+eps01)*(a^2)*z; \n",
" sqrt(w+eps01)*(b^2)*z; sqrt(y+eps11)*(w^2); sqrt(y+eps11)*(b^2); sqrt(y+eps11)*(z^2); sqrt(y+eps11)*(a^2); \n",
" sqrt(y+eps11)*sqrt(w+eps01)*z; sqrt(y+eps11)*sqrt(w+eps01)*a; sqrt(y+eps11)*sqrt(w+eps01)*b ; \n",
" log(w+eps01)*(a^2); log(w+eps01)*(y^2); log(w+eps01)*(b^2); log(w+eps01)*(z^2); log(w+eps01)*(a^2)*z; \n",
" log(h); log(h)*w; log(h)*y; log(h)*a; log(h)*b; log(h)*z; \n",
" log(w+eps01)*(b^2)*z; log(y+eps11)*(w^2); log(y+eps11)*(b^2); log(y+eps11)*(z^2); log(y+eps11)*(a^2); \n",
" log(y+eps11)*log(w+eps01)*z; log(y+eps11)*log(w+eps01)*a; log(y+eps11)*log(w+eps01)*b ; \n",
" (w+eps01)^0.92; (y+eps11)^0.92; ((w+eps01)^0.92)*y; ((w+eps01)^0.92)*a; ((w+eps01)^0.92)*b; \n",
" ((w+eps01)^0.92)*z; ((y+eps11)^0.92)*w; \n",
" ((y+eps11)^0.92)*a; ((y+eps11)^0.92)*b; ((y+eps11)^0.92)*z; ((w+eps01)^0.92)*((y+eps11)^0.92); \n",
" ((w+eps01)^0.92)*a*b; ((w+eps01)^0.92)*y*a; \n",
" ((w+eps01)^0.92)*y*b; ((w+eps01)^0.92)*y*z; ((w+eps01)^0.92)*a*z; ((w+eps01)^0.92)*b*z; \n",
" ((y+eps11)^0.92)*w*a; ((y+eps11)^0.92)*w*b; ((y+eps11)^0.92)*w*z; \n",
" ((y+eps11)^0.92)*a*b; ((y+eps11)^0.92)*a*z; ((y+eps11)^0.92)*b*z; ((y+eps11)^0.92)*w*a*z; \n",
" ((y+eps11)^0.92)*w*b*z; ((y+eps11)^0.92)*a*b*z; \n",
" sqrt(w+eps01)*(1/ylow); log(w+eps01-0.9)*(1/ylow); sqrt(w+eps01)*(ylow); log(w+eps01-0.9)*(ylow); \n",
" sqrt(w+eps01)*(ylow^(-2)); \n",
" log(w+eps01-0.9)*(ylow^(-2)); sqrt(w+eps01)*(1/ylow)*a; log(w+eps01-0.9)*(1/ylow)*a; \n",
" sqrt(w+eps01)*(1/ylow)*b; log(w+eps01-0.9)*(1/ylow)*b; \n",
" sqrt(w+eps01)*(1/ylow)*z; log(w+eps01-0.9)*(1/ylow)*z ]\n",
" return X\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"nreg = length(regs(1,1,1,1,1,1,1))\n",
"coeffs = zeros(nreg, T)\n",
"maxrelerr = zeros(7, T);\n",
"\n",
"for t = T:-1:40\n",
" tic()\n",
" index = 1\n",
"\n",
" dim = int(ygridpoints*wgridpoints*hgridpoints*Qt[t])\n",
" X = zeros(dim, 6)\n",
" Y = zeros(dim, 1)\n",
"\n",
" for w = 1:wgridpoints\n",
" wt = wgrid[w, t]\n",
" for h = 1:hgridpoints\n",
" ht = hgrid[h, t]\n",
" for y = 1:ygridpoints\n",
" for q = 1:Qt[t]\n",
" yt = Ybelief[y, q, t]\n",
" a = S_f[1, q, t]\n",
" b = S_f[2, q, t]\n",
" z = S_f[3, q, t]\n",
" X[index, :] = [wt, ht, yt, a, b, z]\n",
" Y[index, :] = V[w, h, y, q, t]\n",
" index = index+1\n",
" end\n",
" end\n",
" end\n",
" end\n",
"\n",
" X2 = zeros(dim, nreg)\n",
" for i = 1:dim\n",
" (wt, ht, yt, at, bt, zt) = X[i, :]\n",
" X2[i, :] = regs(wt, ht, yt, at, bt, zt, t)\n",
" end\n",
"\n",
" coeffs[:, t] = (X2'*X2)^(-1)*(X2'*Y)\n",
"\n",
" # Goodness of fit\n",
" resid = Y - X2*coeffs[:, t]\n",
" maxval = maximum(resid)\n",
" maxloc = indmax(resid)\n",
" avgerr = mean(resid)\n",
" maxrelerr[1, t] = [X2[maxloc, :]*coeffs[:, t] / Y[maxloc]][1]\n",
" @printf \"Value Function interpolation took %.1f seconds\\n\" toc()\n",
"\n",
" # MAXIMIZATION\n",
"\n",
" coeffst = coeffs[:, t]\n",
" sigmay = stdy[t]\n",
" Ht = H[:, t]\n",
"\n",
" if t < T\n",
" wmin = wgrid[1, t+1]\n",
" else\n",
" wmin = -1\n",
" end\n",
"\n",
" for q = 1:5\n",
" tic()\n",
" muy = [Ht'*S_f[:, q, t]][1]\n",
" for y = 1:ygridpoints\n",
" yt = Ybelief[y, q, t]\n",
" for h = 1:hgridpoints\n",
" ht = hgrid[h, t]\n",
" for w = 1:wgridpoints\n",
" wt = wgrid[w, t] # current wealth\n",
" xt = wt + yt # Cash on Hand\n",
" (a, b, z) = S_f[:, q, t] # Belief about income process\n",
" \n",
" expy = LogNormal(muy, sigmay) # Distribution of expected income\n",
"\n",
" function EVprime(wp, y=yt, h=ht, a=a, b=b, z=z, t=t)\n",
" y_l = y - 3*stdy[t] # Lower bound for expected income\n",
" y_h = y + 3*stdy[t] # Higher bound for expected income\n",
"\n",
" # regs*coeffs gives the interpolated value of the value function\n",
" # pdf(exp, y) is the probability of getting this value\n",
" EVp(y) = regs(wp, y, h, a, b, z, t)'*coeffst*pdf(expy, y) \n",
"\n",
" quadgk(EVp, y_l, y_h)[1][1] # Integrate from low to high income expectation\n",
" end\n",
" \n",
" # Bellman: u(consumption, habits) + beta * E[V'(w')]\n",
" Blmn(wp) = -( u(xt - wp/R, ht) + dbeta * EVprime(wp) )\n",
" \n",
" # Maximize value function over choice of w'\n",
" Optimum = optimize(Blmn, wmin, xt - 0.01)\n",
"\n",
" # Store optimal choice of w' and associated value of the VF\n",
" wprime[w, h, y, q] = Optimum.minimum\n",
" V[w, h, y, q] = V_opt = -( Optimum.f_minimum )\n",
" end\n",
" end\n",
" end\n",
" @printf \"Q is %d out of %d, time elapsed %.2f seconds\\n\" q Qt[t] toc()\n",
" end\n",
"end;"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"32.2547639 seconds\n",
"Value Function interpolation took 32.3 seconds\n"
]
},
{
"output_type": "stream",
"stream": "stderr",
"text": [
"WARNING: Module Distributions not defined on process 3\n",
"fatal error on 3: ERROR: Distributions not defined\n",
" in deserialize at serialize.jl:376\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Worker 3 terminated."
]
},
{
"output_type": "stream",
"stream": "stderr",
"text": [
"WARNING: Module Distributions not defined on process 4\n",
"fatal error on 4: ERROR: Distributions not defined\n",
" in deserialize at serialize.jl:376\n"
]
},
{
"ename": "LoadError",
"evalue": "stream is closed or unusable\nwhile loading In[12], in expression starting on line 5",
"output_type": "pyerr",
"traceback": [
"stream is closed or unusable\nwhile loading In[12], in expression starting on line 5",
"",
" in check_open at stream.jl:294",
" in write at stream.jl:730",
" in takebuf_array at iobuffer.jl:216"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Worker 4 terminated."
]
},
{
"output_type": "stream",
"stream": "stderr",
"text": [
"WARNING: Module Distributions not defined on process 5\n",
"fatal error on 5: ERROR: Distributions not defined\n",
" in deserialize at serialize.jl:376\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"Worker 5 terminated.\n",
"\n"
]
}
],
"prompt_number": 12
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment