{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Replicating Maliar, Maliar, and Valli (2010, JEDC) to solve Krusell and Smith (1998, JPE) model using Julia\n", "\n", "By [Shunsuke Hori](https://github.com/Shunsuke-Hori)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Overview of the notebook\n", "This notebook solves the model of Krusell and Smith (1998, JPE) and succesfully replicating the result of Maliar, Maliar, and Valli (2010, JEDC)\n", "\n", "The solution strategy is as follows\n", "\n", "1. Solve individual problem by Euler equation method with 2D interpolation\n", " - Agents are boundedly rational. In the code, they take into account the information about the mean of capital\n", " - Aggregate law of motion is approximated by log-linear relation, i.e. $\\log(K_{t+1})=B1+B2\\log(K_{t})$ for good aggregate state and $\\log(K_{t+1})=B3+B4\\log(K_{t})$ for bad aggregate state \n", " - If specified, Howard's policy iteration is used\n", "1. Compute the path of aggregate capital using the policy function obtained by VFI\n", " - Simulation is based on Monte Carlo. That is, aggregate technology shocks and idiosyncratic employment shocks are drawn for many agents and many periods. Then, using the LLN, the aggregate capital is computed by aggregating all agents for all period.\n", "1. Update the coefficient of aggregate capital law of motion, $B1$, $B2$, $B3$ and $B4$, by regression\n", "1. Check convergence of $B1$, $B2$, $B3$ and $B4$\n", "\n", "Additionally, this notebook also includes the [result](http://localhost:8888/notebooks/KrusellSmith_Euler.ipynb#Implementation-with-value-function-iteration) with value function iteration as a solution method for individual utility maximization problem.\n", "\n", "NOTE: Regarding interpolation, Krusell and Smith uses various interpolation scheme depending on the purpose, including polynomial interpolation. Maliar, Maliar, and Valli uses spline interpolation in their paper. This notebook only uses linear interpolation because I could not find interpolation package for polynomial and spline 2D-interpolation.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Code to solve models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First thing to do is import some packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "using Interpolations # to use interpolation \n", "using QuantEcon # to use a function, gridmake\n", "using Plots # to plot the result\n", "pyplot()\n", "#plotlyjs()\n", "using Optim # to use minimization routine to maximize RHS of bellman equation\n", "using GLM # to regress\n", "using JLD # to save the result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model Setup\n", "\n", "Functions in this cell are prepared for model parameters and initial guess of the solutions \n", "\n", "- Types\n", " - TransitionMatrix: collection of transition matrix\n", " - KSParameter: collection of model parameters, functional forms, and grids\n", " - KSSolution: collection of solution, which is guess at first\n", "- Methods\n", " - create_transition_matrix: construct an instance of type TransitionMatrix\n", " - KSParameter: construct an instance of type KSParameter \n", " - place_polynominal_grid: create polynominal grid \n", " - r: compute interest rate ( = marginal productivity of capital) \n", " - w: compute wage rate ( = marginal productivity of labor)\n", " - KSSolution: construct an instance of type KSSolution\n", " - Guess can be the load of previous result" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/plain": [ "KSSolution" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Collection of transition matrix\n", "\"\"\"\n", "immutable TransitionMatrix\n", " P::Array{Float64,2} # 4x4 \n", " Pz::Array{Float64,2} # 2x2 aggregate shock\n", " Peps_gg::Array{Float64,2} # 2x2 idiosyncratic shock conditional on good to good\n", " Peps_bb::Array{Float64,2} # 2x2 idiosyncratic shock conditional on bad to bad\n", " Peps_gb::Array{Float64,2} # 2x2 idiosyncratic shock conditional on good to bad\n", " Peps_bg::Array{Float64,2} # 2x2 idiosyncratic shock conditional on bad to good\n", "end\n", "\n", "\"\"\"\n", "Collection of model parameters\n", "\"\"\"\n", "immutable KSParameter\n", " u::Function\n", " beta::Float64\n", " alpha::Float64\n", " delta::Float64\n", " theta::Float64\n", " l_bar::Float64\n", " k_grid::Vector{Float64}\n", " K_grid::Vector{Float64}\n", " z_grid::Vector{Float64}\n", " eps_grid::Vector{Float64}\n", " s_grid::Array{Float64,2}\n", " k_size::Int64\n", " K_size::Int64\n", " z_size::Int64\n", " eps_size::Int64\n", " s_size::Int64\n", " ug::Float64\n", " ub::Float64\n", " TransMat::TransitionMatrix # bunch of transition matrix\n", " mu::Float64\n", "end\n", "\n", "\"\"\" \n", "Create transition matrices for aggregate shock,\n", "idiosyncratic shock, and shock state\n", "\n", "##### Arguments\n", "- ug : unemployment rate in good state\n", "- ub : unemployment rate in bad state\n", "- zg_ave_dur : average duration of good state\n", "- zb_ave_dur : average duration of bad state\n", "- ug_ave_dur : average duration of unemployment in good state\n", "- ub_ave_dur : average duration of unemployment in bad state\n", "- puu_rel_gb2bb : prob. of u to u cond. on g to b relative to that of b to b \n", "- puu_rel_bg2gg : prob. of u to u cond. on b to g relative to that of g to g\n", "\"\"\"\n", "function create_transition_matrix(ug::Float64,ub::Float64,\n", " zg_ave_dur::Float64,zb_ave_dur::Float64,\n", " ug_ave_dur::Float64,ub_ave_dur::Float64,\n", " puu_rel_gb2bb::Float64,puu_rel_bg2gg::Float64)\n", " \n", " # probability of remaining in good state\n", " pgg = 1-1/zg_ave_dur\n", " # probability of remaining in bad state\n", " pbb = 1-1/zb_ave_dur\n", " # probability of changing from g to b\n", " pgb = 1-pgg\n", " # probability of changing from b to g\n", " pbg = 1-pbb \n", " \n", " # prob. of 0 to 0 cond. on g to g\n", " p00_gg = 1-1/ug_ave_dur\n", " # prob. of 0 to 0 cond. on b to b\n", " p00_bb = 1-1/ub_ave_dur\n", " # prob. of 0 to 1 cond. on g to g\n", " p01_gg = 1-p00_gg\n", " # prob. of 0 to 1 cond. on b to b\n", " p01_bb = 1-p00_bb\n", " \n", " # prob. of 0 to 0 cond. on g to b\n", " p00_gb=puu_rel_gb2bb*p00_bb\n", " # prob. of 0 to 0 cond. on b to g\n", " p00_bg=puu_rel_bg2gg*p00_gg\n", " # prob. of 0 to 1 cond. on g to b\n", " p01_gb=1-p00_gb\n", " # prob. of 0 to 1 cond. on b to g\n", " p01_bg=1-p00_bg\n", " \n", " # prob. of 1 to 0 cond. on g to g\n", " p10_gg=(ug - ug*p00_gg)/(1-ug)\n", " # prob. of 1 to 0 cond. on b to b\n", " p10_bb=(ub - ub*p00_bb)/(1-ub)\n", " # prob. of 1 to 0 cond. on g to b\n", " p10_gb=(ub - ug*p00_gb)/(1-ug)\n", " # prob. of 1 to 0 cond on b to g\n", " p10_bg=(ug - ub*p00_bg)/(1-ub)\n", " # prob. of 1 to 1 cond. on g to g\n", " p11_gg= 1-p10_gg\n", " # prob. of 1 to 1 cond. on b to b\n", " p11_bb= 1-p10_bb\n", " # prob. of 1 to 1 cond. on g to b\n", " p11_gb= 1-p10_gb\n", " # prob. of 1 to 1 cond on b to g\n", " p11_bg= 1-p10_bg\n", " \n", " # (g1) (b1) (g0) (b0)\n", " P=[pgg*p11_gg pgb*p11_gb pgg*p10_gg pgb*p10_gb;\n", " pbg*p11_bg pbb*p11_bb pbg*p10_bg pbb*p10_bb;\n", " pgg*p01_gg pgb*p01_gb pgg*p00_gg pgb*p00_gb;\n", " pbg*p01_bg pbb*p01_bb pbg*p00_bg pbb*p00_bb\n", " ]\n", " Pz=[pgg pgb;\n", " pbg pbb]\n", " Peps_gg=[p11_gg p10_gg\n", " p01_gg p00_gg]\n", " Peps_bb=[p11_bb p10_bb\n", " p01_bb p00_bb]\n", " Peps_gb=[p11_gb p10_gb\n", " p01_gb p00_gb]\n", " Peps_bg=[p11_bg p10_bg\n", " p01_bg p00_bg]\n", " TransMat=TransitionMatrix(P,Pz,Peps_gg,Peps_bb,Peps_gb,Peps_bg)\n", " return TransMat\n", "end\n", "\n", "\"\"\"\n", "Creates KSParameter instance\n", "\"\"\"\n", "function KSParameter(;\n", " beta::Float64=0.99,\n", " alpha::Float64=0.36,\n", " delta::Float64=0.025,\n", " theta::Float64=1.0,\n", " k_min::Float64=1e-16,\n", " k_max::Float64=1000.0,\n", " k_size::Int64=100,\n", " K_min::Float64=30.0,\n", " K_max::Float64=50.0,\n", " K_size::Int64=4,\n", " z_min::Float64=0.99,\n", " z_max::Float64=1.01,\n", " z_size::Int64=2,\n", " eps_min::Float64=0.0,\n", " eps_max::Float64=1.0,\n", " eps_size::Int64=2,\n", " ug::Float64=0.04,\n", " ub::Float64=0.1,\n", " zg_ave_dur::Float64=8.0,\n", " zb_ave_dur::Float64=8.0,\n", " ug_ave_dur::Float64=1.5,\n", " ub_ave_dur::Float64=2.5,\n", " puu_rel_gb2bb::Float64=1.25,\n", " puu_rel_bg2gg::Float64=0.75,\n", " mu::Float64=0.0\n", " )\n", " if theta == 1.0\n", " u = (c) -> log(c)\n", " else\n", " u = (c) -> (c^(1.0-theta)-1.0)/(1.0-theta)\n", " end\n", " l_bar=1/(1-ub)\n", " k_grid=place_polynominal_grid(k_min,k_max,k_size,degree=7.0) # individual capital grid\n", " K_grid=collect(linspace(K_min,K_max,K_size)) # aggregate capital grid\n", " z_grid=collect(linspace(z_max,z_min,z_size)) # aggregate technology shock\n", " eps_grid=collect(linspace(eps_max,eps_min,eps_size)) # idiosyncratic employment shock\n", " s_grid=gridmake(z_grid,eps_grid) # shock grid\n", " # collection of transition matrices\n", " TransMat=create_transition_matrix(ug,ub,\n", " zg_ave_dur,zb_ave_dur,\n", " ug_ave_dur,ub_ave_dur,\n", " puu_rel_gb2bb,puu_rel_bg2gg)\n", "\n", " ksp=KSParameter(u,beta,alpha,delta,theta,l_bar,k_grid,K_grid,z_grid,eps_grid,s_grid,\n", " k_size,K_size,z_size,eps_size,z_size*eps_size,ug,ub,TransMat,mu)\n", "\n", " return ksp\n", "end\n", "\n", "function place_polynominal_grid(k_min::Float64,k_max::Float64,k_size::Int64;degree::Float64=0.7)\n", " grid=Array{Float64}(k_size)\n", " for (i, x) in enumerate(linspace(0,0.5,k_size))\n", " grid[i]=(x/0.5)^degree*(k_max-k_min)+k_min\n", " end\n", " return grid\n", "end\n", "\"\"\"\n", "Compute interest rate given aggregate capital, labor, and productivity\n", "\n", "##### Arguments\n", "- alpha : capital share\n", "- z : aggregate shock\n", "- K : aggregate capital\n", "- L : aggregate labor\n", "\"\"\"\n", "r(alpha::Float64,z::Float64,K::Float64,L::Float64)=alpha*z*K^(alpha-1)*L^(1-alpha)\n", "\n", "\"\"\"\n", "Compute wage given aggregate capital, labor, and productivity\n", "\n", "##### Arguments\n", "- alpha : capital share\n", "- z : aggregate shock\n", "- K : aggregate capital\n", "- L : aggregate labor\n", "\"\"\"\n", "w(alpha::Float64,z::Float64,K::Float64,L::Float64)=(1-alpha)*z*K^(alpha)*L^(-alpha) \n", "\n", "\"\"\"\n", "Collection of KS solution\n", "\"\"\"\n", "type KSSolution\n", " k_opt::Array{Float64,3}\n", " value::Array{Float64,3}\n", " B::Vector{Float64}\n", " R2::Vector{Float64}\n", "end\n", "\n", "\"\"\"\n", "Create KSSolution instance\n", "\"\"\"\n", "function KSSolution(\n", " ksp::KSParameter;\n", " load_value::Bool=false,\n", " load_B::Bool=false,\n", " filename::String=\"result.jld\")\n", " if load_value || load_B\n", " result=load(filename)\n", " kss_temp=result[\"kss\"]\n", " end\n", " if load_value\n", " k_opt=kss_temp.k_opt\n", " value=kss_temp.value\n", " else\n", " k_opt=ksp.beta*repeat(ksp.k_grid,outer=[1,ksp.K_size,ksp.s_size])\n", " k_opt=0.9*repeat(ksp.k_grid,outer=[1,ksp.K_size,ksp.s_size])\n", " value=ksp.u.(0.1/0.9*k_opt)/(1-ksp.beta)\n", " end\n", " if load_B\n", " B=kss_temp.B\n", " else\n", " B=[0.0, 1.0, 0.0, 1.0]\n", " end\n", " kss=KSSolution(k_opt,value,B,[0.0,0.0])\n", " return kss\n", "end\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Shock Generation\n", "The functions in this cell are used to draw aggregate and idiosyncratic shocks\n", "- Methods\n", " - get_shock: translate an index shock combination grid to shock value\n", " - generate_shocks: draw aggregate and idiosyncratic shocks, which is the main function in this cell\n", " - draw_eps_shock (local function inside generate_shocks): draw idiosyncratic shocks given aggregate shocks" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "generate_shocks" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "\"\"\" \n", "Translate shock index into shock value\n", "\n", "##### Arguments\n", "- s_grid : shock grid\n", "- s_i : shock index\n", "\"\"\"\n", "get_shock(s_grid::Array{Float64,2},s_i::Integer) = \n", " s_grid[s_i,1], s_grid[s_i,2] \n", "\n", "\"\"\"\n", "Generate aggregate and idiosyncratic shock\n", "\n", "##### Arguments\n", "- ksp : instance of KSParameter type\n", "- z_shock_size : size of aggregate shock\n", "- population : size idiosyncratic shock in one period\n", "\"\"\"\n", "function generate_shocks(ksp::KSParameter;\n", " z_shock_size::Int64=11000,population::Int64=10000)\n", " \n", " \n", " \"\"\" \n", " Draw idiosyncratic shock given previous idiosyncratic shock and \n", " transition matrix.\n", " The transition matrix must be consistent with aggregate shock\n", "\n", " ##### Arguments\n", " - eps_shock : preallocated vector. current period shock is stored in it\n", " - eps_shock_before : previous period idiosyncratic shock\n", " - Peps : transition matrix of the current period\n", " \"\"\"\n", " function draw_eps_shock(eps_shock_before::Vector{Float64},\n", " Peps::Array{Float64,2})\n", " \n", " eps_shocks = similar(eps_shock_before)\n", " # loop over entire population\n", " for i=1:length(eps_shocks)\n", " rand_draw=rand()\n", " eps_shocks[i]=ifelse(eps_shock_before[i]==1.0,\n", " Float64(Peps[1,1]>rand_draw), # if employed before\n", " Float64(Peps[2,1]>rand_draw)) # if unemployed before\n", " end\n", " return eps_shocks\n", " end\n", " \n", " # unpack parameters\n", " Peps_gg=ksp.TransMat.Peps_gg\n", " Peps_bg=ksp.TransMat.Peps_bg\n", " Peps_gb=ksp.TransMat.Peps_gb\n", " Peps_bb=ksp.TransMat.Peps_bb\n", " \n", " # draw aggregate shock\n", " z_shock=simulate(MarkovChain(ksp.TransMat.Pz,ksp.z_grid),z_shock_size)\n", " \n", " ### Let's draw individual shock ###\n", " eps_shock=Array{Float64}(z_shock_size,population) # preallocation\n", " \n", " # first period\n", " rand_draw=rand(population)\n", " if z_shock[1]==ksp.z_grid[1] # if good\n", " eps_shock[1,:].=Int64.(rand_draw.>ksp.ug) # if draw is higher, become employed\n", " elseif z_shock[1]==ksp.z_grid[2] # if bad\n", " eps_shock[1,:].=Int64.(rand_draw.>ksp.ub) # if draw is higher, become employed\n", " else\n", " error(\"the value of z_shock[1] (=$(z_shock[1])) is strange\")\n", " end\n", " \n", " # from second period ...\n", " for t=2:z_shock_size\n", " if z_shock[t]==ksp.z_grid[1] && z_shock[t-1]==ksp.z_grid[1] # if g to g\n", " eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_gg)\n", " elseif z_shock[t]==ksp.z_grid[1] && z_shock[t-1]==ksp.z_grid[2] # if b to g\n", " eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_bg)\n", " elseif z_shock[t]==ksp.z_grid[2] && z_shock[t-1]==ksp.z_grid[1] # if g to b\n", " eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_gb)\n", " elseif z_shock[t]==ksp.z_grid[2] && z_shock[t-1]==ksp.z_grid[2] # if b to b\n", " eps_shock[t,:]=draw_eps_shock(eps_shock[t-1,:],Peps_bb)\n", " else\n", " error(\"the value of z_shock[t] (=$(z_shock[t])) is strange\")\n", " end\n", " end\n", " \n", " # adjustment\n", " for t=1:z_shock_size\n", " n_e=Int64(sum(eps_shock[t,:]))\n", " er_ideal = ifelse(z_shock[t] == ksp.z_grid[1],\n", " 1.0-ksp.ug, 1.0-ksp.ub)\n", " gap = Int64(round(er_ideal*population)) - n_e\n", " if gap > 0\n", " for j=1:gap\n", " become_employed_i = \n", " sample(find(x-> isapprox(x,ksp.eps_grid[2]), eps_shock[t,:]))\n", " eps_shock[t,become_employed_i] = ksp.eps_grid[1]\n", " end\n", " elseif gap < 0\n", " for j=1:(-gap)\n", " become_unemployed_i = \n", " sample(find(x-> isapprox(x,ksp.eps_grid[1]), eps_shock[t,:]))\n", " eps_shock[t,become_unemployed_i] = ksp.eps_grid[2]\n", " end\n", " end \n", " end\n", " \n", " return z_shock, eps_shock \n", "end\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Functions for Euler equation method\n", "\n", "Functions in this cell is used to solve individual problem by Euler equation method\n", "\n", "- Methods\n", " - compute_Kp_L: compute approximated $K'$ and $L$, which depend on current aggregate shock\n", " - compute_expectation_FOC: compute expectation term in Euler equation\n", " - euler_method!: find optimal policy by Euler equation methods, namely, iteration of Euler equation\n", " \n", " $$\\left(c\\right)^{-\\theta}=\\beta E\\left[\\left(c'\\right)^{-\\theta}(1-\\delta+r')\\right]$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "euler_method!" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Compute next period aggregate capital and labor\n", "\n", "##### Arguments\n", "- K : Current aggregate capital\n", "- s_i : current shock index\n", "- B : coefficient of ALM for capital\n", "\"\"\"\n", "function compute_Kp_L(K::AbstractFloat,s_i::Integer,B::Vector{Float64},ksp::KSParameter)\n", " Kp,L=ifelse(ksp.s_grid[s_i,1]==ksp.z_grid[1],\n", " (exp(B[1]+B[2]*log(K)),ksp.l_bar*(1-ksp.ug)),\n", " (exp(B[3]+B[4]*log(K)),ksp.l_bar*(1-ksp.ub)))\n", " Kp = ifelse(Kpksp.K_grid[end],ksp.K_grid[end],Kp)\n", " return Kp, L\n", "end\n", "\n", "\"\"\"\n", "Compute expectation term in Euler equation\n", "\n", "##### Arguments\n", "- kp : next period individual capital holding\n", "- Kp : next period aggregate capital\n", "- s_i : current state of shock\n", "- ksp : KSParameter instance\n", "\"\"\"\n", "function compute_expectation_FOC(kp::Float64,\n", " Kp::Float64,\n", " s_i::Int64,\n", " ksp::KSParameter)\n", " alpha, theta, delta, l_bar, mu, P =\n", " ksp.alpha, ksp.theta, ksp.delta, ksp.l_bar, ksp.mu, ksp.TransMat.P\n", " expec = 0.0\n", " for s_n_i = 1:ksp.s_size\n", " zp, epsp = ksp.s_grid[s_n_i,1], ksp.s_grid[s_n_i,2]\n", " Kpp, Lp = compute_Kp_L(Kp,s_n_i,kss.B,ksp)\n", " rn=r(alpha,zp,Kp,Lp)\n", " kpp=interpolate((ksp.k_grid,ksp.K_grid),\n", " kss.k_opt[:,:,s_n_i],Gridded(Linear()))\n", " cp = (rn+1-delta)*kp+\n", " w(alpha,zp,Kp,Lp)*(epsp*l_bar+mu*(1.0-epsp))-kpp[kp,Kp]\n", " expec = expec + P[s_i,s_n_i]*(cp)^(-theta)*(1-delta+rn)\n", " end \n", " return expec\n", "end\n", "\n", "\"\"\"\n", "Solve individual problem by Euler equation method\n", "\n", "##### Arguments\n", "- ksp : KSParameter instance\n", "- kss : KSSolution instance\n", "- n_iter : maximum number of iteration of Euler equation method\n", "- tol : tolerance of policy function convergence\n", "- update_k : degree of update of policy function\n", "\"\"\"\n", "function euler_method!(ksp::KSParameter,\n", " kss::KSSolution;\n", " n_iter::Integer=30000,\n", " tol::AbstractFloat=1e-8,\n", " update_k::AbstractFloat=1e-8)\n", " println(\"solving individual problem by Euler equation method\")\n", " alpha, beta, delta, theta, l_bar, mu = \n", " ksp.alpha, ksp.beta, ksp.delta, ksp.theta, ksp.l_bar, ksp.mu\n", " k_grid, k_size = ksp.k_grid, ksp.k_size\n", " K_grid, K_size = ksp.K_grid, ksp.K_size\n", " s_grid, s_size = ksp.s_grid, ksp.s_size\n", " k_min, k_max = minimum(k_grid), maximum(k_grid)\n", " counter=0\n", " k_opt_n=similar(kss.k_opt)\n", " while true\n", " counter += 1\n", " dif_k=0.0\n", " for s_i = 1:s_size\n", " z, eps = s_grid[s_i,1], s_grid[s_i,2]\n", " for K_i = 1:K_size\n", " K = K_grid[K_i]\n", " Kp, L = compute_Kp_L(K,s_i,kss.B,ksp)\n", " for k_i = 1:k_size\n", " k=k_grid[k_i]\n", " wealth = (r(alpha,z,K,L)+1-delta)*k+\n", " w(alpha,z,K,L)*(eps*l_bar+mu*(1.0-eps))\n", " expec=compute_expectation_FOC(kss.k_opt[k_i,K_i,s_i],Kp,s_i,ksp)\n", " cn=(beta*expec)^(-1.0/theta)\n", " k_opt_n[k_i,K_i,s_i] = wealth-cn\n", " k_opt_n[k_i,K_i,s_i] = ifelse(k_opt_n[k_i,K_i,s_i]>k_max,k_max,k_opt_n[k_i,K_i,s_i])\n", " k_opt_n[k_i,K_i,s_i] = ifelse(k_opt_n[k_i,K_i,s_i]=n_iter\n", " println(\"Euler method failed to converge with $counter iterations (dif =$dif_k)\")\n", " break\n", " end\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate path of aggregate capital\n", "- Mathods\n", " - simulate_aggregate_path!: simulating aggregate path of capital using shocks drawn" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/plain": [ "simulate_aggregate_path!" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Simulate aggregate capital's path using policy function \n", "and generated aggregate and idiosyncratic shock\n", "\n", "##### Arguments\n", "- ksp : KSParameter instance\n", "- z_shocks : aggregate shocks\n", "- eps_shocks : idiosyncratic shocks\n", "- k_population : initial capital holding of all agents\n", "\"\"\"\n", "function simulate_aggregate_path!(ksp::KSParameter,kss::KSSolution,\n", " z_shocks::Vector{Float64},eps_shocks::Array{Float64,2},\n", " k_population::Vector{Float64},K_ts::Vector{Float64})\n", " \n", " println(\"simulating aggregate path ... please wait ... \")\n", " \n", " T=length(z_shocks) # simulated duration\n", " N=size(eps_shocks,2) # number of agents\n", " \n", " \n", " # loop over T periods\n", " for (t,z) = enumerate(z_shocks)\n", " K_ts[t]=mean(k_population) # current aggrgate capital\n", " \n", " # s_i_base takes 1 when good and 2 when bad \n", " s_i_base=ifelse(z==ksp.z_grid[1],1,2) \n", " \n", " # loop over individuals\n", " for (i,k) in enumerate(k_population)\n", " eps = eps_shocks[t,i] # idiosyncratic shock\n", " s_i=s_i_base+2*(1-Int64(eps)) # transform (z,eps) to s_i\n", " # obtain next capital holding by interpolation\n", " \n", " itp_pol=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,s_i],Gridded(Linear()))\n", " k_population[i]=itp_pol[k,K_ts[t]]\n", " end\n", " end\n", " \n", " return nothing\n", "end\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Obtaining aggregate law of motion coefficient\n", "\n", "The functions in this cell are used to obtain the coefficient of approximate aggregate capital law of motion (ALM)\n", "\n", "- Methods\n", " - regress_ALM: regress aggregate capital law of motion and obtain approximate ALM\n", " - find_ALM_coef!: main iteration. update ALM coefficient until convergence" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "find_ALM_coef! (generic function with 1 method)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Obtain new aggregate law of motion coefficients \n", "using the aggregate capital flaw\n", "\n", "##### Arguments\n", "- ksp : KSParameter instance\n", "- z_shock : aggregate shocks\n", "- K_ts : aggregate capital flaw\n", "- n_discard : number of discarded samples\n", "\"\"\"\n", "function regress_ALM!(ksp::KSParameter,kss::KSSolution,\n", " z_shock::Vector{Float64},K_ts::Vector{Float64};\n", " n_discard::Int64=100)\n", " z_grid=ksp.z_grid\n", " n_g=count(i->(i==z_grid[1]),z_shocks[n_discard+1:end-1])\n", " n_b=count(i->(i==z_grid[2]),z_shocks[n_discard+1:end-1])\n", " B_n=Vector{Float64}(4)\n", " x_g=Vector{Float64}(n_g)\n", " y_g=Vector{Float64}(n_g)\n", " x_b=Vector{Float64}(n_b)\n", " y_b=Vector{Float64}(n_b)\n", " i_g=0\n", " i_b=0\n", " for t = n_discard+1:length(z_shocks)-1\n", " if z_shocks[t]==z_grid[1]\n", " i_g=i_g+1\n", " x_g[i_g]=log(K_ts[t])\n", " y_g[i_g]=log(K_ts[t+1])\n", " else\n", " i_b=i_b+1\n", " x_b[i_b]=log(K_ts[t])\n", " y_b[i_b]=log(K_ts[t+1])\n", " end\n", " end\n", " resg=lm(hcat(ones(n_g,1),x_g),y_g)\n", " resb=lm(hcat(ones(n_b,1),x_b),y_b)\n", " kss.R2=[r2(resg),r2(resb)]\n", " B_n[1],B_n[2]=coef(resg)\n", " B_n[3],B_n[4]=coef(resb)\n", " dif_B=maximum(abs(B_n-kss.B))\n", " println(\"difference of ALM coefficient is $dif_B and B =$B_n\")\n", " return B_n, dif_B\n", "end\n", "\n", "function find_ALM_coef!(\n", " ksp::KSParameter,\n", " kss::KSSolution,\n", " z_shocks::Vector{Float64},\n", " eps_shocks::Array{Float64,2};\n", " tol_ump::AbstractFloat=1e-8,\n", " max_iter_ump::Integer=100,\n", " Howard_on::Bool=true,\n", " Howard_n_iter::Integer=20,\n", " tol_B::AbstractFloat=1e-8,\n", " max_iter_B::Integer=20,\n", " update_B::AbstractFloat=0.3,\n", " T_discard::Integer=100,\n", " print_skip_VFI::Integer=10,\n", " method::Symbol=:Euler,\n", " update_k::AbstractFloat=0.3)\n", "\n", " K_ts=similar(z_shocks)\n", " # populate initial capital holdings\n", " k_population=37.9893*ones(size(eps_shocks,2))\n", " counter_B=0\n", " while true\n", " counter_B=counter_B+1\n", " println(\" --- Iteration over ALM coefficient: $counter_B ---\")\n", "\n", " # solve individual problem\n", " if method == :VFI\n", " solve_bellman!(ksp,kss,\n", " tol=tol_ump,\n", " max_iter=max_iter_ump,\n", " Howard=Howard_on,\n", " Howard_n_iter=Howard_n_iter,\n", " print_skip=print_skip_VFI)\n", " elseif method == :Euler\n", " euler_method!(ksp,kss,\n", " n_iter=max_iter_ump,\n", " tol=tol_ump,\n", " update_k=update_k)\n", " end\n", " \n", "\n", " # compute aggregate path of capital\n", " simulate_aggregate_path!(ksp,kss,z_shocks,eps_shocks,k_population,K_ts)\n", " \n", " # obtain new ALM coefficient by regression\n", " B_n,dif_B = regress_ALM!(ksp,kss,z_shocks,K_ts,n_discard=T_discard)\n", "\n", " # check convergence\n", " if dif_B < tol_B\n", " println(\"-----------------------------------------------------\")\n", " println(\"ALM coefficient successfully converged : dif =$dif_B\")\n", " println(\"-----------------------------------------------------\")\n", " kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B\n", " break\n", " elseif counter_B == max_iter_B\n", " println(\"----------------------------------------------------------------\")\n", " println(\"Iteration over ALM coefficient reached its maximum ($max_iter_B)\")\n", " println(\"----------------------------------------------------------------\")\n", " kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B\n", " break\n", " end\n", " \n", " # Update B\n", " kss.B .= update_B .* B_n .+ (1-update_B) .* kss.B\n", " end\n", " return K_ts\n", "end\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting\n", "The function is used to plot the path of true path of aggregate capital and approximated one\n", "- Methods\n", " - plot_ALM: plot the path of true path of aggregate capital and approximated one" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "plot_ALM" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Plot true and approximated ALM of capital \n", "\n", "##### Arguments\n", "- ksp.z_grid : aggregate shock grid\n", "- z_shocks : aggregate shock\n", "- kss.B : ALM coefficient\n", "- K_ts : actual path of capital\n", "\"\"\"\n", "function plot_ALM(z_grid::Vector{Float64},\n", " z_shocks::Vector{Float64},\n", " B::Vector{Float64},\n", " K_ts::Vector{Float64};\n", " T_discard=1000)\n", "\n", " K_ts_approx = similar(K_ts) # preallocation\n", "\n", " # compute approximate ALM for capital\n", " K_ts_approx[T_discard]=K_ts[T_discard]\n", "\n", " for t=T_discard:length(z_shocks)-1\n", " if z_shocks[t] == z_grid[1]\n", " K_ts_approx[t+1]=exp(B[1]+B[2]*log(K_ts_approx[t]))\n", " elseif z_shocks[t] == z_grid[2]\n", " K_ts_approx[t+1]=exp(B[3]+B[4]*log(K_ts_approx[t]))\n", " end\n", " end\n", "\n", " plot(T_discard+1:length(K_ts),K_ts[T_discard+1:end],lab=\"true\",color=:red,line=:solid)\n", " plot!(T_discard+1:length(K_ts),K_ts_approx[T_discard+1:end],lab=\"approximation\",color=:blue,line=:dash)\n", " title!(\"aggregate law of motion for capital\")\n", " \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation by Euler method" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, construct model instance ksp and initial guess of the solution inside kss \n", "\n", "(Grid size inconsistency is also checked, which may return error when exiting result is loaded by load_value=true)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# instance of KSParameter\n", "ksp=KSParameter(k_size=100,\n", " k_min=1e-16)\n", "# instance of KSSolution\n", "kss=KSSolution(ksp,load_value=false,load_B=false)\n", "if size(kss.k_opt,1) != length(ksp.k_grid)\n", " error(\"loaded data is inconsistent with k_size\")\n", "end\n", "if size(kss.k_opt,2) != length(ksp.K_grid)\n", " error(\"loaded data is inconsistent with K_size\")\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's draw the shock for stochastic simulation of aggregate law of motion" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# generate shocks\n", "srand(123)\n", "z_shocks,eps_shocks =generate_shocks(ksp;\n", " z_shock_size=11000,population=5000);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, the following cell solves the model with Euler equation method" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " --- Iteration over ALM coefficient: 1 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1868 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.4556643772779668 and B = [0.455664,0.877296,0.438162,0.880418]\n", " --- Iteration over ALM coefficient: 2 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 2005 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.01194706494540887 and B = [0.137682,0.963462,0.119502,0.96685]\n", " --- Iteration over ALM coefficient: 3 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1479 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.012429438071340365 and B = [0.149423,0.960332,0.135204,0.962694]\n", " --- Iteration over ALM coefficient: 4 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1516 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.007084614314802812 and B = [0.147807,0.960707,0.134031,0.962974]\n", " --- Iteration over ALM coefficient: 5 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1159 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.004062551435912554 and B = [0.146911,0.960937,0.133124,0.963222]\n", " --- Iteration over ALM coefficient: 6 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1044 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.002383876077763747 and B = [0.146451,0.961055,0.132628,0.963361]\n", " --- Iteration over ALM coefficient: 7 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1065 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0014303729909590501 and B = [0.146213,0.961116,0.132363,0.963436]\n", " --- Iteration over ALM coefficient: 8 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1050 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0008768443163971185 and B = [0.146088,0.961148,0.132222,0.963476]\n", " --- Iteration over ALM coefficient: 9 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 1017 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0005485109058530613 and B = [0.146023,0.961164,0.132146,0.963498]\n", " --- Iteration over ALM coefficient: 10 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 974 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.00034954334404080356 and B = [0.145988,0.961172,0.132105,0.96351]\n", " --- Iteration over ALM coefficient: 11 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 926 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0002264320715512913 and B = [0.14597,0.961177,0.132083,0.963517]\n", " --- Iteration over ALM coefficient: 12 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 874 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.00014875045461057446 and B = [0.14596,0.961179,0.132071,0.96352]\n", " --- Iteration over ALM coefficient: 13 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 820 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 9.885788332636425e-5 and B = [0.145955,0.96118,0.132065,0.963522]\n", " --- Iteration over ALM coefficient: 14 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 765 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 6.631382514887818e-5 and B = [0.145952,0.96118,0.132062,0.963524]\n", " --- Iteration over ALM coefficient: 15 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 708 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.48072540583655e-5 and B = [0.145951,0.961181,0.13206,0.963524]\n", " --- Iteration over ALM coefficient: 16 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 651 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.0443087124187862e-5 and B = [0.14595,0.961181,0.132059,0.963525]\n", " --- Iteration over ALM coefficient: 17 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 594 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.076821253785277e-5 and B = [0.145949,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 18 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 537 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.4209519261981773e-5 and B = [0.145949,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 19 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 480 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 9.741695670306694e-6 and B = [0.145949,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 20 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 424 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 6.687466508698003e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 21 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 370 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.594411339209348e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 22 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 317 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.1576579837266916e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 23 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 266 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.1703967988118134e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 24 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 219 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.4916336017467557e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 25 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 176 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.024884198480569e-6 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 26 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 137 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.039676926667848e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 27 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 103 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.834097304673435e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 28 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 74 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.319196261453161e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 29 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 50 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.2794120549396446e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 30 ---\n", "solving individual problem by Euler equation method\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Euler method converged with 37 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.5634870192959838e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 31 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 27 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.0716623247142287e-7 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 32 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 20 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.337721727451729e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 33 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 18 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 5.00449526541491e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 34 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 15 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.408943094473926e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 35 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 13 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.3183727299036505e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 36 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 10 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.5783928553059212e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 37 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 8 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.0743600203921844e-8 and B = [0.145948,0.961181,0.132058,0.963525]\n", " --- Iteration over ALM coefficient: 38 ---\n", "solving individual problem by Euler equation method\n", "Euler method converged with 7 iterations\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.295628601244886e-9 and B = [0.145948,0.961181,0.132058,0.963525]\n", "-----------------------------------------------------\n", "ALM coefficient successfully converged : dif = 7.295628601244886e-9\n", "-----------------------------------------------------\n", "16333.330134 seconds (20.04 G allocations: 9.153 TB, 27.86% gc time)\n" ] } ], "source": [ "# find ALM coefficient\n", "@time K_ts = find_ALM_coef!(ksp,kss,z_shocks,eps_shocks,\n", " tol_ump=1e-8,max_iter_ump=10000,\n", " Howard_on=true,Howard_n_iter=20,\n", " tol_B=1e-8, max_iter_B=50,update_B=0.3,\n", " T_discard=1000,print_skip_VFI=10,\n", " method=:Euler, update_k=0.7);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare the true aggreate law of motion for capital and approximated one with figure and regression" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_ALM(ksp.z_grid,z_shocks,\n", " kss.B,K_ts)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Approximated aggregate capital law of motion\n", "log(K_{t+1})=0.14594821741362846+0.9611811624862514log(K_{t}) in good time (R2 = 0.9999988558436316)\n", "log(K_{t+1})=0.13205800455894173+0.9635249205659238log(K_{t}) in bad time (R2 = 0.9999972131809477)\n" ] } ], "source": [ "#kss.B # Regression coefficient\n", "println(\"Approximated aggregate capital law of motion\")\n", "println(\"log(K_{t+1})=$(kss.B[1])+$(kss.B[2])log(K_{t}) in good time (R2 =$(kss.R2[1]))\")\n", "println(\"log(K_{t+1})=$(kss.B[3])+$(kss.B[4])log(K_{t}) in bad time (R2 = $(kss.R2[2]))\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The approximated law of motion of capital is very close to the true one, which implies that assuming agents are partially rational is not bad idea since the difference of their actions are negligible.\n", "\n", "Note: The mean of capital, about 40, is sufficiently close to Maliar, Maliar, Valli but higher than Krusell-Smith." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": true }, "outputs": [], "source": [ "save(\"result_Euler.jld\",\"kss\",kss)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean of capital implied by regression is 40.049373762168585\n", "mean of capital implied by simulation is 40.160799927814296\n" ] } ], "source": [ "# Compute mean of capital implied by regression\n", "mc=MarkovChain(ksp.TransMat.Pz)\n", "sd=stationary_distributions(mc)[1]\n", "logKg=kss.B[1]/(1-kss.B[2])\n", "logKb=kss.B[3]/(1-kss.B[4])\n", "meanK_reg=exp(sd[1]*logKg+sd[2]*logKb)\n", "meanK_sim=mean(K_ts[1001:end])\n", "println(\"mean of capital implied by regression is$meanK_reg\")\n", "println(\"mean of capital implied by simulation is $meanK_sim\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figures in Krusell-Smith" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's prepare some functions for figures in Krusell-Smith" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "plot_Fig2 (generic function with 1 method)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function plot_Fig1(ksp,kss,K_ts)\n", " B=kss.B\n", " K_min, K_max = minimum(K_ts), maximum(K_ts)\n", " K_lim=collect(linspace(K_min,K_max,100))\n", " Kp_g=exp(B[1]+B[2]*log(K_lim))\n", " Kp_b=exp(B[3]+B[4]*log(K_lim))\n", " \n", " p=plot(K_lim,Kp_g,linestyle=:solid,lab=\"Good\")\n", " plot!(p,K_lim,Kp_b,linestyle=:solid,lab=\"Bad\")\n", " plot!(p,K_lim,K_lim,color=:black,linestyle=:dash,lab=\"45 degree\",width=0.5)\n", " title!(p,\"FIG1: Tomorrow's vs. today's aggregate capital\")\n", " p\n", "end\n", "\n", "function plot_Fig2(ksp,kss,K_eval_point)\n", " k_lim=collect(linspace(0,80,1000))\n", " itp_e=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,1],Gridded(Linear()))\n", " itp_u=interpolate((ksp.k_grid,ksp.K_grid),kss.k_opt[:,:,3],Gridded(Linear()))\n", " \n", " kp_e(k)=itp_e[k,K_eval_point]\n", " kp_u(k)=itp_u[k,K_eval_point]\n", " \n", " p=plot(k_lim,kp_e.(k_lim),linestyle=:solid,lab=\"employed\")\n", " plot!(p,k_lim,kp_u.(k_lim),linestyle=:solid,lab=\"unemployed\")\n", " plot!(p,k_lim,k_lim,color=:black,linestyle=:dash,lab=\"45 degree\",width=0.5)\n", " title!(p,\"FIG2: Individual policy function \\n at K=$K_eval_point when good state\")\n", " p\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, plot the replication figure:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure 1" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_Fig1(ksp,kss,K_ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure 2" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_Fig2(ksp,kss,40)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Both figures are replicated well. \n", "\n", "Note that the mean of capital is approximately 40 in this replication, which is different from Krusell-Smith but same as Maliar, Maliar, Valli. Therefore, Figure 1 is plotted around 40 and policy function for Figure 2 is evaluated at K=40" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Solution with value function iteration\n", "In this section, each agent's utility maximization problem is solved by value function iteration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's prepare functions for value function iteration\n", "## Bellman equation and VFI\n", "The functions in this cell are used to solve individual household problem by VFI\n", "- Methods\n", " - rhs_bellman: evaluate RHS of bellman equation given the guess of value function\n", " - compute_expectation_VFI: compute expectation term of Bellman equation given current shock state\n", " - solve_bellman_once!: maximize RHS of Bellman equation for an agent of state (k_i,K_i,s_i)\n", " - solve_bellman!: maximize RHS of Bellman equation for all state until convergence\n", " - iterate_policy!: iterating policy and compute values under the policy. used for Howard" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "code_folding": [] }, "outputs": [ { "data": { "text/plain": [ "iterate_policy!" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\"\"\"\n", "Compute right hand side of bellman equation\n", "\n", "##### Arguments\n", "- kp : next period capital\n", "- ksm : KSModel instance\n", "- k : current individual capital\n", "- K : current aggregate capital\n", "- L : current labor\n", "- zeps_i : \n", "\"\"\"\n", "function rhs_bellman(ksp::KSParameter,\n", " kp::AbstractFloat,value::Array{Float64,3},\n", " k::AbstractFloat,K::AbstractFloat,s_i::Integer)\n", " u,s_grid,beta,alpha,l_bar,delta, mu =\n", " ksp.u, ksp.s_grid, ksp.beta, ksp.alpha, ksp.l_bar, ksp.delta, ksp.mu\n", " z, eps = get_shock(s_grid,s_i)\n", " Kp,L = compute_Kp_L(K,s_i,kss.B,ksp) # Next period aggregate capital and current aggregate labor\n", " c = (r(alpha,z,K,L)+1-delta)*k+\n", " w(alpha,z,K,L)*(eps*l_bar+(1.0-eps)*mu)-kp # current consumption \n", " expec = compute_expectation(kp,Kp,value,s_i,ksp)\n", " return u(c)+beta*expec\n", "end\n", "\n", "\"\"\"\n", "Compute expectation of next period value\n", "\n", "##### Arguments\n", "- kp : next period individual capital\n", "- Kp : next period aggregate capital\n", "- value : given value \n", "- s_i : current shock state\n", "- ksp : KSParameter instance\n", "\"\"\"\n", "function compute_expectation(\n", " kp::AbstractFloat, # next period indicidual capital\n", " Kp::AbstractFloat, # next period aggragte capital\n", " value::Array{Float64,3}, # next period value\n", " s_i::Integer, # index of current state,\n", " ksp::KSParameter \n", " )\n", " k_grid, K_grid = ksp.k_grid, ksp.K_grid # unpack grid\n", " beta, P = ksp.beta, ksp.TransMat.P # unpack parameters\n", " \n", " # compute expectations by summing up\n", " expec=0.0\n", " for s_n_i=1:ksp.s_size\n", " value_itp=interpolate((k_grid,K_grid),value[:,:,s_n_i],Gridded(Linear()))\n", " expec = expec + P[s_i,s_n_i]*value_itp[kp,Kp]\n", " end\n", " return expec\n", "end\n", "\"\"\"\n", "Solve bellman equation for all states once\n", "\n", "##### Arguments\n", "- k : individual capital\n", "- K : aggregate capital\n", "- s_i : shock state index\n", "- ksp : KSParameter\n", "- kss : KSSolution\n", "\"\"\"\n", "function solve_bellman_once!(\n", " k_i::Integer,\n", " K_i::Integer,\n", " s_i::Integer,\n", " ksp::KSParameter,\n", " kss::KSSolution,\n", " )\n", " # obtain minimum and maximum of grid\n", " k_min, k_max = ksp.k_grid[1], ksp.k_grid[end]\n", " \n", " # unpack parameters\n", " alpha,delta,l_bar, mu = \n", " ksp.alpha, ksp.delta, ksp.l_bar, ksp.mu\n", " \n", " # obtain state value\n", " k=ksp.k_grid[k_i] # obtain individual capital value\n", " K=ksp.K_grid[K_i] # obtain aggregate capital value\n", " z, eps = get_shock(ksp.s_grid,s_i) # obtain shock value\n", " \n", " Kp,L=compute_Kp_L(K,s_i,kss.B,ksp) # next aggregate capital and current aggregate labor\n", " # if kp>k_c_pos, consumption is negative \n", " k_c_pos=(r(alpha,z,K,L)+1-delta)*k+\n", " w(alpha,z,K,L)*(eps*l_bar+(1-eps)*mu)\n", " obj(kp)=-rhs_bellman(ksp,kp,kss.value,k,K,s_i) # objective function\n", " res=optimize(obj,k_min,min(k_c_pos,k_max)) # maximize value\n", " # obtain result\n", " kss.k_opt[k_i,K_i,s_i]=Optim.minimizer(res) \n", " kss.value[k_i,K_i,s_i]=-Optim.minimum(res)\n", " return nothing\n", "end\n", "\n", "\"\"\"\n", "Solve bellman equation for all states until convergence\n", "\n", "##### Arguments\n", "- ksp : KSParameter\n", "- kss : KSSolution\n", "- tol : tolerance of value function difference\n", "- max_iter : maximum number of iteration\n", "\"\"\"\n", "function solve_bellman!(\n", " ksp::KSParameter,\n", " kss::KSSolution;\n", " tol::AbstractFloat=1e-8,\n", " max_iter::Integer=100,\n", " Howard::Bool=false,\n", " Howard_n_iter::Integer=20,\n", " print_skip::Integer=10)\n", " counter_VFI=0 # counter\n", " while true\n", " counter_VFI += 1\n", " value_old=copy(kss.value) # guessed value\n", " # maximize value for all state\n", " [solve_bellman_once!(k_i,K_i,s_i,ksp,kss)\n", " for k_i in 1:ksp.k_size, K_i in 1:ksp.K_size, s_i in 1:ksp.s_size]\n", " # Howard's policy iteration\n", " !Howard || iterate_policy!(ksp,kss,n_iter=Howard_n_iter)\n", " # difference of guessed and new value\n", " dif=maxabs(value_old-kss.value)\n", " # print covergence process\n", " !(counter_VFI % print_skip ==0) ||\n", " println(\"VFI iteration $counter_VFI : dif =$dif\")\n", " # if difference is sufficiently small\n", " if dif= max_iter\n", " println(\"VFI reached its maximum iteration : $max_iter\")\n", " break\n", " end\n", " end\n", "end\n", "\n", "\"\"\"\n", "Iterate the value function fixing the policy function\n", "\n", "##### Arguments\n", "- ksp : KSParameter instance\n", "- kss : KSSolution instance\n", "- n_iter : number of iterations\n", "\"\"\"\n", "function iterate_policy!(ksp::KSParameter,\n", " kss::KSSolution;n_iter::Int64=20)\n", " value=similar(kss.value)\n", " for i=1:n_iter\n", " # update value using policy\n", " value .= \n", " [rhs_bellman(ksp,\n", " kss.k_opt[k_i,K_i,s_i],kss.value,\n", " ksp.k_grid[k_i],ksp.K_grid[K_i],s_i)\n", " for k_i in 1:ksp.k_size,\n", " K_i in 1:ksp.K_size,\n", " s_i in 1:ksp.s_size]\n", " kss.value.=copy(value)\n", " end\n", " return nothing\n", "end\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Implementation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's skip the following steps in this section to save computational time.\n", "- consturuction of ksp instance since it is same\n", "- consturuction of kss instance to use the previous result as initial guess of the solution\n", "- draws of the shocks to use same ones" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, instead of constructing kss again, obtain value from the policy function derived by Euler method:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": true }, "outputs": [], "source": [ "iterate_policy!(ksp,kss,n_iter=30);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, the model is solved by VFI in the next cell" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " --- Iteration over ALM coefficient: 1 ---\n", "VFI iteration 10 : dif = 4.364130337367499\n", "VFI iteration 20 : dif = 0.5155086466810985\n", "VFI iteration 30 : dif = 0.06101406344856741\n", "VFI iteration 40 : dif = 0.007215234574232454\n", "VFI iteration 50 : dif = 0.000852061281875649\n", "VFI iteration 60 : dif = 0.00010046707080846318\n", "VFI iteration 70 : dif = 1.1826375214241125e-5\n", "VFI iteration 80 : dif = 1.3896425912207633e-6\n", "VFI iteration 90 : dif = 1.6298264426950482e-7\n", "VFI iteration 100 : dif = 1.9079607227467932e-8\n", " ** VFI converged successfully!! dif = 8.085862646112218e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.02551266404864358 and B = [0.130832,0.964894,0.106545,0.970074]\n", " --- Iteration over ALM coefficient: 2 ---\n", "VFI iteration 10 : dif = 0.04110720950990299\n", "VFI iteration 20 : dif = 0.004799456575142358\n", "VFI iteration 30 : dif = 0.0005606688813486471\n", "VFI iteration 40 : dif = 6.551400201715296e-5\n", "VFI iteration 50 : dif = 7.655140279894113e-6\n", "VFI iteration 60 : dif = 8.94252366379078e-7\n", "VFI iteration 70 : dif = 1.0442380471431534e-7\n", "VFI iteration 80 : dif = 1.2187229003757238e-8\n", " ** VFI converged successfully!! dif = 9.830898761720164e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.010571478106761914 and B = [0.138905,0.962971,0.113833,0.968363]\n", " --- Iteration over ALM coefficient: 3 ---\n", "VFI iteration 10 : dif = 0.0004685646893847206\n", "VFI iteration 20 : dif = 5.4314097894803126e-5\n", "VFI iteration 30 : dif = 6.4477072214685904e-6\n", "VFI iteration 40 : dif = 7.622459747835819e-7\n", "VFI iteration 50 : dif = 8.993168876259006e-8\n", "VFI iteration 60 : dif = 1.0597830168990185e-8\n", " ** VFI converged successfully!! dif = 8.557094588468317e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.006449099771311431 and B = [0.140108,0.962642,0.114784,0.968106]\n", " --- Iteration over ALM coefficient: 4 ---\n", "VFI iteration 10 : dif = 0.0004316023300816596\n", "VFI iteration 20 : dif = 5.057148842979586e-5\n", "VFI iteration 30 : dif = 6.002880127198296e-6\n", "VFI iteration 40 : dif = 7.094776037774864e-7\n", "VFI iteration 50 : dif = 8.366657766600838e-8\n", "VFI iteration 60 : dif = 1.0639979564075475e-8\n", " ** VFI converged successfully!! dif = 7.953985914355144e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0037177793034470913 and B = [0.140774,0.962464,0.11558,0.967893]\n", " --- Iteration over ALM coefficient: 5 ---\n", "VFI iteration 10 : dif = 0.0001699352112325414\n", "VFI iteration 20 : dif = 1.9677074277524298e-5\n", "VFI iteration 30 : dif = 2.2995970141437283e-6\n", "VFI iteration 40 : dif = 2.691489555672888e-7\n", "VFI iteration 50 : dif = 3.161997597089794e-8\n", " ** VFI converged successfully!! dif = 8.755819180805702e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0021012252914709983 and B = [0.141146,0.962364,0.116081,0.967759]\n", " --- Iteration over ALM coefficient: 6 ---\n", "VFI iteration 10 : dif = 0.00016965897827958543\n", "VFI iteration 20 : dif = 1.848068291110394e-5\n", "VFI iteration 30 : dif = 2.05708909106761e-6\n", "VFI iteration 40 : dif = 2.325381274204119e-7\n", "VFI iteration 50 : dif = 2.647169594638399e-8\n", " ** VFI converged successfully!! dif = 8.945391982706496e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0012679245607524897 and B = [0.141284,0.962326,0.116284,0.967703]\n", " --- Iteration over ALM coefficient: 7 ---\n", "VFI iteration 10 : dif = 4.9276457787073014e-5\n", "VFI iteration 20 : dif = 4.933374896154419e-6\n", "VFI iteration 30 : dif = 5.170394956621749e-7\n", "VFI iteration 40 : dif = 5.5685518418613356e-8\n", " ** VFI converged successfully!! dif = 9.463576589041622e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0007726858753848709 and B = [0.141347,0.962309,0.116399,0.967672]\n", " --- Iteration over ALM coefficient: 8 ---\n", "VFI iteration 10 : dif = 3.44231706321807e-5\n", "VFI iteration 20 : dif = 3.568709871615283e-6\n", "VFI iteration 30 : dif = 3.8495443277497543e-7\n", "VFI iteration 40 : dif = 4.244049023327534e-8\n", " ** VFI converged successfully!! dif = 9.12990572032868e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.00048209680162884794 and B = [0.141369,0.962302,0.116458,0.967656]\n", " --- Iteration over ALM coefficient: 9 ---\n", "VFI iteration 10 : dif = 1.758010006369659e-5\n", "VFI iteration 20 : dif = 1.7869241446533124e-6\n", "VFI iteration 30 : dif = 1.9031887177334283e-7\n", "VFI iteration 40 : dif = 2.07742800739652e-8\n", " ** VFI converged successfully!! dif = 8.598362910561264e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.0003055378304484685 and B = [0.141377,0.9623,0.11649,0.967647]\n", " --- Iteration over ALM coefficient: 10 ---\n", "VFI iteration 10 : dif = 1.023216435669383e-5\n", "VFI iteration 20 : dif = 1.03666968698235e-6\n", "VFI iteration 30 : dif = 1.102799274121935e-7\n", "VFI iteration 40 : dif = 1.2026134754705708e-8\n", " ** VFI converged successfully!! dif = 9.64348600973608e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.00019619175347125595 and B = [0.141377,0.9623,0.116508,0.967642]\n", " --- Iteration over ALM coefficient: 11 ---\n", "VFI iteration 10 : dif = 6.056440042812028e-6\n", "VFI iteration 20 : dif = 6.10131223766075e-7\n", "VFI iteration 30 : dif = 6.479382363977493e-8\n", " ** VFI converged successfully!! dif = 8.801350759313209e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 0.00012730360336368762 and B = [0.141376,0.962301,0.116518,0.96764]\n", " --- Iteration over ALM coefficient: 12 ---\n", "VFI iteration 10 : dif = 3.566279190181376e-6\n", "VFI iteration 20 : dif = 3.556169758667238e-7\n", "VFI iteration 30 : dif = 3.754672661671066e-8\n", " ** VFI converged successfully!! dif = 9.877851425699191e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 8.313458897782344e-5 and B = [0.141374,0.962301,0.116524,0.967638]\n", " --- Iteration over ALM coefficient: 13 ---\n", "VFI iteration 10 : dif = 2.2697702775076323e-6\n", "VFI iteration 20 : dif = 2.2698196744386223e-7\n", "VFI iteration 30 : dif = 2.4034420675889123e-8\n", " ** VFI converged successfully!! dif = 9.869268069451209e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 5.463237199027082e-5 and B = [0.141373,0.962302,0.116527,0.967637]\n", " --- Iteration over ALM coefficient: 14 ---\n", "VFI iteration 10 : dif = 1.2852221971115796e-6\n", "VFI iteration 20 : dif = 1.2533530480141053e-7\n", "VFI iteration 30 : dif = 1.3040221347182523e-8\n", " ** VFI converged successfully!! dif = 8.32721980259521e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.63293125280606e-5 and B = [0.141371,0.962302,0.116529,0.967636]\n", " --- Iteration over ALM coefficient: 15 ---\n", "VFI iteration 10 : dif = 9.561501315147325e-7\n", "VFI iteration 20 : dif = 9.57708721216477e-8\n", "VFI iteration 30 : dif = 1.017394879454514e-8\n", " ** VFI converged successfully!! dif = 8.143899776769103e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.4483813326064974e-5 and B = [0.14137,0.962302,0.116531,0.967636]\n", " --- Iteration over ALM coefficient: 16 ---\n", "VFI iteration 10 : dif = 5.508533149622963e-7\n", "VFI iteration 20 : dif = 5.388682211560081e-8\n", " ** VFI converged successfully!! dif = 8.823121788736898e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.636877019500771e-5 and B = [0.14137,0.962302,0.116532,0.967636]\n", " --- Iteration over ALM coefficient: 17 ---\n", "VFI iteration 10 : dif = 4.021055701741716e-7\n", "VFI iteration 20 : dif = 4.010468046544702e-8\n", " ** VFI converged successfully!! dif = 8.30476665214519e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.0984150278420257e-5 and B = [0.141369,0.962302,0.116532,0.967636]\n", " --- Iteration over ALM coefficient: 18 ---\n", "VFI iteration 10 : dif = 2.3281569383470924e-7\n", "VFI iteration 20 : dif = 2.2501069452118827e-8\n", " ** VFI converged successfully!! dif = 9.03918362382683e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.2048826815218625e-6 and B = [0.141369,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 19 ---\n", "VFI iteration 10 : dif = 1.230218344971945e-7\n", "VFI iteration 20 : dif = 1.0416044915473321e-8\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ " ** VFI converged successfully!! dif = 8.250196970038814e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.988788248755371e-6 and B = [0.141369,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 20 ---\n", "VFI iteration 10 : dif = 1.0984092568833148e-7\n", "VFI iteration 20 : dif = 1.089637180484715e-8\n", " ** VFI converged successfully!! dif = 8.687948138685897e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.290747322109988e-6 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 21 ---\n", "VFI iteration 10 : dif = 8.936746098697768e-8\n", "VFI iteration 20 : dif = 9.178563686873531e-9\n", " ** VFI converged successfully!! dif = 9.178563686873531e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.158295437842961e-6 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 22 ---\n", "VFI iteration 10 : dif = 3.8132498048071284e-8\n", " ** VFI converged successfully!! dif = 9.028951808431884e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.5509704746619057e-6 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 23 ---\n", "VFI iteration 10 : dif = 1.3471407100951183e-8\n", " ** VFI converged successfully!! dif = 8.787935712462058e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.0163098228543888e-6 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 24 ---\n", "VFI iteration 10 : dif = 4.011184273622348e-8\n", " ** VFI converged successfully!! dif = 8.319318567373557e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 6.248587137436257e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 25 ---\n", "VFI iteration 10 : dif = 2.6097382033185568e-8\n", " ** VFI converged successfully!! dif = 8.564882136852248e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.880477443818364e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 26 ---\n", "VFI iteration 10 : dif = 2.0421111912583e-8\n", " ** VFI converged successfully!! dif = 8.35751734484802e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.843728103413268e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 27 ---\n", "VFI iteration 10 : dif = 3.812928639490565e-8\n", " ** VFI converged successfully!! dif = 8.57767190609593e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 2.919172302912054e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 28 ---\n", "VFI iteration 10 : dif = 1.3114174635120435e-8\n", " ** VFI converged successfully!! dif = 8.181245902960654e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.6344598458006843e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 29 ---\n", "VFI iteration 10 : dif = 5.4228507906373125e-8\n", " ** VFI converged successfully!! dif = 9.760071861819597e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 9.524729271959131e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 30 ---\n", "VFI iteration 10 : dif = 3.021312977580237e-8\n", " ** VFI converged successfully!! dif = 8.285155672638211e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.177423369530977e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 31 ---\n", "VFI iteration 10 : dif = 2.392872033851745e-8\n", " ** VFI converged successfully!! dif = 8.194376732717501e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.0872591149624355e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 32 ---\n", "VFI iteration 10 : dif = 5.144397619005758e-8\n", " ** VFI converged successfully!! dif = 9.24671894608764e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.570228442410171e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 33 ---\n", "VFI iteration 10 : dif = 2.3067457277647918e-8\n", " ** VFI converged successfully!! dif = 9.729632211019634e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.874108731709129e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 34 ---\n", "VFI iteration 10 : dif = 7.674799462620285e-8\n", "VFI iteration 20 : dif = 8.989559319161344e-9\n", " ** VFI converged successfully!! dif = 8.989559319161344e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.674765916561462e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 35 ---\n", "VFI iteration 10 : dif = 3.5518326058081584e-8\n", " ** VFI converged successfully!! dif = 9.848633908404736e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.1407666089535695e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 36 ---\n", "VFI iteration 10 : dif = 5.78716878862906e-8\n", " ** VFI converged successfully!! dif = 8.39969516164274e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 3.834874243158204e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 37 ---\n", "VFI iteration 10 : dif = 8.088653657978284e-8\n", "VFI iteration 20 : dif = 9.574307568982476e-9\n", " ** VFI converged successfully!! dif = 9.574307568982476e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.068997078165992e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 38 ---\n", "VFI iteration 10 : dif = 6.087861947889905e-8\n", " ** VFI converged successfully!! dif = 8.969834652816644e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.422016042227028e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 39 ---\n", " ** VFI converged successfully!! dif = 9.783235555005376e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.2206419956750647e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 40 ---\n", "VFI iteration 10 : dif = 7.822649195077247e-8\n", "VFI iteration 20 : dif = 9.270820555684622e-9\n", " ** VFI converged successfully!! dif = 9.270820555684622e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.8181412958506726e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 41 ---\n", "VFI iteration 10 : dif = 1.760491841196199e-8\n", " ** VFI converged successfully!! dif = 9.28258714338881e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 5.6862735645091256e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 42 ---\n", "VFI iteration 10 : dif = 5.217941634327872e-8\n", " ** VFI converged successfully!! dif = 9.463349215366179e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 9.779224406647469e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 43 ---\n", "VFI iteration 10 : dif = 2.28354224418581e-8\n", " ** VFI converged successfully!! dif = 9.649710364101338e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 4.013947332848211e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 44 ---\n", " ** VFI converged successfully!! dif = 7.641517640877282e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 8.247222312018909e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 45 ---\n", "VFI iteration 10 : dif = 5.345765430320171e-8\n", " ** VFI converged successfully!! dif = 9.702773695607902e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 5.5271162438530475e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 46 ---\n", "VFI iteration 10 : dif = 6.705204214085825e-8\n", " ** VFI converged successfully!! dif = 9.735970252222614e-9\n", "simulating aggregate path ... please wait ... \n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "difference of ALM coefficient is 5.700652201678924e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 47 ---\n", "VFI iteration 10 : dif = 2.999195203301497e-8\n", " ** VFI converged successfully!! dif = 8.3674080997298e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.5016738530437834e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 48 ---\n", "VFI iteration 10 : dif = 1.08337872006814e-8\n", " ** VFI converged successfully!! dif = 8.718018307263264e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 1.779748049768326e-7 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 49 ---\n", "VFI iteration 10 : dif = 1.6313094874931267e-8\n", " ** VFI converged successfully!! dif = 8.512699878338026e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 6.603349556044691e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", " --- Iteration over ALM coefficient: 50 ---\n", "VFI iteration 10 : dif = 7.444867833328317e-8\n", "VFI iteration 20 : dif = 8.854442512529204e-9\n", " ** VFI converged successfully!! dif = 8.854442512529204e-9\n", "simulating aggregate path ... please wait ... \n", "difference of ALM coefficient is 7.825109485382065e-8 and B = [0.141368,0.962303,0.116533,0.967635]\n", "----------------------------------------------------------------\n", "Iteration over ALM coefficient reached its maximum (50)\n", "----------------------------------------------------------------\n", "26895.406260 seconds (27.08 G allocations: 12.690 TB, 25.15% gc time)\n" ] } ], "source": [ "# find ALM coefficient\n", "@time K_ts = find_ALM_coef!(ksp,kss,z_shocks,eps_shocks,\n", " tol_ump=1e-8,max_iter_ump=10000,\n", " Howard_on=true,Howard_n_iter=20,\n", " tol_B=1e-8, max_iter_B=50,update_B=0.3,\n", " T_discard=1000,print_skip_VFI=10,\n", " method=:VFI, update_k=0.7);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The following same exercises show that the main result is same as before" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_ALM(ksp.z_grid,z_shocks,\n", " kss.B,K_ts)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Approximated aggregate capital law of motion\n", "log(K_{t+1})=0.14136817405539825+0.9623027152453734log(K_{t}) in good time (R2 = 0.9999992137407112)\n", "log(K_{t+1})=0.11653324573318116+0.9676352636916012log(K_{t}) in bad time (R2 = 0.9999982189177621)\n" ] } ], "source": [ "#kss.B # Regression coefficient\n", "println(\"Approximated aggregate capital law of motion\")\n", "println(\"log(K_{t+1})=$(kss.B[1])+$(kss.B[2])log(K_{t}) in good time (R2 =$(kss.R2[1]))\")\n", "println(\"log(K_{t+1})=$(kss.B[3])+$(kss.B[4])log(K_{t}) in bad time (R2 = $(kss.R2[2]))\")" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true }, "outputs": [], "source": [ "save(\"result_VFI.jld\",\"kss\",kss)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean of capital implied by regression is 39.46272307237643\n", "mean of capital implied by simulation is 39.70193508699299\n" ] } ], "source": [ "# Compute mean of capital implied by regression\n", "mc=MarkovChain(ksp.TransMat.Pz)\n", "sd=stationary_distributions(mc)[1]\n", "logKg=kss.B[1]/(1-kss.B[2])\n", "logKb=kss.B[3]/(1-kss.B[4])\n", "meanK_reg=exp(sd[1]*logKg+sd[2]*logKb)\n", "meanK_sim=mean(K_ts[1001:end])\n", "println(\"mean of capital implied by regression is$meanK_reg\")\n", "println(\"mean of capital implied by simulation is \$meanK_sim\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figures in Krusell-Smith" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure 1" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_Fig1(ksp,kss,K_ts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure 2" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot_Fig2(ksp,kss,40)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, the figures are succesfully replicated." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.1", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.1" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "31px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": false, "toc_section_display": "block", "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }