次のような動学的最適化問題を考える．
消費からの効用の現在価値の総和を最大にするような消費経路は何か．
数式で書くと： $$ \begin{align*} &\max_{\{c_t\}_{t=0}^{\infty}} \sum_{t=0}^{\infty} \beta^t u(c_t) \\\\ &\ \text{ s.t. }\ k_{t+1} = f(k_t) - c_t,\quad \text{$k_0$: given} \end{align*} $$
目的関数の最大値を初期条件 $k_0$ の関数として $v^*(k_0)$ と書く
(最適解の存在を保証する条件が満たされていると暗黙に仮定)．
この関数 $v^*$ を最適価値関数 (optimal value function)，あるいは単に価値関数という．
In [1]:
using QuantEcon, QuantEcon.Models
using PyPlot
具体的な関数形として
を用いる．
(実は，このケースは手で解ける．)
また $\beta = 0.95$ とする．
In [2]:
# Production function
a = 0.65
f(k) = k^a
# Instantaneous utility function
u(c) = log(c)
# Discount factor
b = 0.95
Out[2]:
モデル上は状態 (資本量) $k$ は連続変数だが，数値計算のためには離散化しないといけない．
In [3]:
# Grid
# grid_min = 1e-6
grid_max = 2
grid_size = 150
Out[3]:
QuantEcon.Models の GrowthModel タイプのインスタンスを作る：
In [4]:
gm = GrowthModel(f, b, u, grid_max, grid_size)
Out[4]:
bellman_operator メソッドは (その名の通り) Bellman operator を生成する：
(モデル上は $v$ は連続変数についての関数だが，ここでは長さ grid_size の配列．
内部で線形補間を行っている．)
In [5]:
T(v) = bellman_operator(gm, v)
Out[5]:
Value iteration によって価値関数 $v^*$ を求める．
compute_fixed_point は初期関数 v_0 に operator T を何度も (最大 max_iter 回) 施す．
変化分が err_tol (デフォルト値 1e-3) 以下になったら止まる．
In [6]:
v_0 = zeros(grid_size)
v_star = compute_fixed_point(T, v_0, max_iter=200)
Out[6]:
価値関数 $v^*$ の様子：
In [7]:
fig, ax = subplots()
ax[:set_title]("Optimal value function")
ax[:set_xlabel]("\$k\$")
ax[:set_ylabel]("\$v^*(k)\$")
ax[:set_ylim](-40, -30)
ax[:plot](gm.grid, v_star)
Out[7]:
最適政策 $\sigma$ を求める．
bellman_operator は ret_policy=true とすると max の代わりに argmax を返す．
In [8]:
sigma = bellman_operator(gm, v_star, ret_policy=true)
Out[8]:
In [9]:
fig, ax = subplots()
ax[:set_title]("Optimal policy function")
ax[:set_xlabel]("\$k\$")
ax[:set_ylabel]("\$\\sigma(k)\$")
ax[:plot](gm.grid, sigma)
Out[9]:
最後に，最適政策 $\sigma$ の下での資本量の推移を見る．
(sigma は長さ grid_size の配列．
補間しないといけない．)
In [10]:
using Grid: CoordInterpGrid, BCnan, InterpLinear
sigma_interpolated =
CoordInterpGrid(gm.grid, sigma, BCnan, InterpLinear)
series_length = 25
k = Array(Float64, series_length)
k[1] = 0.1
for t=1:series_length-1
k[t+1] = gm.f(k[t]) - sigma_interpolated[k[t]]
end
In [11]:
fig, ax = subplots()
ax[:set_title]("Capital accumulation")
ax[:set_xlabel]("\$t\$")
ax[:set_ylabel]("\$k_t\$")
ax[:set_ylim](0.1, 0.3)
ax[:plot](0:series_length-1, k, "o-")
Out[11]:
In [ ]: