9.7.2 - Githubissues

IoriS commented 7 years ago

小住さんと同様の方法で行ったのですが、policyが異常に小さくなる上にVの値が負になってしまいました。sminを0にしても解決しなかったのですがどうしたらよろしいでしょうか…

using QuantEcon
using BasisMatrices
using Plots
using Optim

struct RenewableResourceManagement
    alpha::Float64
    beta::Float64
    gamma::Float64
    kappa::Float64
    delta::Float64
    s_vec::Vector{Float64}
end

n = 8
smin = 6
smax = 9
fspace = Basis(ChebParams(n,smin,smax))
snodes = nodes(fspace)

RRM = RenewableResourceManagement(4.0,1.0,0.5,0.2,0.9,snodes[1])

function update_bellman!(RRM::RenewableResourceManagement,V::Vector)
    a,b,g,k,d = RRM.alpha, RRM.beta, RRM.gamma, RRM.kappa, RRM.delta
    V_func = LinInterp(RRM.s_vec, V)
    V_new = similar(V)
    policy = similar(V)

    for (i, s) in enumerate(RRM.s_vec)
        objective(x) =  (x^(1-g))/(1-g) -k*x +d*(V_func.(a*(s-x)-0.5*b*(s-x)^2))
        opt = optimize(objective,1e-10,s)
        V_new[i] = -opt.minimum
        policy[i] = opt.minimizer
    end
    return V_new,policy
end

V = Vector{Float64}(length(RRM.s_vec))
for i in 1:length(RRM.s_vec)
    V[i] = 6
end

V_comp = similar(V)
policy = similar(V)
tol = sqrt(eps())
max_iter = 1000
V_error = 1.0
i = 1
while V_error > tol && i <= max_iter
    V_comp ,policy = update_bellman!(RRM,V)
    V_error = maximum(abs, V_comp - V)
    copy!(V,V_comp)
    i += 1
end

oyamad commented 7 years ago

optimize は最小化なのに最大化したいものを与えているということはないですか．

IoriS commented 7 years ago

objectiveの関数を-にするのを忘れていただけでした・・・ありがとうございます

IoriS commented 7 years ago

LQモデルに当てはめる方なんですが、transition functionのs_{t+1} = alpha (s_t-x_t)-0.5beta(s_t-xt)^2をlectureのx{t+1}=Ax_t+Bu_tの形に変形する部分が上手くできないのですがどうしたらよいでしょうか・・・

oyamad commented 7 years ago

Transition function の方は1次式で置き換えるということなので，単に2次の項を無視してはどうですか．

IoriS commented 7 years ago

alpha(s_t-x_t)の部分だけ取り出すということでしょうか

oyamad commented 7 years ago

そうです．

IoriS commented 7 years ago

この後u=x-xstarのようにしてlectureのAとBを出してR=0,Q=1としてlectureの例と同じようにRfを設定するので大丈夫ですか…？

oyamad commented 7 years ago

