9.7.1 lqapproxについて

m21kosumi commented 7 years ago

lqapproxがどういったメソッドで成り立ち、どのように今回のモデルに適用できるのかがよくわかりません。オリジナルのMATLAB版に掲載されているlqapprox.mも参照しましたが、理解できない部分が複数あります。

まず、lqapprox.m内で使われているdpssというfunctionではsteady state valuesを計算しているようですが、dpss.m内の以下のコードの二行目はどういった操作を表しているのでしょうか。

theta=[s0(:);x0(:);lambda0(:)]; theta=real(broyden('dpssres',theta,func,params,e,delta,n,m)); sstar=theta(1:n)'; xstar=theta(n+1:n+m)'; lambdastar=theta(n+m+1:end)';
lqapprox内には関数fとgのderivativeを求める部分がありますが、これは教科書p.250のfunctionに記載されているf,gの３つのoutで代用できるのでしょうか。できない場合は、model functionに実装された３つのoutはどこで使用すればよいのでしょうか。
lqapprox内にはQZ decompositionとRiccati equation function iterationの二種類のmethodが実装されていますが、両者の違いは何でしょうか。

oyamad commented 7 years ago

よく読んでみないとちょっとよくわからないのですが，とりあえず lqapprox は無視して，value iteration を適当な initial value function から始めてみたらどうなりますか．

m21kosumi commented 7 years ago

そちらの方は、前の二人の発表者と同様にbellmanをupdateするメソッドで以下のように実装してあります。実行すると投資量xが極めて小さくなる問題が起きていますが…。

struct GrowthModel alpha::Float64 beta::Float64 gamma::Float64 sigma::Float64 delta::Float64 #discount factor s_vec::Vector{Float64} end

n = 10 smin = 5 smax = 10 fspace = Basis(ChebParams(n, smin, smax))

snodes = nodes(fspace)

GM = GrowthModel(0.2, 0.5, 0.9,0.1,0.9, snodes[1])

function update_bellman!(GM::GrowthModel,V::Vector) a,b,g,sigma,d = GM.alpha, GM.beta, GM.gamma, GM.sigma, GM.delta V_func = LinInterp(GM.s_vec, V) V_new = similar(V) x_opt = similar(V)

nshocks = 3
e,w = qnwlogn(nshocks,-(sigma^2)/2,sigma^2)

for (s_idx,s) in enumerate(GM.s_vec)
    objective(x) = -(((s - x)^(1 - a)/(1 - a)) + d * dot(V_func.(g*x + e*(x^b)), w))
    opt = optimize(objective, 1e-10, s)
    V_new[s_idx] = - opt.minimum
    x_opt[s_idx] = opt.minimizer
end
return V_new,x_opt
end

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

V_computed = similar(V) x_opt = similar(V)

for i in 1:100
    V_computed, x_opt = update_bellman!(GM,V)
end

これを実行すると、x_optは次のような値になります。

10-element Array{Float64,1}: 1.00002e-10 1.00002e-10 1.0e-10
1.00002e-10 1.00003e-10 1.00002e-10 1.0e-10
1.0e-10
1.0e-10
1.00002e-10

教科書中ではChebyshevとLQ approxinmationを結果を比較する図が載っているので、その図を出すためにはlqapproxを実装する必要があるように見えます。

oyamad commented 7 years ago

e にはどういう値が入っているのですか．

m21kosumi commented 7 years ago

eは [0.836771, 0.995012, 1.18318] になっています。

m21kosumi commented 7 years ago

update_bellman!はobjectiveでV_funcを使っていることに問題がありそうです。もう一度考え直してみます。

oyamad commented 7 years ago

メインのループ

for i in 1:100
    V_computed, x_opt = update_bellman!(GM,V)
end

で V を更新していないのが問題のように見えますが，copy!(V, V_computed) を加えるとどうなりますか．

m21kosumi commented 7 years ago

加えてみましたが、x_optにはあまり変化がありません。

10-element Array{Float64,1}: 1.00005e-10 1.00008e-10 1.00014e-10 1.00006e-10 1.00006e-10 1.00011e-10 1.00013e-10 1.00014e-10 1.00009e-10 1.00015e-10

