{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using BenchmarkTools\n", "using Plots" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "VFI (generic function with 1 method)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "type Model\n", " # primitive parameter\n", " beta::Float64 #subjective discount factor\n", " sigma::Float64 # relative riskb aversion\n", " delta::Float64 #depriciation rate\n", " alpha::Float64 # capital share\n", "\n", " # discretize asset space\n", " agrid::Vector{Float64}\n", "end\n", "\n", "\n", "function argmax(mat)\n", " values, indices = findmax(mat,2)\n", " return ind2sub(size(mat),vec(indices))[2]\n", "end\n", "\n", "\n", "function VFI(m::Model; max_iter::Int=1000, tol::Float64=1e-5)\n", " const penalty = -999999999.9\n", " const na = size(m.agrid, 1)\n", " \n", " #initialize value function and so on\n", " v = zeros(na, na) # temp value function\n", " c = zeros(na, na) # consuption matrix\n", " util = zeros(na, na) # utility matrix\n", " \n", " v0 = zeros(na, 1) # initial guess of value function\n", " Tv = zeros(na, 1) # update value function\n", " pol_a = zeros(na, 1) # policy function\n", " \n", " #create consuption and utility matrix \n", " for i in 1:na\n", " for j in 1:na\n", " @inbounds c[j,i] = m.agrid[j]^m.alpha + (1-m.delta) * m.agrid[j] - m.agrid[i]\n", " @inbounds util[j,i] =(c[j,i]^(1-m.sigma)) / (1-m.sigma)\n", " if c[j,i] <= 0\n", " @inbounds util[j,i] = penalty # penalty\n", " end\n", " end\n", " end\n", " \n", " # value function iteration\n", " err = 0.0\n", " for t in 1:max_iter\n", " # calculate temp value function\n", " for i in 1:na\n", " for j in 1:na\n", " @inbounds v[j, i] = util[j, i] + m.beta * v0[i]\n", " end\n", " end\n", " \n", " Tv = maximum(v, 2) # obtain new value funtion\n", " err = maximum(abs.(Tv - v0)) # update error\n", " v0 = Tv # update value function\n", " \n", " if err < tol\n", " # obtain policy function\n", " a_index = argmax(v)\n", " for i in 1:na\n", " @inbounds pol_a[i] = m.agrid[a_index[i]]\n", " end\n", " break\n", " end\n", " end\n", " \n", " if err >= tol\n", " println(\"VFI does not converge in $(max_iter) times\")\n", " end\n", " \n", " return(m.agrid, v0, pol_a)\n", "end" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Model(0.95, 2.0, 0.1, 0.33, [0.316086, 0.340205, 0.364324, 0.388443, 0.412562, 0.436681, 0.4608, 0.484919, 0.509038, 0.533157 … 6.10465, 6.12877, 6.15289, 6.17701, 6.20113, 6.22524, 6.24936, 6.27348, 6.2976, 6.32172])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta = 0.95 #subjective discount factor\n", "sigma = 2.0 # relative riskb aversion\n", "delta = 0.1 #depriciation rate\n", "alpha = 0.33 # capital share\n", "\n", "# Steady state\n", "aterm = 1.0/beta -(1.0 -delta)\n", "kstar = alpha/aterm\n", "kstar = kstar^(1.0/(1.0-alpha))\n", "\n", "amin = 0.1 * kstar\n", "amax = 2 * kstar\n", "na = 250\n", "\n", "model = Model(beta, sigma, delta, alpha, linspace(amin, amax, na))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "([0.316086, 0.340205, 0.364324, 0.388443, 0.412562, 0.436681, 0.4608, 0.484919, 0.509038, 0.533157 … 6.10465, 6.12877, 6.15289, 6.17701, 6.20113, 6.22524, 6.24936, 6.27348, 6.2976, 6.32172], [-22.9349; -22.7696; … ; -15.6818; -15.6718], [0.484919; 0.509038; … ; 5.95994; 5.95994])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "agrid, v0, pol_a = VFI(model)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "BenchmarkTools.Trial: \n", " memory estimate: 2.87 MiB\n", " allocs estimate: 1828\n", " --------------\n", " minimum time: 24.128 ms (0.00% GC)\n", " median time: 25.455 ms (0.00% GC)\n", " mean time: 26.098 ms (0.77% GC)\n", " maximum time: 33.859 ms (4.05% GC)\n", " --------------\n", " samples: 192\n", " evals/sample: 1" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@benchmark VFI(model)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", " \n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(agrid, pol_a, color=\"blue\", linewidth=1.5, label=\"Policy Function\")\n", "plot!(agrid, agrid, color=\"red\", linewidth=1.5, label=\"45 degree Line\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ " \n", " \n" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(agrid, v0, color=\"blue\", linewidth=1.5, label=\"Value Function\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Variables:\n", " #self#::#VFI\n", " m::Model\n", "\n", "Body:\n", " begin \n", " return $(Expr(:invoke, MethodInstance for #VFI#1(::Int64, ::Float64, ::Function, ::Model), :(Main.#VFI#1), 1000, 1.0e-5, :(#self#), :(m)))\n", " end::Tuple{Array{Float64,1},Array{Float64,2},Array{Float64,2}}\n" ] } ], "source": [ "@code_warntype VFI(model)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0.028092 seconds (1.84 k allocations: 2.868 MiB)\n" ] }, { "data": { "text/plain": [ "([0.316086, 0.340205, 0.364324, 0.388443, 0.412562, 0.436681, 0.4608, 0.484919, 0.509038, 0.533157 … 6.10465, 6.12877, 6.15289, 6.17701, 6.20113, 6.22524, 6.24936, 6.27348, 6.2976, 6.32172], [-22.9349; -22.7696; … ; -15.6818; -15.6718], [0.484919; 0.509038; … ; 5.95994; 5.95994])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@time VFI(model)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.2", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.2" } }, "nbformat": 4, "nbformat_minor": 2 }