まずはコードを書いて動かしてみる． (https://github.com/OyamaZemi/MirandaFackler.notebooks/issues/5#issuecomment-333311225 のように，control variable をもう少しずらさないといけない気がするが．)

IoriS commented 7 years ago

u(x)をxstarの周りで2次までで展開したものを(x-z)^2+~の形にしてlectureのuをu=x-zとしてやったのですが明らかにおかしい結果が出てしまいました・・・

a,b,g,k,d = RRM.alpha, RRM.beta, RRM.gamma, RRM.kappa, RRM.delta
sstar = (a^2 - 1/(d^2))/(2*b)
xstar = sstar - (d*a - 1)/(d*b)
ax = (xstar^(g+1))*(xstar^(-g)-k+0.5*g*xstar^(-g))/g

Q = 1.0
R = zeros(2,2)
A = [a -a*ax
        0.0 1.0]
B = [-a; 0]
lq = LQ(Q,R,A,B;bet=d)

s0 = [0.0;1.0]

sp, xp = compute_sequence(lq,s0)

stock = vec(sp[1,:])

101-element Array{Float64,1}: 0.0
-29.6962
-148.481
-623.62
-2524.18
-10126.4
-40535.3
-1.62171e5 -6.48713e5 -2.59488e6 -1.03796e7 -4.15183e7 -1.66073e8 ⋮
-3.79244e54 -1.51698e55 -6.0679e55 -2.42716e56 -9.70864e56 -3.88346e57 -1.55338e58 -6.21353e58 -2.48541e59 -9.94165e59 -3.97666e60 -1.59066e61

oyamad commented 7 years ago

alpha(s_t-x_t)の部分だけ取り出すということでしょうか

すみません，これは「そうです」ではありませんでしたね．s - x を (s - sstar) - (x - xstar) + (sstar - xstar) と書き換えて，2次式を展開して整理して，2次の項 ((s - sstar)^2, (x - xstar)^2, (s - sstar) (x - xstar)) を無視する，ということでした．

IoriS commented 7 years ago

shadow priceをプロットするのにBasisMatricesのfunevalを使う必要がありそうなのですが、CompEconの方だとfuneval(c, fspace, x, order)でorderを1にすると補間の関数を1回微分してくれるそうなのですがBasisMatrices版だとorderは0以外受け付けない旨のエラーが出てしまうのですが何か代わりの方法ありませんでしょうか（BasisMatricesの元のコード見に行ってもよくわかりませんでした…） ''' veval = funeval(c, fspace, RRM.s_vec,1)

passing order as integer only allowed for 0. Try calling the version where order is a matrix '''

ghost commented 7 years ago

まだちゃんと読んではいませんが、evalbase でスプラインなら、基底関数の導関数を作れる気がします

ghost commented 7 years ago

ちょっと試したけどうまくいなかったので、evalbase で本当に微分できるか自信ないです…

BasisMarices がだめそうなら、微分してくれるパッケージ(使ったことないけど ForwardDiff.jl とか？)を使ってみるとか？

oyamad commented 7 years ago

Documentation を読むと evalbase でできそうだけど，うまくいかなかったというのはどううまくいかなかったの？

ghost commented 7 years ago

evalbase を使ったら、基底関数っぽい行列が出てきたので、それを Φ.vals[1] の代わりに用いて、あとは微分しない場合と全く同様に計算してみました
そうして、funeval で出てきた値をプロットしてみたところ、適当に簡単な関数をとっていたのですが、その導関数とは異なるグラフを描きました

oyamad commented 7 years ago

具体的にどういうコードを実行した？

ghost commented 7 years ago

basis = Basis(SplineParams(5, -2, 2, 3))
xs,_ = nodes(basis)
a = evalbase(SplineParams(5, -2, 2, 3), xs, 1)
ys = xs .^2
cd = a \ ys
plot(collect(-2:0.1:2), funeval(cd, basis, collect(-2:0.1:2)))

と真の関数を y=x^2 にしたのですが、その導関数がうまくプロットされませんでした

ys の値が元の関数のままなのがいけないのかと思って、

ys2 = xs .*2
cd2 = a \ ys2
plot(collect(-2:0.1:2), funeval(cd2, basis, collect(-2:0.1:2)))

としても、思ったような図になりませんでした

oyamad commented 7 years ago

これでどうですか：

f(x) = x^3
a, b = -2, 2
x = linspace(a, b, 50)

basis = Basis(SplineParams(5, a, b, 3))
S, (xgrid,) = nodes(basis)
Phi = BasisMatrix(basis, Expanded(), S, 0)
y = f.(xgrid)
c = Phi.vals[1] \ y
interp0 = funeval(c, basis, x)
plot(x, interp0)

order = 1
B1 = evalbase(basis.params[1], x, order)
interp1 = B1 * c
plot(x, interp1)

order = 2
B2 = evalbase(basis.params[1], x, order)
interp2 = B2 * c
plot(x, interp2)

ghost commented 7 years ago

ありがとうございます。正しく微分できているようです

係数は微分しないときと共通で、それに evalbase で作った基底関数の行列を乗じれば、補間関数の導関数を評価することができるのですね

IoriS commented 7 years ago

自分もevalbaseで目的の関数が得られました、ありがとうございます

LQの方なんですが色々考えたのですがやはりu(x)に0<x<sの制約を付ける方法がわかりません…

oyamad commented 7 years ago

LQ は制約なし (無視) です．

IoriS commented 7 years ago

i=s-xとしてu(x)をiとsについて変形してこれを最小化するように設定したのですがやはりうまくいきませんでした…

a,b,g,k,d = RRM.alpha, RRM.beta, RRM.gamma, RRM.kappa, RRM.delta
sstar = (a^2 - 1/(d^2))/(2*b)
xstar = sstar - (d*a - 1)/(d*b)
istar = xstar - sstar
c1 = xstar^(-g)-k
c2 = 0.5*(-g*xstar^(-g-1))
c3 = (xstar^(1-g))/(1-g) -(k+c1)*xstar +c2*xstar^2
c4 = (c1-2*c2*xstar)/(2*c2)

Q = c2
R = [c2 0
        0 0]
N = [-c2 0]
A = [0 c4*(a-b*istar)+0.5*b*istar^2
        0.0 1.0]
B = [a-b*istar; 0]
C = zeros(2,2)
lq = LQ(Q,R,A,B,C,N;bet=d)

s0 = [0.0;1.0]
sp, up = compute_sequence(lq,s0)
stock = vec(sp[1,:])

101-element Array{Float64,1}:
      0.0       
    -62.4493    
   -492.656     
  -3456.3       
 -23872.5       
     -1.64518e5 
     -1.13341e6 
     -7.80797e6 
     -5.37883e7 
     -3.70542e8 
     -2.55262e9 
     -1.75847e10
     -1.21139e11
      ⋮         
     -4.1761e75 
     -2.87687e76
     -1.98185e77
     -1.36527e78
     -9.4052e78 
     -6.47914e79
     -4.46341e80
     -3.07479e81
     -2.11819e82
     -1.4592e83 
     -1.00522e84
     -6.92488e84

IoriS commented 7 years ago

先ほど伺った方法を試してみたのですがやはり変な結果になりました…（一応notebook上げました）

IoriS commented 7 years ago

今計算したらRがnonnegative definiteでなければならないという条件をクリアしてなかったので一般にこの式変形でLQに落とし込めるというわけではないのかもしれません・・・

oyamad commented 7 years ago

実装上どうなっているかわからないけど，数学的には f_ss が非正なら大丈夫なはずです．

IoriS commented 7 years ago

f(s,x)=log(x^a-1),g(s,x)=xの例でstateとcontrolをずらさない方法を試したのをノートブックの末尾に追加しましたが、やはり上手くいっていないようです…

oyamad commented 6 years ago

A の方も変更しないといけないのを忘れてませんか．

OyamaZemi / MirandaFackler.notebooks

9.7.2 #7