oyamad commented 7 years ago

smin を小さい値にしてみてください．

m21kosumi commented 7 years ago

sminを0に設定したところ、次のようにx_optが出力されました！グラフを描画したところ、教科書のLQ approximationの図とほぼ一致しているように見えました。

10-element Array{Float64,1}: 0.0615535 0.538067 1.3794
2.41191
3.49467
4.50866
5.48781
6.35387
6.93494
7.25209

oyamad commented 7 years ago

smin=5 は大きすぎて，「任意の initial value function V に対してうまくいく」というわけにはいかない，ただし V が最初から答えに十分近ければうまくいく，ということのようですね．

oyamad commented 7 years ago

lqapprox の方は，定常状態の周りで効用関数は2次まで，生産関数は1次まで展開したもので置き換えることで強引に linear-quadratic のモデルだと思うことにして，あとは LQ Dynamic Programming Problems のやり方・関数で代用してみてください．

m21kosumi commented 7 years ago

今回の例をLQのモデルに当てはめるところで詰まってしまいました。 A,B,Cにあたる行列は(たぶん)生成できましたが、R,Qをどのように設定すればよいのかがわかりません。 value functionをx,sについてそれぞれ整理し、xstar,sstarとの偏差が出るようなquadraticな形に変形してR,Qを求める、ということでしょうか。 (それで解釈が正しければ、この式変形は私の手には負えそうにありません…)

oyamad commented 7 years ago

そう言わずにやってみましょう．

Utility flow の関数を u(c) として，それを定常状態の消費 cstar (= sstar - xstar) の周りで2次の多項式で置き換える： u(c) = u(cstar) + u'(cstar) (c - cstar) + (1/2) u''(cstar) (c - cstar)．
上式を平方完成して △ (c - cstar - ○)^2 + (定数項) の形にする．
c - cstar - ○ を control variable (lecture でいうところの u) とみなす．
定数項は最大化元に影響を与えないので無視してよい (value function の値はその分ずれるので注意)．
これで Q が定まるはず (R はゼロ)．
State transition の方に定数項が出てきたら，state variable の定義を工夫して調整する (lecture 参照)．

これでやってみてください．

m21kosumi commented 7 years ago

なんとかLQモデルらしきものを組み立てましたが、いずれかの(もしくは全ての)設定が間違っているようで、errorが出て動きません。

GM::GrowthModel a,b,gamma,sigma,d = GM.alpha, GM.beta, GM.gamma, GM.sigma, GM.delta

nshocks = 3
e,w = qnwlogn(nshocks,-(sigma^2)/2,sigma^2)

 A = [b*dot(e*(xstar),w) gamma*xstar+(1+b)*dot(e*(xstar^b),w)+gamma;
     　0                 1]
 B = [-(b*dot(e*(xstar),w));
     　0]
 C = zeros(1,2)

 R = zeros(2,2)
 Q = [-1/2*a*((sstar -xstar)^(-a-1))]

 lq = LQ(Q,R,A,B,C; bet=b)

 x0 = [5.0;1.0]

 compute_sequence(lq, x0)

以下がerror全文です。

MethodError: no method matching getZ(::Array{Float64,1}, ::Float64, ::Float64) Closest candidates are: getZ(::Float64, ::Float64, ::Float64) at C:\Users\mao21.julia\v0.6\QuantEcon\src\matrix_eqn.jl:192 getZ(::Float64, ::Float64, ::Union{Array{T,1} where T, Array{T,2} where T}) at C:\Users\mao21.julia\v0.6\QuantEcon\src\matrix_eqn.jl:207 getZ(::Array{T,2} where T, ::Float64, ::Array{T,2} where T) at C:\Users\mao21.julia\v0.6\QuantEcon\src\matrix_eqn.jl:222

 Stacktrace:
  [1] #solve_discrete_riccati#105(::Float64, ::Int64, ::Function, ::Array{Float64,2}, ::Array{Float64,1}, 
 ::Array{Float64,2}, ::Array{Float64,1}, ::RowVector{Float64,Array{Float64,1}}) at 
 C:\Users\mao21\.julia\v0.6\QuantEcon\src\matrix_eqn.jl:123
  [2] stationary_values!(::QuantEcon.LQ) at C:\Users\mao21\.julia\v0.6\QuantEcon\src\lqcontrol.jl:202
  [3] compute_sequence(::QuantEcon.LQ, ::Array{Float64,1}, ::Int64) at 
 C:\Users\mao21\.julia\v0.6\QuantEcon\src\lqcontrol.jl:306
  [4] compute_sequence(::QuantEcon.LQ, ::Array{Float64,1}) at 
 C:\Users\mao21\.julia\v0.6\QuantEcon\src\lqcontrol.jl:305
  [5] include_string(::String, ::String) at .\loading.jl:515

