Commit 250d204c authored by Michel Juillard's avatar Michel Juillard
Browse files

adding updated a, Pinf and Pstar to Diffuse Filter

parent 9bd9e1cd
......@@ -223,6 +223,27 @@ function get_QQ!(c::AbstractMatrix{T}, a::AbstractMatrix{T}, b::AbstractMatrix{T
mul!(c, work, transpose(a))
end
# att = a + K'*v
function get_updated_a!(att::AbstractVector{T}, a::AbstractVector{T}, K::AbstractMatrix{T}, v::AbstractVector{T}) where T <: AbstractFloat
copy!(att, a)
mul!(att, transpose(K), v, 1.0, 1.0)
end
# Ptt = P - K'*Z*P
function get_updated_Ptt!(Ptt::AbstractMatrix{T}, P::AbstractMatrix{T}, K::AbstractMatrix{T}, ZP::AbstractMatrix{T}) where T <: AbstractFloat
copy!(Ptt, P)
mul!(Ptt, transpose(K), ZP, -1.0, 1.0)
end
# Pstartt = Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar %(5.14) DK(2012)
function get_updated_Pstartt!(Pstartt::AbstractMatrix{T}, Pstar::AbstractMatrix{T}, ZPstar::AbstractMatrix{T},
Kinf::AbstractMatrix{T}, ZPinf::AbstractMatrix{T},
Kstar::AbstractMatrix{T}, vPinftt::AbstractMatrix{T}, PTmp::AbstractMatrix{T}) where T <: AbstractFloat
copy!(Pstartt, Pstar)
mul!(Pstartt, transpose(ZPstar), Kinf, -1.0, 1.0)
mul!(Pstartt, transpose(ZPinf), Kstar, -1.0, 1.0)
end
# v = y - Z*a -- basic
function get_v!(v::AbstractVector{T}, y::AbstractVecOrMat{T}, z::AbstractVecOrMat{T}, a::AbstractVector{T}, iy::U, ny::U) where {T <: AbstractFloat, U <: Integer}
copyto!(v, 1, y, iy, ny)
......@@ -381,6 +402,12 @@ function update_a!(a1::AbstractArray{U}, a::AbstractArray, d::AbstractArray{U},
a1 .+= d
end
# a = d + T*att
function update_a!(a1::AbstractVector{U}, d::AbstractVector{U}, T::AbstractMatrix{U}, att::AbstractVector{U}) where U <: AbstractFloat
copy!(a1, d)
mul!(a1, T, att)
end
# a = d + a + K'*v
function filtered_a!(a1::AbstractArray{U}, a::AbstractArray{U}, d::AbstractArray{U}, K::AbstractArray{U}, v::AbstractArray{U}, work::AbstractArray{U}) where U <: AbstractFloat
a1 .= a
......@@ -440,29 +467,35 @@ function filtered_P!(P1::AbstractArray{U}, P::AbstractArray{U}, K::AbstractArray
gemm!('T', 'N', -1.0, K, ZP, 1.0, P1)
end
# P = T*P*T'+ QQ
function update_P!(P::AbstractArray{U}, P1::AbstractArray{U}, T::AbstractArray{U}, QQ::AbstractArray{U}, Ptmp::AbstractArray{U}) where U <: AbstractFloat
mul!(Ptmp, T, P1)
# P = T*Ptt*T' + QQ
function update_P!(P::AbstractMatrix{U}, T::AbstractMatrix{U}, Ptt::AbstractMatrix{U}, QQ::AbstractMatrix{U}, Ptmp::AbstractMatrix{U}) where U <: AbstractFloat
mul!(Ptmp, Ptt, transpose(T))
copy!(P, QQ)
gemm!('N', 'T', 1.0, Ptmp, T, 1.0, P)
mul!(P, T, Ptmp, 1.0, 1.0)
end
# Pinf = T*Pinftt*T'
function update_P!(P::AbstractMatrix{U}, T::AbstractMatrix{U}, Ptt::AbstractMatrix{U}, Ptmp::AbstractMatrix{U}) where U <: AbstractFloat
mul!(Ptmp, Ptt, transpose(T))
mul!(P, T, Ptmp)
end
# Pstar = T*(Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar)*T'+QQ; %(5.14) DK(2012)
function update_Pstar!(Pstar1, Pstar, T, ZPinf, ZPstar, Kinf, Kstar, QQ, PTmp)
function update_Pstar!(Pstar, T, ZPinf, ZPstar, Kinf, Kstar, QQ, PTmp)
copy!(PTmp, Pstar)
mul!(PTmp, transpose(ZPstar), Kinf, -1.0, 1.0)
mul!(PTmp, transpose(ZPinf), Kstar, -1.0, 1.0)
copy!(Pstar, PTmp)
mul!(PTmp, T, Pstar)
copy!(Pstar1, QQ)
mul!(Pstar1, PTmp, transpose(T), 1.0, 1.0)
copy!(Pstar, QQ)
mul!(Pstar, PTmp, transpose(T), 1.0, 1.0)
end
# Pinf = T*(Pinf-Pinf*Z'*Kinf)*T'; %(5.14) DK(2012)
function update_Pinf!(Pinf1, Pinf, T, ZPinf, Kinf, PTmp)
function update_Pinf!(Pinf, T, ZPinf, Kinf, PTmp)
mul!(Pinf, transpose(ZPinf), Kinf, -1.0, 1.0)
mul!(PTmp, T, Pinf)
mul!(Pinf1, PTmp, transpose(T))
mul!(Pinf, PTmp, transpose(T))
end
# r_{t-1} = Z_t'*iF_t*v_t + L_t'r_t
......
......@@ -9,7 +9,9 @@ function kalman_filter!(Y::AbstractArray{U},
R::AbstractArray{U},
Q::AbstractArray{U},
a::AbstractArray{U},
att::AbstractArray{U},
P::AbstractArray{U},
Ptt::AbstractArray{U},
start::V,
last::V,
presample::V,
......@@ -51,9 +53,11 @@ function kalman_filter!(Y::AbstractArray{U},
vR = changeR ? view(R, :, :, t) : view(R, :, :)
vQ = changeQ ? view(Q, :, :, t) : view(Q, :, :)
va = changeA ? view(a, :, t) : view(a, :)
vatt = changeA ? view(att, :, t) : view(att, :)
va1 = changeA ? view(a, :, t + 1) : view(a, :)
vd = changeD ? view(d, :, t) : view(d, :)
vP = changeP ? view(P, :, :, t) : view(P, :, :)
vPtt = changeP ? view(Ptt, :, :, t) : view(Ptt, :, :)
vP1 = changeP ? view(P, :, :, t + 1) : view(P, :, :)
if changeR || changeQ
get_QQ!(ws.QQ, vR, vQ, ws.RQ)
......@@ -91,18 +95,24 @@ function kalman_filter!(Y::AbstractArray{U},
# K = iF*Z*P
get_K!(vK, vZP, vcholF)
end
# a = d + T(a + K'*v)
update_a!(va1, va, vd, vK, vv, ws.a1, vT)
# att = a + K'*v
get_updated_a!(vatt, va, vK, vv)
# a = d + T*att
update_a!(va1, vd, vT, vatt)
if !steady
# P = T*(P - K'*Z*P)*T'+ QQ
copy!(vP1, vP)
update_P!(vP1, vT, ws.QQ, vK, vZP, ws.PTmp)
copy!(ws.oldP, vP)
# Ptt = P - K'*Z*P
get_updated_Ptt!(vPtt, vP, vK, vZP)
# P = T*Ptt*T + QQ
update_P!(vP1, vT, vPtt, ws.QQ, ws.PTmp)
ws.oldP .-= vP1
if norm(ws.oldP) < ns*eps()
steady = true
else
copy!(ws.oldP, vP)
end
elseif t > 1
copy!(vP1, vP)
vPtt1 = view(ws.Ptt, :, : , t-1)
copy!(vPtt, vPtt1)
end
end
t += 1
......@@ -191,8 +201,11 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
R::AbstractArray{U},
Q::AbstractArray{U},
a::AbstractArray{U},
att::AbstractArray{U},
Pinf::AbstractArray{U},
Pinftt::AbstractArray{U},
Pstar::AbstractArray{U},
Pstartt::AbstractArray{U},
start::V,
last::V,
presample::V,
......@@ -235,11 +248,14 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
get_QQ!(ws.QQ, vR, vQ, ws.RQ)
end
va = changeA ? view(a, :, t) : view(a, :)
vatt = changeA ? view(att, :, t) : view(att, :)
va1 = changeA ? view(a, :, t + 1) : view(a, :)
vd = changeD ? view(d, :, t) : view(d, :)
vPinf = changePinf ? view(Pinf, :, :, t) : view(Pinf, :, :)
vPinftt = changePinf ? view(Pinftt, :, :, t) : view(Pinftt, :, :)
vPinf1 = changePinf ? view(Pinf, :, :, t + 1) : view(Pinf, :, :)
vPstar = changePstar ? view(Pstar, :, :, t) : view(Pstar, :, :)
vPstartt = changePstar ? view(Pstartt, :, :, t) : view(Pstartt, :, :)
vPstar1 = changePstar ? view(Pstar, :, :, t + 1) : view(Pstar, :, :)
vv = view(ws.v, 1:ndata)
vvH = view(vH, pattern, pattern)
......@@ -269,7 +285,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
get_cholF!(vcholH, vvH)
cholHset = true
end
ws.lik[t] += ndata*l2pi + univariate_step(Y, t, c, vZsmall, vvH, d, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, ws, pattern)
ws.lik[t] += ndata*l2pi + univariate_step(Y, t, vc, vZsmall, vvH, vd, vT, ws.QQ, va, vPinf, vPstar, diffuse_kalman_tol, kalman_tol, ws, pattern)
end
else
ws.lik[t] = ndata*l2pi + log(det_from_cholesky(vcholF))
......@@ -282,12 +298,18 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
# note that there is a typo in DK (2003) with "+ Kinf" instead of "- Kinf",
# but it is correct in their appendix
get_Kstar!(vKstar, vZsmall, vPstar, vFstar, vKinf, vcholF)
# Pstar = T*(Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar)*T'+QQ; %(5.14) DK(2012)
update_Pstar!(vPstar1, vPstar, vT, vZPinf, vZPstar, vKinf, vKstar, ws.QQ, ws.PTmp)
# Pinf = T*(Pinf-Pinf*Z'*Kinf)*T'; %(5.14) DK(2012)
update_Pinf!(vPinf1, vPinf, vT, vZPinf, vKinf, ws.PTmp)
# a = d + T*(a+Kinf*v); %(5.13) DK(2012)
update_a!(va1, va, vd, vKinf, vv, ws.a1, vT)
get_updated_a!(vatt, va, vKinf, vv)
# a = d + T*att
update_a!(va1, vd, vT, vatt)
# Pinf_tt = Pinf - Kinf'*Z*Pinf %(5.14) DK(2012)
get_updated_Ptt!(vPinftt, vPinf, vKinf, vZPinf)
# Pinf = T*Ptt*T
update_P!(vPinf1, vT, vPinftt, ws.PTmp)
# Pstartt = Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar %(5.14) DK(2012)
get_updated_Pstartt!(vPstartt, vPstar, vZPstar, vKinf, vZPinf,
vKstar, vPinftt, ws.PTmp)
# Pinf = T*Ptt*T + QQ
update_P!(vPstar1, vT, vPstartt, ws.QQ, ws.PTmp)
end
t += 1
end
......@@ -303,8 +325,11 @@ function diffuse_kalman_filter!(Y::AbstractArray{U},
R::AbstractArray{U},
Q::AbstractArray{U},
a::AbstractArray{U},
att::AbstractArray{U},
Pinf::AbstractArray{U},
Pinftt::AbstractArray{U},
Pstar::AbstractArray{U},
Pstartt::AbstractArray{U},
start::V,
last::V,
presample::V,
......@@ -317,8 +342,8 @@ function diffuse_kalman_filter!(Y::AbstractArray{U},
nobs = last - start + 1
get_QQ!(ws.QQ, R, Q, ws.RQ)
lik_cst = (nobs - presample)*ny*log(2*pi)
t = diffuse_kalman_filter_init!(Y, c, Z, H, d, T, R, Q, a, Pinf, Pstar, start, last, presample, tol, ws, data_pattern)
kalman_filter!(Y, c, Z, H, d, T, R, Q, a, Pstar, t + 1, last, presample, ws, data_pattern)
t = diffuse_kalman_filter_init!(Y, c, Z, H, d, T, R, Q, a, att, Pinf, Pinftt, Pstar, Pstartt, start, last, presample, tol, ws, data_pattern)
kalman_filter!(Y, c, Z, H, d, T, R, Q, a, att, Pstar, Pstartt, t + 1, last, presample, ws, data_pattern)
vlik = view(ws.lik, start + presample:last)
return -0.5*sum(vlik)
end
......@@ -332,8 +357,11 @@ function diffuse_kalman_filter!(Y::AbstractArray{U},
R::AbstractArray{U},
Q::AbstractArray{U},
a::AbstractArray{U},
att::AbstractArray{U},
Pinf::AbstractArray{U},
Pinftt::AbstractArray{U},
Pstar::AbstractArray{U},
Pstartt::AbstractArray{U},
start::V,
last::V,
presample::V,
......@@ -344,5 +372,6 @@ function diffuse_kalman_filter!(Y::AbstractArray{U},
m, n = size(Y)
full_data_pattern = [collect(1:m) for i = 1:n]
diffuse_kalman_filter!(Y, c, Z, H, d, T, R, Q, a, Pinf, Pstar, start, last, presample, tol, ws, full_data_pattern)
lik = diffuse_kalman_filter!(Y, c, Z, H, d, T, R, Q, a, att, Pinf, Pinftt, Pstar, Pstartt, start, last, presample, tol, ws, full_data_pattern)
return lik
end
......@@ -207,7 +207,7 @@ function kalman_filter_2!(Y::AbstractArray{U},
# P_{t|t} = P_t - K'*Z*P_t
filtered_P!(vPt, vP, vK, vZP, ws.PTmp)
# P_{t+1} = T*P_{t|t}*T'+ QQ
update_P!(vP1, vPt, vT, ws.QQ, ws.PTmp)
update_P!(vP1, vT, vPt, ws.QQ, ws.PTmp)
#=
ws.oldP .-= vP1
......
......@@ -312,10 +312,10 @@ function univariate_step(Y, t, c, Z, H, d, T, QQ, a, Pinf, Pstar, diffuse_kalman
end
llik = 0.0
l2pi = log(2*pi)
ndata = size(pattern)
ndata = length(pattern)
for i=1:ndata
Zi = view(ws.Zstar, pattern[i], :)
v = get_v(ws.ystar, c, ws.Zstar, a, pattern[i])
v = get_v!(ws.ystar, c, ws.Zstar, a, pattern[i])
Fstar = get_Fstar!(Zi, Pstar, H[i], ws.uKstar)
Finf = get_Finf!(Zi, Pstar, ws.uKstar)
# Conduct check of rank
......@@ -326,7 +326,7 @@ function univariate_step(Y, t, c, Z, H, d, T, QQ, a, Pinf, Pstar, diffuse_kalman
# Also the test for singularity is better set coarser for Finf than for Fstar for the same reason
if Finf > diffuse_kalman_tol # F_{\infty,t,i} = 0, use upper part of bracket on p. 175 DK (2012) for w_{t,i}
ws.Kinf_Finf .= ws.uKinf./Finf
a .+= prediction_error.*ws.Kinf_Finf
a .+= v .* ws.Kinf_Finf
# Pstar = Pstar + Kinf*(Kinf_Finf'*(Fstar/Finf)) - Kstar*Kinf_Finf' - Kinf_Finf*Kstar'
ger!( Fstar/Finf, ws.uKinf, ws.Kinf_Finf, Pstar)
ger!( -1.0, ws.uKstar, ws.Kinf_Finf, Pstar)
......@@ -346,8 +346,10 @@ function univariate_step(Y, t, c, Z, H, d, T, QQ, a, Pinf, Pstar, diffuse_kalman
copy!(ws.a1, d)
mul!(ws.a1, T, a, 1.0, 1.0)
a .= ws.a1
mul!(ws.PTmp, T, P)
copy!(P, RQR)
mul!(P, ws.PTmp, T', 1.0, 1.0)
mul!(ws.PTmp, T, Pinf)
mul!(Pinf, ws.PTmp, transpose(T))
mul!(ws.PTmp, T, Pstar)
copy!(Pstar, QQ)
mul!(Pstar, ws.PTmp, transpose(T), 1.0, 1.0)
return llik + 2*log(detLTcholH)
end
......@@ -63,10 +63,13 @@ H_0 = copy(H)
H = zeros(ny, ny) + I(ny)
ws1 = KalmanLikelihoodWs(ny, ns, np, nobs)
P = zeros(ns, ns, nobs+1)
s = copy(s_0)
llk_1 = kalman_filter!(y, zeros(ny), Z, H, zeros(ns), T, R, Q, s, P, 1, nobs, 0, ws1, full_data_pattern)
P = zeros(ns, ns, nobs + 1)
s = zeros(ns, nobs + 1)
s[:, 1] .= s_0
Ptt = zeros(ns, ns, nobs)
stt = zeros(ns, nobs)
llk_1 = kalman_filter!(y, zeros(ny), Z, H, zeros(ns), T, R, Q, s, stt, P, Ptt, 1, nobs, 0, ws1, full_data_pattern)
@test P[:, :, 2] R*Q*R'
P = zeros(ns, ns)
......@@ -101,9 +104,12 @@ end
ws1 = KalmanLikelihoodWs(ny, ns, np, nobs)
P = zeros(ns, ns, nobs+1)
s = copy(s_0)
s = zeros(ns, nobs + 1)
s[:, 1] .= s_0
Ptt = zeros(ns, ns, nobs)
stt = zeros(ns, nobs)
llk_1 = kalman_filter!(y, zeros(ny), Z, H, zeros(ns), T, R, Q, s, P, 1, nobs, 0, ws1, full_data_pattern)
llk_1 = kalman_filter!(y, zeros(ny), Z, H, zeros(ns), T, R, Q, s, stt, P, Ptt, 1, nobs, 0, ws1, full_data_pattern)
@test P[:, :, 2] R*Q*R'
P = zeros(ns, ns)
......@@ -197,11 +203,14 @@ end
c = zeros(ny)
d = zeros(ns)
s = copy(s_0)
stt = similar(s)
P = copy(P_0)
Ptt = similar(P)
nobs1 = 1
ws1 = KalmanLikelihoodWs{Float64, Int64}(ny, ns, np, nobs1)
kalman_filter!(y, c, Z, H, d, T, R, Q, s, P, 1, nobs1, 0, ws1, full_data_pattern)
kalman_filter!(y, c, Z, H, d, T, R, Q, s, stt, P, Ptt, 1, nobs1, 0, ws1, full_data_pattern)
cs = zeros(ny, nobs)
Zs = zeros(ny, ns, nobs)
......@@ -211,7 +220,9 @@ end
Rs = zeros(ns, np, nobs)
Qs = zeros(np, np, nobs)
ss = zeros(ns, nobs + 1)
stt = zeros(ns, nobs)
Ps = zeros(ns, ns, nobs + 1)
Ptt = zeros(ns, ns, nobs)
for i = 1:nobs
cs[:, i] = c
......@@ -226,7 +237,7 @@ end
ss[:, 1] = s_0
Ps[:, :, 1] = P_0
kalman_filter!(y, cs, Zs, Hs, ds, Ts, Rs, Qs, ss, Ps, 1, nobs1, 0, ws1, full_data_pattern)
kalman_filter!(y, cs, Zs, Hs, ds, Ts, Rs, Qs, ss, stt, Ps, Ptt, 1, nobs1, 0, ws1, full_data_pattern)
@test ss[:, nobs1+1] s
@test Ps[:, : , nobs1+1] P
......@@ -283,7 +294,11 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
Pstar = copy(Pstar_0)
copy!(ws4.QQ, R*Q*R')
t = KalmanFilterTools.diffuse_kalman_likelihood_init!(Y, Z, H, T, ws4.QQ, a, Pinf, Pstar, 1, nobs, 1e-8, ws4)
t = KalmanFilterTools.diffuse_kalman_likelihood_init!(Y, Z, H, T,
ws4.QQ, a,
Pinf, Pstar,
1, nobs,
1e-8, ws4)
llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
# Dynare returns minus log likelihood
......@@ -293,10 +308,15 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
c = zeros(ny)
d = zeros(ns)
aa = repeat(a_0, 1, nobs)
PPinf = repeat(Pinf_0, 1, 1, nobs)
PPstar = repeat(Pstar_0, 1, 1, nobs)
t1 = KalmanFilterTools.diffuse_kalman_filter_init!(Y, c, Z, H, d, T, R, Q, aa, PPinf, PPstar, 1, nobs, 0, 1e-8, ws5, full_data_pattern)
aa = repeat(a_0, 1, nobs+1)
att = similar(aa)
PPinf = repeat(Pinf_0, 1, 1, nobs+1)
Pinftt = similar(PPinf)
PPstar = repeat(Pstar_0, 1, 1, nobs+1)
Pstartt = similar(PPstar)
t1 = KalmanFilterTools.diffuse_kalman_filter_init!(Y, c, Z, H, d, T, R, Q, aa, att,
PPinf, Pinftt, PPstar, Pstartt,
1, nobs, 0, 1e-8, ws5, full_data_pattern)
@test t1 == t
@test llk_3 -0.5*sum(ws5.lik[1:t1])
@test aa[:, t1 + 1] vars["a"]
......@@ -306,7 +326,8 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
a = copy(a_0)
Pinf = copy(Pinf_0)
Pstar = copy(Pstar_0)
t = KalmanFilterTools.diffuse_kalman_likelihood_init!(Y, z, H, T, ws4.QQ, a, Pinf, Pstar, 1, nobs, 1e-8, ws4)
t = KalmanFilterTools.diffuse_kalman_likelihood_init!(Y, z, H, T, ws4.QQ, a, Pinf,
Pstar, 1, nobs, 1e-8, ws4)
llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
# Dynare returns minus log likelihood
......@@ -317,7 +338,8 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
aa = repeat(a_0, 1, nobs)
PPinf = repeat(Pinf_0, 1, 1, nobs)
PPstar = repeat(Pstar_0, 1, 1, nobs)
t1 = KalmanFilterTools.diffuse_kalman_filter_init!(Y, c, z, H, d, T, R, Q, aa, PPinf, PPstar, 1, nobs, 0, 1e-8, ws5, full_data_pattern)
t1 = KalmanFilterTools.diffuse_kalman_filter_init!(Y, c, z, H, d, T, R, Q, aa, att, PPinf, Pinftt,
PPstar, Pstartt, 1, nobs, 0, 1e-8, ws5, full_data_pattern)
@test t1 == t
@test llk_3 -0.5*sum(ws5.lik[1:t1])
@test aa[:, t1 + 1] vars["a"]
......@@ -330,9 +352,12 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
aa = repeat(a_0, 1, nobs + 1)
Pinf = copy(Pinf_0)
Pinftt = similar(Pinf)
Pstar = copy(Pstar_0)
llk_4a = diffuse_kalman_filter!(Y, c, Z, H, d, T, R, Q, aa, Pinf, Pstar, 1, nobs, 0, 1e-8, ws5)
Pstartt = similar(Pstar)
llk_4a = diffuse_kalman_filter!(Y, c, Z, H, d, T, R, Q, aa, att, Pinf, Pinftt, Pstar, Pstartt,
1, nobs, 0, 1e-8, ws5)
@test llk_4a llk_4
a = copy(a_0)
......@@ -344,7 +369,7 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
aa = repeat(a_0, 1, nobs + 1)
Pinf = copy(Pinf_0)
Pstar = copy(Pstar_0)
llk_5a = diffuse_kalman_filter!(Y, c, z, H, d, T, R, Q, aa, Pinf, Pstar, 1, nobs, 0, 1e-8, ws5)
llk_5a = diffuse_kalman_filter!(Y, c, z, H, d, T, R, Q, aa, att, Pinf, Pinftt, Pstar, Pstartt, 1, nobs, 0, 1e-8, ws5)
@test llk_5a llk_5
a = copy(a_0)
......
......@@ -153,6 +153,24 @@ ZW = rand(ny, ny)
KalmanFilterTools.get_M!(M, cholF, ZW)
@test M -inv(F)
# Ptt = P - K'*Z*P
Ptt = randn(ns, ns)
P - randn(ns, ns)
K = randn(ny, ns)
ZP = randn(ny, ns)
get_updated_Ptt!(Ptt, P, K, ZP)
@test Ptt P - transpose(K)*ZP
# Pstartt = Pstar-Pstar*Z'*Kinf-Pinf*Z'*Kstar %(5.14) DK(2012)
Pstartt = randn(ns, ns)
Pstar = randn(ns, ns)
ZPstar = randn(ny, ns)
Kinf = randn(ny, ns)
ZPinf = randn(ny, ns)
Kstar = randn(ny, ns)
KalmanFilterTools.get_updated_Pstartt!(Pstartt, Pstar, ZPstar, Kinf, ZPinf, Kstar)
@test Pstartt Pstar - transpose(ZPstar)*Kinf - transpose(ZPinf)*Kstar
# V_t = P_t - P_t*N_{t-1}*P_t
V = Matrix{Float64}(undef, ns, ns)
#P = Matrix{Float64}(undef, ns, ns)
......@@ -170,10 +188,18 @@ RQ = zeros(ns, np)
KalmanFilterTools.get_QQ!(QQ, R, Q, RQ)
@test QQ R*Q*R'
#v = Y[:,t] - Z*a
a = rand(ns)
# att = a + K'*v
att = randn(ns)
a = randn(ns)
K = randn(ny, ns)
v = rand(ny)
y = rand(ny, 1)
KalmanFilterTools.get_updated_a!(att, a, K, v)
@test att a + transpose(K)*v
#v = Y[:,t] - Z*a
a = randn(ns)
v = randn(ny)
y = randn(ny, 1)
KalmanFilterTools.get_v!(v, y, Z, a, 1, ny)
@test v y[:,1] - Z*a
......@@ -239,6 +265,14 @@ vT = view(T, :, :)
KalmanFilterTools.update_a!(va1, va, vd, vK, vv, a1, vT)
@test va1 vd + vT*(va + vK'*vv)
# a = d + T*att
a = randn(ns)
d = randn(ns)
T = randn(ns, ns)
att = randn(ns)
KalmanFilterTools.update_a!(a, d, T, att)
@test a d + T*att
# K = K + Z*W*M*W'
K = randn(ny, ns)
K_0 = copy(K)
......@@ -287,6 +321,14 @@ Ptmp = similar(P)
KalmanFilterTools.update_P!(vP1, vT, QQ, vK, vZP, PTmp)
@test vP1 vT*(vP - vK'*vZP)*vT' + QQ
# P = T*Ptt*T' + QQ
P = randn(ns, ns)
T = randn(ns ns)
Ptt = randn(ns, ns)
QQ = randn(ns, ns)
KalmanFilterTools.update_P!(P, T, Ptt, QQ)
@test P T*Ptt*transpose(T) + QQ
# r_{t-1} = Z_t'*iF_t*v_t + L_t'r_t
r = randn(ns)
r1 = similar(r)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment