{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 動的計画法を速く解こう！\n",
"\n",
"** とくぷ(Twitter:@tokutokupooh, GitHub:tokupooh)**\n",
"\n",
"[6日目](https://docs.google.com/document/d/160xW45I2v-MUkJF8TyaXhENCLJA243WOGdjK-G22UQM/edit)→ 7日目\n",
"\n",
"この記事は[みゅーもり Advent Calendar](https://adventar.org/calendars/3429)の7日目の記事になります。\n",
"\n",
"本記事では様々な場面で用いられている動的計画法を具体例に基づき、解法について概説します。\n",
"\n",
"\n",
"具体的には\n",
"1. Value Function Iteration(VFI)\n",
"2. Endogenous Grid Method(EGM)\n",
"\n",
"について仕組みと実装を紹介します。\n",
"\n",
"\n",
"例として取り上げるモデルについては詳細は省きますが、\n",
"さらに知りたい方については\n",
"[Ljungqvist and Sargent(2018)](https://mitpress.mit.edu/books/recursive-macroeconomic-theory-fourth-edition)\n",
"の17章がおすすめです。\n",
"\n",
"\n",
"前提知識としては\n",
"- 簡単な微分の知識\n",
"- 動的計画法とはなにか\n",
"\n",
"のみを仮定します。\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## モデル\n",
"\n",
"この節では、本記事で注目するモデルについて概説します。\n",
"\n",
"\n",
"次のような最大化問題を考えます。\n",
"\n",
"### 定式化\n",
"\\begin{equation*}\n",
"\\max_{\\{c_t,a_{t+1}\\}_{t=0}^{\\infty}}\\sum_{t=0}^{\\infty} \\beta^t u(c_t)\n",
"\\end{equation*}\n",
"subject to \n",
"\\begin{equation*}\n",
"c_t + a_{t+1} = (1+r)a_t + wz_t\n",
"\\end{equation*}\n",
"\\begin{equation*}\n",
"a_{t+1} \\ge 0\n",
"\\end{equation*}\n",
"\\begin{equation*}\n",
"\\log z_{t+1} = \\rho \\log z_t + \\epsilon_t\n",
"\\end{equation*}\n",
"\n",
"- $u$: 効用関数\n",
"\n",
"- $c$: 消費\n",
"- $a$: 資産\n",
"- $\\beta$: 割引因子\n",
"- $r$: 利子率\n",
"- $w$: 賃金\n",
"- $z$: 個人の生産性\n",
"- $\\rho$: 生産性ショックの持続性\n",
"- $\\epsilon$: 平均0,分散$\\simga^2$の正規分布に従うランダムな要素 \n",
"\n",
"\n",
"直観的にはこのモデルは次のように考えることができます。\n",
"\n",
"無限期間生きる家計が生涯における割引期待効用を最大にするように消費、貯蓄の計画を決定します。\n",
"\n",
"家計は毎期一定の所得を手に入れることができますが、そうした所得はランダムに決まる自らの生産性に依存します。\n",
"\n",
"所得がランダムに決まるので、所得が少ない時期が長く続く場合は借金をして消費を増やそうとするかもしれません。\n",
"しかし借金がここでは許されていないので、手に入れた所得をすべて消費に回します。\n",
"\n",
"### ベルマン方程式\n",
"\n",
"以上のような流列問題はベルマン方程式を用いて定式化しなおすことができます。ベルマン方程式は次のようになります。\n",
"\n",
"\\begin{equation*}\n",
"v(a,z) = \\max_{a'}\\left[u(c) + \\beta E v(a',z') \\right]\n",
"\\end{equation*}\n",
"subject to\n",
"\n",
"\\begin{equation*}\n",
"c + a' = (1+r)a + wz\n",
"\\end{equation*}\n",
"\\begin{equation*}\n",
"a_{t+1} \\ge 0\n",
"\\end{equation*}\n",
"\\begin{equation*}\n",
"\\log z' = \\rho \\log z + \\epsilon\n",
"\\end{equation*}\n",
"\n",
"\n",
"### Parametarization\n",
"\n",
"ここでこのモデルを計算するためにパラメータを特定化しましょう。\n",
"\n",
"- $u(c) = \\frac{c^{1-\\sigma}}{1-\\sigma}$\n",
"- $\\beta=0.96$\n",
"- $\\sigma=3.0$: 相対的リスク回避度\n",
"- $\\phi =0.0$: 借入制約\n",
"\n",
"AR(1)プロセスは取り扱いのしやすいマルコフプロセスによって近似します。\n",
"詳細についてはTauchen(1986), Rouwenhorst(1995)をご覧ください。\n",
"\n",
"ここでは生産性のグリッドを$z$,生産性の推移行列を$\\Pi$とし、\n",
"AR(1)プロセスが以下のマルコフプロセスで記述できるとします。\n",
"\n",
"- $z = [0.2,1.0]$\n",
"- $\\Pi = \\begin{bmatrix}0.7 & 0.3\\\\ 0.1 & 0.9\\end{bmatrix}$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Value Function Iteration\n",
"\n",
"VFIは縮小写像の定理より収束することが知られており、安定的に動的計画法を解くことができます。\n",
"\n",
"ただしモデルの構造などを利用していないことから計算速度は遅いという欠点があります。\n",
"\n",
"\n",
"### アイディア\n",
"\n",
"\n",
"Value functionを予想したうえで、その予想値のもとでValue functionを繰り返し更新していくことで解を求めます。\n",
"\n",
"Value functionはベルマン方程式を利用して次のように更新していきます。\n",
"\n",
"\\begin{equation*}\n",
"v^{(n)}(a,z) = \\max_{a'} \\left[u(c) + \\beta Ev^{(n-1)}(a',z')\\right]\n",
"\\end{equation*}\n",
"\n",
"今回はVFIの仕組みを概観することが目的なので、最も単純なVFIについて紹介します。\n",
"\n",
"### アルゴリズム\n",
"\n",
"1. 今期の資産$a$と来期の資産$a'$をそれぞれ$n$個に区切る $\\{a_i,a'_j\\}_{i,j=1}^n$\n",
"2. 価値関数の初期値を定める\n",
"3. 収束の基準を定める\n",
"4. 効用関数を各グリッドポイントにおいて計算する\n",
"5. 各状態変数のグリッドについて $v^{(n)}(a,z) = \\max_{a'} \\left[u(c) + \\beta Ev^{(n-1)}(a',z')\\right]$ を計算する\n",
"6. $v^{(n)}$と$v^{(n-1)}$の誤差が十分小さくなれば、計算を終わらせる。そうでなければ価値関数の初期値を$v^{(n)}$と設定し、5に戻る。\n",
"\n",
"\n",
"以下ではアルゴリズムについて補足を行います。\n",
"\n",
"**step 4**\n",
"\n",
"効用関数を各グリッドポイントで計算しますが、ここでは$a,a',z$のすべての組について計算する必要があります。\n",
"\n",
"そのため効用関数の行列は$n\\times n \\times 2$となります。\n",
"\n",
"**step 5**\n",
"\n",
"ステップ5では具体的には状態変数の組について\n",
"\\begin{equation*}\n",
"v^{(n)}(a_i,z_h)= \\max_{a_j}\\{u((1+r)a_i + wz_h- a'_j)+ \\beta \\sum_{l=1}^{2} \\Pi_{jl}v^{(n-1)}(a' _j z'_l)\\} \n",
"\\end{equation*}\n",
"を計算していきます。\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Plots.GRBackend()"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#=利用するパッケージ=#\n",
"using LinearAlgebra\n",
"using BenchmarkTools\n",
"using Plots\n",
"gr() "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"struct Model{TI<:Integer, TF<:AbstractFloat,TV<:Vector,TM<:Matrix}\n",
" \n",
" # EGM, VFIともに同じ設定で解くので共通部分をまとめて使いやすくする\n",
"\n",
" β::TF # 割引因子\n",
" σ::TF # 相対的リスク回避度\n",
" ϕ::TF # 正の値を取る借入制約(どれだけ借りられるかを示す)\n",
" \n",
" a2grid::TV #Array{TF,1} # 資産グリッド\n",
" na::TI # 資産グリッドのサイズ\n",
" \n",
" zgrid::TV #Array{TF,1} # ショックのリッド\n",
" nz::TI # ショックグリッドのサイズ\n",
" Π::TM #Array{TF,2} # 推移行列\n",
"\n",
" maxiter::TI # 反復上限回数\n",
" tol::TF # 誤差許容度\n",
"\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Construct (generic function with 1 method)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function Construct(;β=0.96,\n",
" σ=3.0,\n",
" ϕ=0.0,\n",
" a2min=0.0,\n",
" a2max=10.0, \n",
" na=500) \n",
" \n",
" # 借入制約を含んだグリッドを作る\n",
" # ここでは借入制約付近にグリッドポイントが多くなるように作成する\n",
" dist_grid = (a2max-(a2min- ϕ))/(na-1).^2 # グリッドポイント間の距離\n",
" indices = collect(Int64,1:na) # グリッドのインデックス\n",
" a2grid = (a2min-ϕ) .+ dist_grid*(indices.-1).^2 # 資産グリッドの形成\n",
" \n",
" # ショックのグリッドを作る\n",
" zgrid = [0.2, 1.0] # 生産性ショック\n",
" nz = length(zgrid)\n",
" Π =[0.7 0.3; 0.1 0.9] # 推移行列 \n",
" \n",
" # 反復に関する設定\n",
" maxiter = 10000 # 反復上限回数\n",
" tol = 1e-13 # 誤差許容度\n",
" \n",
" return Model(β,σ,ϕ,a2grid,na,zgrid,nz,Π,maxiter,tol) \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"Model{Int64,Float64,Array{Float64,1},Array{Float64,2}}(0.96, 3.0, 0.0, [0.0, 4.01605e-5, 0.000160642, 0.000361444, 0.000642568, 0.00100401, 0.00144578, 0.00196786, 0.00257027, 0.003253 … 9.64253, 9.68193, 9.72141, 9.76096, 9.8006, 9.84032, 9.88012, 9.92, 9.95996, 10.0], 500, [0.2, 1.0], 2, [0.7 0.3; 0.1 0.9], 10000, 1.0e-13)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = Construct() # パラメータなどを含んだオブジェクトを作った"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"code_folding": []
},
"outputs": [
{
"data": {
"text/plain": [
"VFI (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@inline function VFI(m::Model)\n",
" \n",
" # 要素価格\n",
" r = 0.03 # 利子率\n",
" w = 1.0 # 賃金\n",
" \n",
" # ペナルティ\n",
" neg = -1e10\n",
" # 効用関数\n",
" @inline u(c,σ) = c > 1e-10 ? (c^(1.0-σ))/(1.0-σ) : neg\n",
" \n",
" \n",
" #=利用する変数の初期化=#\n",
" util = zeros(m.na,m.na,m.nz)\n",
" v0 = ones(m.na,m.nz) #価値関数の初期予想\n",
" Tv = similar(v0) #アップデートされた価値関数\n",
" DEV = similar(v0) # \\beta * E V(a',z')の値を入れる行列\n",
" pol_a = similar(v0) # 資産に関する政策関数\n",
" pol_c = similar(v0) # 消費に関する政策関数\n",
" \n",
" # 効用行列の計算 step 4\n",
" @inbounds for z in eachindex(m.zgrid)\n",
" @inbounds for j in eachindex(m.a2grid)\n",
" @inbounds for i in eachindex(m.a2grid)\n",
" c = w*m.zgrid[z] + (1.0+r)*m.a2grid[i] - m.a2grid[j]\n",
" util[i,j,z] = u(c, m.σ)\n",
" end\n",
" end\n",
" end\n",
" \n",
" #= VFIの実行 =#\n",
" for VFI_it in 1: m.maxiter \n",
" \n",
" #=step 5=#\n",
" DEV = m.β*v0*m.Π' #\\beta * E V(a',z')の部分を先に計算しておく\n",
" @inbounds for z in eachindex(m.zgrid)\n",
" @inbounds for i in eachindex(m.a2grid)\n",
" vmin = neg\n",
" @inbounds for j in eachindex(m.a2grid)\n",
" v = util[i,j,z] + DEV[j,z]\n",
" if v > vmin # 計算のたびに大きいものを代入していく 最適化ステップ\n",
" vmin = v\n",
" Tv[i,z] = vmin\n",
" pol_a[i,z] = m.a2grid[j]\n",
" end\n",
" end\n",
" # 政策関数を予算制約から計算\n",
" pol_c[i,z] = w*m.zgrid[z] +(1.0+r)*m.a2grid[i] - pol_a[i,z]\n",
" end\n",
" end\n",
" \n",
" #=step 6 =#\n",
" err = norm(Tv-v0, Inf) # 誤差の確認\n",
" v0 = copy(Tv) # アップデートしたものを予想のものに代入\n",
" if err < m.tol #十分に誤差が小さくなったかを確認\n",
" break\n",
" end\n",
" \n",
" if VFI_it == m.maxiter\n",
" println(\"The model does not converge...\")\n",
" end\n",
" end \n",
" \n",
" return v0, pol_a,pol_c\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([-61.5264 -26.669; -61.5212 -26.6687; … ; -11.7195 -10.4462; -11.6891 -10.4265], [0.0 0.451243; 0.0 0.451243; … ; 9.369 10.0; 9.40783 10.0], [0.2 0.548757; 0.200041 0.548798; … ; 1.08976 1.25876; 1.09217 1.3])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v, pol_a_VFI, pol_c_VFI = VFI(m)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: \n",
" memory estimate: 26.55 MiB\n",
" allocs estimate: 2936\n",
" --------------\n",
" minimum time: 1.503 s (0.32% GC)\n",
" median time: 1.552 s (0.74% GC)\n",
" mean time: 1.578 s (1.81% GC)\n",
" maximum time: 1.704 s (1.05% GC)\n",
" --------------\n",
" samples: 4\n",
" evals/sample: 1"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark VFI(m)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"labels=[\"Bad State\",\"Good State\"]\n",
"plot(m.a2grid,pol_a_VFI, label=labels,legend=:bottomright)\n",
"plot!(m.a2grid,m.a2grid, label =\"45 degree\",legend=:bottomright)\n",
"title!(\"Policy function for asset using VFI\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"このようにVFIでは約1.6秒でこのような政策関数を得ることが出来ました。 \n",
"\n",
"政策関数のプロットをみると生産性が低く、資産が少ない場合は来期に向けた貯蓄をしない様子が見受けられます。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Endogenous Grid Method\n",
"\n",
"\n",
"次にEndogenous Grid Method(EGM)について紹介します。\n",
"\n",
"VFIは操作変数に関して最適化を行うステップが存在し、この部分の計算量が状態変数の数に応じて非常に多くなってしまいます。\n",
"\n",
"特に操作変数を連続値として解くようなステップではnumerical solverを状態変数のグリッドポイントごとに使います。\n",
"\n",
"そのためVFIよりも計算時間がかかることが予想されます。\n",
"\n",
"しかしモデルの構造を十分に生かすことでnumerical solverの利用を回避し、計算速度を大幅に上げることができます。\n",
"\n",
"### アイディア\n",
"\n",
"前述のモデルの1階条件から求められるオイラー方程式は次のようになります。\n",
"\n",
"\\begin{equation*}\n",
"u'(c_t)= \\beta(1+r)E_t u'(c_{t+1})\n",
"\\end{equation*}\n",
"\n",
"仮定より$u'$は逆関数が取れるので、$c_t = u^{'-1}\\left( \\beta (1+r)E_t u'(c_{t+1})\\right)$とすることができます。\n",
"\n",
"ここからわかることは$c_{t+1}$の値が分かっていれば自動的に$c_t$が決定されるということです。\n",
"\n",
"したがってEGMではt+1期についてのグリッドを作ることがポイントになっています。\n",
"\n",
"\n",
"\n",
"### アルゴリズム\n",
"\n",
"1. 来期期の資産,a'をn個に区切る。\n",
"2. 消費に関する政策関数,$c' = g_0(a',z')$を予想する。\n",
"3. 収束の基準を定める。\n",
"4. $c = u^{'-1}\\left( \\beta (1+r)E u'(g_0(a',z'))\\right)$を利用して今期の消費$\\text{ctil} $を求める\n",
"5. 4.で求めた消費と$a'$のグリッドをもとに$a$のグリッドを求める\n",
"6. 4,5から得られる消費関数$\\text{ctil}= g(a,z) $を$a'$で評価しなおす。\n",
"7. 消費関数の誤差が十分小さくなれば、計算を終わらせる。そうでなければ消費関数をアップデートし、4に戻る。\n",
"\n",
"以下ではアルゴリズムについていくつかの補足を行います。\n",
"\n",
"説明のために a'のグリッドを$a'=\\{a'_1,...,a'_{na}\\}$,zのグリッドを$z= \\{z_1,...,z_{nz}\\}$とします。\n",
"\n",
"\n",
"**step 4**\n",
"\n",
"オイラー方程式から今期の消費量を求めますが、具体的には\n",
"\n",
"\\begin{equation*}\n",
"\\text{ctil}_{ij} = u^{'-1}\\left(\\beta(1+r) \\sum_{l = 1}^{nz} \\Pi_{jl} u'(g_0(a'_i,z_l))\\right)\n",
"\\end{equation*}\n",
"\n",
"とすることで各要素を計算します。\n",
"\n",
"**step 5**\n",
"\n",
"ステップ5で求められた今期の資産を$a^\\ast$とします。\n",
"\n",
"予算制約から$a^\\ast$の各要素は次のように求めることができます。\n",
"\n",
"\\begin{equation*}\n",
"a^\\ast_{ij} = \\frac{a'_i + ctil_{ij} - wz_j}{1+r}\n",
"\\end{equation*}\n",
"\n",
"ここで分かるように$a^\\ast$は行列となります。VFIではt期の資産のグリッドはベクトルですから、この点は大きく異なります。\n",
"\n",
"ここで$a'_1= -\\phi$であるので、$a^\\ast_{1j}$は借入制約がバインドするような最大の今期の資産となります。\n",
"\n",
"\n",
"**step 6**\n",
"\n",
"ステップ4,5から消費関数は$\\text{ctil}_{ij}= g(a^\\ast_{ij},z_j)$として得られます。\n",
"\n",
"ただしここで得られた消費関数は$a^\\ast$と$z$上で定義されていますが、オリジナルの消費関数$g_0(a'_i,z_j)$は$a'$と$z$上で定義されています。\n",
"\n",
"そこで$g(a^\\ast_{ij},z_j)$を$a'$と$z$上で線形補間をすることで評価し直し、消費関数をアップデートします。($g_1(a'_i,z_j)$とする。)\n",
"\n",
"ただし借入制約の存在には注意が必要です。借入制約の存在から下限については線形の外挿をするのではなく、ひと手間必要になります。\n",
"\n",
"\n",
"ステップ5で紹介したように$a^\\ast_{1,j}$は借入制約がバインドする時の今期の資産の量です。\n",
"\n",
"そこで$a'_i \\le a^\\ast_{1,j}$となるような$a'$のグリッドポイントについては次のように消費関数を求めます。\n",
"\\begin{equation*}\n",
"g_1(a'_i,z') = (1+r)a'_i + wz_j - a'_1\n",
"\\end{equation*}\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"EGM (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@inline function EGM(m::Model)\n",
"\n",
" #=便利な関数の設定=#\n",
" @inline u′(c, σ) = (1.0/c)^(σ) # 限界効用関数\n",
" @inline u′_inv(val, σ) = (1.0/val)^(1.0/σ) # 限界効用の逆関数\n",
" \n",
" # 要素価格\n",
" r = 0.03 # 利子率\n",
" w = 1.0 # 賃金\n",
" \n",
" # 簡素化のために設定\n",
" #a2mat = repeat(m.a2grid,1,m.nz)\n",
" #zmat = repeat(m.zgrid', m.na)\n",
"\n",
" #= 初期化 =#\n",
" a1= zeros(m.na,m.nz) # 内生的に決定される今期の資産のグリッド\n",
" wealth = (1.0+r)*m.a2grid .+ w*m.zgrid' #資産に関する政策関数のために必要\n",
" ctil =similar(a1) # 内生的に決まる消費\n",
"\n",
" pol_a = similar(a1) # 資産に関する消費関数\n",
" pol_c = ones(m.na,m.nz) # 消費に関する政策関数の初期予想\n",
" pol_c1 = similar(a1) # updateする消費に関する政策関数\n",
" \n",
" for iter in 1:m.maxiter\n",
" # オイラー方程式から内生的にcを求める\n",
" ctil .= u′_inv.((u′.(pol_c, m.σ)*m.β*(1.0+r)*m.Π'), m.σ) # step4\n",
" # 予算制約を利用して内生的にaのグリッドを作る\n",
" @. a1 = (m.a2grid + ctil - w*m.zgrid')/(1.0+r) # step 5\n",
" \n",
" # 政策関数の更新 \n",
" pol_c1 = policyinterp(m,r,w,a1,ctil,pol_c1) # step 6\n",
" err = norm(pol_c1-pol_c, Inf) # 誤差の計算\n",
" pol_c =copy(pol_c1) # 政策関数のアップデート\n",
"\n",
" if err < m.tol # step 7\n",
" pol_a = wealth - pol_c # 資産に関する政策関数を求める\n",
" break\n",
" end\n",
" \n",
" if iter == m.maxiter\n",
" println(\"The model does not converge...\")\n",
" end\n",
" end\n",
" \n",
" return pol_a, pol_c\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"policyinterp (generic function with 1 method)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@inline function policyinterp(m::Model, r0, w0, a1,\n",
" pol_c, pol_c_u)\n",
" \n",
" #政策関数を線形補間する\n",
" @inbounds for (i_z, v_z) in enumerate(m.zgrid)\n",
" @inbounds for (i_a, v_a) in enumerate(m.a2grid)\n",
" if v_a <= a1[1, i_z] # 借入制約にかかっている時はこのように消費を決定する\n",
" pol_c_u[i_a, i_z] = (1.0+r0)*v_a + w0*v_z - m.a2grid[1]\n",
"\n",
" elseif v_a >= a1[end, i_z] # 外挿は線形\n",
" d =(v_a-a1[end, i_z])*(pol_c[end, i_z] - pol_c[end-1, i_z])/(a1[end, i_z]-a1[end-1, i_z])\n",
" pol_c_u[i_a,i_z] = pol_c[end, i_z] + d\n",
" \n",
" else #定義域内でも線形補間\n",
" @views ind = searchsortedlast(a1[:,i_z],v_a)\n",
" ind1 = ind + 1\n",
" d = (v_a-a1[ind, i_z])*(pol_c[ind1, i_z] - pol_c[ind, i_z])/(a1[ind1, i_z]-a1[ind, i_z])\n",
" pol_c_u[i_a, i_z] = pol_c[ind, i_z] + d\n",
" end \n",
" end\n",
" end\n",
" return pol_c_u\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([0.0 0.448097; 0.0 0.448132; … ; 9.3614 10.0786; 9.40087 10.1182], [0.2 0.551903; 0.200041 0.55191; … ; 1.09736 1.18016; 1.09913 1.18184])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pol_a_EGM, pol_c_EGM = EGM(m)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"BenchmarkTools.Trial: \n",
" memory estimate: 19.78 MiB\n",
" allocs estimate: 2554\n",
" --------------\n",
" minimum time: 154.325 ms (0.00% GC)\n",
" median time: 159.805 ms (2.69% GC)\n",
" mean time: 168.083 ms (5.19% GC)\n",
" maximum time: 322.242 ms (38.72% GC)\n",
" --------------\n",
" samples: 30\n",
" evals/sample: 1"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@benchmark EGM(m)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"labels=[\"Bad State\",\"Good State\"]\n",
"plot(m.a2grid,pol_a_EGM, label=labels,legend=:bottomright)\n",
"plot!(m.a2grid,m.a2grid, label =\"45 degree\",legend=:bottomright)\n",
"title!(\"Policy function for asset using EGM\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"EGMの場合は約0.17秒で上のプロットのような政策関数が得られることが分かります。\n",
"これは先ほどのVFIに比べて**9倍**高速になっています。\n",
"\n",
"VFIで求められた政策関数とほぼ概形は同じで、資産が少ないうえに生産性が低い場合には貯蓄を行わないような状況になっています。\n",
"\n",
"EGMで得られた消費関数は線形補間を行っていることから滑らかになっています。\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 結論\n",
"\n",
"本記事ではベルマン方程式を解く際に使われる2つの手法について概説しました。 \n",
"\n",
"VFIは非常に安定的な解法ですが、モデルの構造をうまく利用するEGMを用いることで\n",
"同様の答えを**9倍**速く手に入れることができます。\n",
"\n",
"EGMのアイディアは理解するまでに時間がかかるかもしれませんが、コーディングは非常に簡単でおすすめです。\n",
"\n",
"ぜひみなさんもEGMを利用することで、人生を豊かにしてみてください。\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.2",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.2"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": true,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}