oyamad commented 7 years ago

入力データの次元がどれか合っていないということなのでしょうね．

ちなみに，state transition の1次式での近似は，9.1節によると epsilon は epsilon_star で固定してランダムネスはないものとしてやっているようです (したがって qnwlogn とか使わない)．

LQ での state と control は何にしましたか．

m21kosumi commented 7 years ago

stateはs_{t+1} = gamma(s_t-c_t)+estar(s_t-c_t)^bで、controlは(c - cstar + cstar/alpha)になっています。

oyamad commented 7 years ago

B と C の次元 (size) は同じでないといけないと思います．
Q をスカラーにしてみてください．
bet=d でしょう．

oyamad commented 7 years ago

間違いがわかりました．Control variable を c (+ 調整項) にしたのがまずかった．LQ では制約条件も無視するので，stock をぜんぶ消費するのが最適になってしまうに決まってますね．x (+ 調整項) を control variable にしないといけない．

m21kosumi commented 7 years ago

u(s,x)を(sstar,xstar)について展開・平方完成してstateとcontrol variableを導出したあと、それに合わせてstate transitionを調整する、という理解で合っていますか？

oyamad commented 7 years ago

合っています．

m21kosumi commented 7 years ago

u(s,x)を展開してみたのですが、この際に出てくるsとxの積の項はどのように処理したらよいでしょうか…？

oyamad commented 7 years ago

Cross term の係数行列は optional argument N で加えられます (documentation)．

m21kosumi commented 7 years ago

次のように設定して動かしてみましたが、正しく動いているように見えません…。 conpute_sequenceも動かしてみましたが、sを全て投資に回してしまっているようです。

GM::GrowthModel a,b,gamma,sigma,d = GM.alpha, GM.beta, GM.gamma, GM.sigma, GM.delta

cstar = sstar - xstar 
m1 = gamma+b*(xstar)^(b-1)
m2 = 1/2*a*(cstar^(-a-1))

A = [0 gamma*xstar+(xstar^b)-sstar+(cstar/a)*(m1 - 1);
     0                 1]
B = [m1;
     0]
C = zeros(2,1)

R = [m2 0;
     0  1]
Q = m2
N = [-m2 0]

lq = LQ(Q,R,A,B,C,N; bet=d)

P, F, d = stationary_values(lq)

([0.0 0.0; 0.0 10.0], [-1.0 0.0], 0.0)

s0 = [5.0;1.0] slq,xlq = compute_sequence(lq,x0)

 ([5.0 6.55971 … 4.75631e5 5.2848e5; 1.0 1.0 … 1.0 1.0], [5.0 6.55971 … 4.28067e5 4.75631e5], [0.0 0.0 … 
 0.0 0.0; 0.0 0.0 … 0.0 0.0])

m21kosumi commented 7 years ago

sとxの積の部分にstateとcontrol variableの定数部分を加味するのを忘れていたのが原因のようです。もう一度考え直してみます。

m21kosumi commented 7 years ago

u(s,x)を a(s-s-b)^2 + c(x-x-d)^2 + 2e(s-s-b)(x-x-d) と置いて展開し、a,b,c,d,eを求めることでstateとcontrol variableを導出しようとしたのですが、a,c,eを求めsとxの係数について連立方程式を立てたところ、bとdの部分が消えてしまい解が求まりません。なにか他に式変形の方法はありますか…？

oyamad commented 7 years ago

それを展開して2次の Taylor 展開の式 (例えばp.226の最初の式の右辺) と係数比較するので合ってます．(1次の項 f_s^* (s - s^*), f_x^* (x - x^*) を忘れてませんか．)

m21kosumi commented 7 years ago

一次の項も計算には含まれています。もう一度始めから計算しなおしてみましたが、やはりsの係数についての等式とxの係数についての等式が、どちらもcstar = alpha*(b-d)になってしまいます。

oyamad commented 7 years ago

2a = f_ss, 2c = f_xx, 2e = f_sx で，y = (-b, -d) とすると Hf y = Df (ただし，Df は f の (s, x) での gradient vector，Hf は f の (s, x) での Hessian matrix) となりませんか．

m21kosumi commented 7 years ago

そうなりますが、a=cのためHessian matrixが実質的に[2a 2e; 2e 2a]に、またf_s = -f_xのためgradient vectorも実質的に(f_s^* -f_s^*)になり、yについて解くことができません。

oyamad commented 7 years ago

ああなるほど．このアプローチ (s-s, x-x をさらに s-s-b, x-x-d とずらす) が根本的によくなかったですね．すみません．

State variable は (s, 1) と次元を増やすことになるので，それを使うべきでしたね． (s 1) R (s 1)' + Q x^2 + (s 1) N x (R は 2x2, Q はスカラー，N は 2x1) として (ただし，s-s, x-x を見やすくするため単に s, x と書いている)，2次の展開式と係数比較すればよさそうです． R = [(1/2)*f_ss (1/2)*f_s; (1/2)*f_s f], Q = (1/2)*f_xx, N = [f_sx; f_x] (ただし，それぞれ s-s, x-x で評価する) という感じになりませんか (係数は，ずれているかもしれない)．

m21kosumi commented 7 years ago

前回のゼミを踏まえてLQモデルを作り直しているのですが、A,B,Cはどう設定すればよいのでしょうか。 (というより、そもそもs_tをs_t+1とみてg(x)を代入する、という理解で合っていますか…？)

GM::GrowthModel a,b,gamma,sigma,d = GM.alpha, GM.beta, GM.gamma, GM.sigma, GM.delta

estar = 1 
xstar = ((1 - d*gamma)/(d*b))^(1/(b - 1))
sstar = gamma*xstar + xstar^b
cstar = sstar - xstar

u = ((gamma-1)*xstar + xstar^b)/(1-a)
ux = (((gamma-1)*xstar + xstar^b)^(-a))*(gamma-1+b*(xstar^(b-1)))
uxx = ((-a*(gamma-1+b*(xstar^(b-1))))^2)*((gamma-1)*xstar+xstar^b)^(-a-1)
      +(((gamma-1)*xstar+xstar^b)^(-a))*b*(b-1)*(xstar^(b-2))

A = [1 0;
    0 0]
B = zeros(2,1)
C = zeros(2,1)

R = [0 0;
     0 -u]
Q = -1/2*uxx
N = [0 -1/2*ux]

lq = LQ(Q,R,A,B,C,N; bet=d)

oyamad commented 7 years ago

一つのやり方は

g(s, x) = x
f(s, x) = u((gamma * s + s^beta) - x)

とするものです．

oyamad commented 6 years ago

QuantEcon.jl の stationary_values(::LQ) はうまくいかないので，この notebook の Value Function Iteration の section のようにやってみてください．

上の g(s, x) = x のやり方でも，いままでの g(s, x) = gamma * (s - x) + (s - x)^beta のやり方でも，どちらでもよいです．

oyamad commented 6 years ago

@m21kosumi @IoriS QuantEcon.jl の solve_discrete_riccati の bug fix をしたので，それをインストールして stationary_values(::LQ) を試してみてください．~~まだマージされてなくて fix-riccati ブランチにあります．~~

EDIT: マージされたので

Pkg.checkout("QuantEcon", "master")

とするだけでOK．

(もとのリリース版に戻すときは

Pkg.free("QuantEcon")

とする．)

m21kosumi commented 6 years ago

lqapproxに用いるstationary_values自体に問題があったことが原因でした。

OyamaZemi / MirandaFackler.notebooks

9.7.1 lqapproxについて #5