設定は7.6.2 Asset Replacementを引き継いでいる。
毎年初めにkeepかreplaceか選択する。a年経過後のoutputがq(a)、replacement costがk
kが$k_{t+1} = \bar{k} + γ(k_t- \bar{k}) + ε$で決定される
state variables: k, a
action variables: x (1:keep, 2:replace)
Reward functionは、 $f(p, a, x) = pq(a)$ ($x = keep$)、$f(p, a, x) = pq(0)-c$ ($x = replace$)
Value functionは、 $V(p, a) = \max \{pq(a) + δE_ε V(h(p, ε), a+1), pq(0) - c + E_ε V(h(p, ε), 1) \}$
In [1]:
using QuantEcon
using Plots
In [2]:
function q(a)
return 50 - 2.5 * a - 2.5 * a^2
end
Out[2]:
In [3]:
struct AssetReplacement
price::Float64
kbar::Float64
gamma::Float64
abar::Int
sigma::Float64
delta::Float64
k_vec::Vector{Float64}
function AssetReplacement(price=2,
kbar=100,
gamma=0.5,
abar=6,
sigma=15,
delta=0.9,
k_vec=collect(linspace(30, 190, 100)))
return new(price, kbar, gamma, abar, sigma, delta, k_vec)
end
end
In [4]:
function update_bellman!(ar::AssetReplacement, Vs, Vs_new)
price, kbar, gamma, abar, sigma, delta = ar.price, ar.kbar, ar.gamma, ar.abar, ar.sigma, ar.delta
V_funcs = [LinInterp(ar.k_vec, Vs[:, a]) for a in 1:abar]
e, w = qnwnorm(5, 0, sigma)
policy = similar(Vs, Int)
for a in 1: abar-1
for (k_idx, k) in enumerate(ar.k_vec)
V_keep = price * q(a) + delta * dot(V_funcs[a+1].(e .+ kbar .+ (gamma * (k-kbar))), w)
V_replace = price * q(0) - k + delta * dot(V_funcs[1].(e .+kbar .+ (gamma * (k-kbar))), w)
Vs_new[k_idx, a] = max(V_keep, V_replace)
if V_keep > V_replace
policy[k_idx, a] = 1
else
policy[k_idx, a] = 2
end
end
end
for (k_idx, k) in enumerate(ar.k_vec)
Vs_new[k_idx, abar] = price * q(0) - k + delta * dot(V_funcs[1].(e .+kbar .+ (gamma * (k-kbar))), w)
policy[k_idx, abar] = 2
end
return Vs_new, policy
end
Out[4]:
In [5]:
ar = AssetReplacement()
Out[5]:
In [6]:
Vs = ones(length(ar.k_vec), ar.abar);
In [7]:
Vs_computed = similar(Vs)
tol=sqrt(eps())
max_iter=500
error = tol + 1
i = 0
resid = similar(Vs)
while error > tol && i < max_iter
update_bellman!(ar, Vs, Vs_computed)
copy!(resid, Vs-Vs_computed)
error = maximum(abs, resid)
copy!(Vs, Vs_computed)
i += 1
end
In [8]:
_,policy = update_bellman!(ar, Vs, similar(Vs))
Out[8]:
In [9]:
policy
Out[9]:
In [10]:
Vs
Out[10]:
In [11]:
plotlyjs()
Out[11]:
In [12]:
plot(ar.k_vec, [Vs[:, a] for a in 1:ar.abar-1], xlims=(0, 200), label=["age = 1" "age = 2" "age = 3" "age = 4" "age = 5"])
Out[12]:
In [13]:
plot(ar.k_vec, [policy[:, a] for a in 1:ar.abar-1], xlims=(0, 200), label=["age = 1" "age = 2" "age = 3" "age = 4" "age = 5"])
Out[13]: