Skip to content
GitLab
Commit 5d975ceb authored by MichelJuillard's avatar MichelJuillard
Browse files

fixed diffuse filter and smoother

parent 3e0e216c
Hide whitespace changes
Inline Side-by-side
src/kalman_base.jl
View file @ 5d975ceb
......@@ -206,7 +206,7 @@ function get_L!(L::AbstractArray{U}, T::AbstractArray{U}, K::AbstractArray{U}, z
for j = 1:length(z)
zj = z[j]
for i = 1:size(L, 1)
L[i, zj] += K[i, j]
L[i, zj] -= K[i, j]
end
end
end
......@@ -238,8 +238,12 @@ function get_L!(L::AbstractArray{U}, T::AbstractArray{U}, K::AbstractArray{U}, z
end
# L1_t = - KDK*Z (DK 5.12)
function get_L1!(L1::AbstractMatrix{T}, KDK::AbstractMatrix{T}, Zsmall::AbstractVector{U}) where {T <: AbstractFloat, U <: Integer}
vL1 = view(L1, :, Zsmall)
function get_L1!(L1::AbstractMatrix{T}, KDK::AbstractMatrix{T}, Z::AbstractMatrix{T}) where T <: AbstractFloat
mul!(L1, KDK, Z, -1.0, 0.0)
end
function get_L1!(L1::AbstractMatrix{T}, KDK::AbstractMatrix{T}, z::AbstractVector{U}) where {T <: AbstractFloat, U <: Integer}
vL1 = view(L1, :, z)
vL1 .= -KDK
end
......@@ -625,25 +629,25 @@ function update_r!(r1::AbstractVector{T}, z::AbstractVector{U}, iFv::AbstractVec
gemm!('T', 'N', 1.0, L, r, 1.0, r1)
end
# rstar_{t-1} = Z_t'*iFinf_t*v_t + Linf_t'rstar_t + Lstar_t'*rinf_t (DK 5.21)
function update_r!(rstar::AbstractVector{T}, Z::AbstractMatrix{T},
iFv::AbstractVector{T}, Linf::AbstractMatrix{T},
rstar1::AbstractVector{T}, Lstar::AbstractMatrix{T},
rinf1::AbstractVector{T}) where T <: AbstractFloat
mul!(rstar, transpose(Z), iFv)
mul!(rstar, transpose(Linf), rstar1, 1.0, 1.0)
mul!(rstar, transpose(Lstar), rinf1, 1.0, 1.0)
end
# rstar_{t-1} = z_t'*iFinf_t*v_t + Linf_t'rstar_t + Lstar_t'*rinf_t (DK 5.21)
function update_r!(rstar::AbstractVector{T}, z::AbstractVector{U},
iFv::AbstractVector{T}, Linf::AbstractMatrix{T},
rstar1::AbstractVector{T}, Lstar::AbstractMatrix{T},
rinf1::AbstractVector{T}) where {T <: AbstractFloat, U <: Integer}
vrstar = view(rstar, z)
vrstar .= iFv
mul!(rstar, transpose(Linf), rstar1, 1.0, 1.0)
mul!(rstar, transpose(Lstar), rinf1, 1.0, 1.0)
# r1_{t-1} = Z_t'*iFinf_t*v_t + L0_t'r1_t + L1_t'*0_t (DK 5.21)
function update_r!(r1::AbstractVector{T}, Z::AbstractMatrix{T},
iFv::AbstractVector{T}, L0::AbstractMatrix{T},
r1_1::AbstractVector{T}, L1::AbstractMatrix{T},
r0_1::AbstractVector{T}) where T <: AbstractFloat
mul!(r1, transpose(Z), iFv)
mul!(r1, transpose(L0), r1_1, 1.0, 1.0)
mul!(r1, transpose(L1), r0_1, 1.0, 1.0)
end
# r1_{t-1} = z_t'*iFinf_t*v_t + L0_t'r1_t + L1_t'*r0_t (DK 5.21)
function update_r!(r1::AbstractVector{T}, z::AbstractVector{U},
iFv::AbstractVector{T}, L0::AbstractMatrix{T},
r1_1::AbstractVector{T}, L1::AbstractMatrix{T},
r0_1::AbstractVector{T}) where {T <: AbstractFloat, U <: Integer}
vr1 = view(r1, z)
vr1 .= iFv
mul!(r1, transpose(L0), r1_1, 1.0, 1.0)
mul!(r1, transpose(L1), r0_1, 1.0, 1.0)
end
# W = T(W - K'*iF*Z*W)
......
src/kalman_filter.jl
View file @ 5d975ceb
......@@ -8,7 +8,7 @@ struct KalmanFilterWs{T, U} <: KalmanWs{T, U}
QQ::Matrix{T}
v::Matrix{T}
F::Matrix{T}
cholF::Matrix{T}
cholF::Array{T}
cholH::Matrix{T}
iF::Array{T}
iFv::Array{T}
......@@ -49,7 +49,7 @@ struct KalmanFilterWs{T, U} <: KalmanWs{T, U}
RQ = Matrix{T}(undef, ns, np)
QQ = Matrix{T}(undef, ns, ns)
F = Matrix{T}(undef, ny, ny)
cholF = Matrix{T}(undef, ny, ny)
cholF = Array{T}(undef, ny, ny, nobs)
cholH = Matrix{T}(undef, ny, ny)
iF = Array{T}(undef, ny, ny, nobs)
a1 = Vector{T}(undef, ns)
......@@ -159,11 +159,11 @@ function kalman_filter!(Y::AbstractArray{U},
end
viFv = changeiFv ? view(ws.iFv, 1:ndata, t) : view(ws.iFv, 1:ndata)
vv = view(ws.v, 1:ndata)
vv = view(ws.v, 1:ndata, t)
vF = view(ws.F, 1:ndata, 1:ndata)
vvH = view(vH, pattern, pattern)
vZP = view(ws.ZP, 1:ndata, :)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, 1)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
# v = Y[:,t] - c - Z*a
get_v!(vv, Y, vc, vZsmall, va, t, pattern)
......@@ -225,7 +225,7 @@ struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U}
QQ::Matrix{T}
RQ::Matrix{T}
c::Vector{T}
v::Vector{T}
v::Matrix{T}
F::Matrix{T}
iF::Matrix{T}
iFv::Matrix{T}
......@@ -257,7 +257,7 @@ struct DiffuseKalmanFilterWs{T, U} <: KalmanWs{T, U}
QQ = Matrix{T}(undef, ns, ns)
RQ = Matrix{T}(undef, ns, np)
c = Vector{T}(undef, ny)
v = Vector{T}(undef, ny)
v = Matrix{T}(undef, ny, nobs)
F = Matrix{T}(undef, ny, ny)
iF = Matrix{T}(undef, ny,ny )
iFv = Matrix{T}(undef, ny, nobs)
......@@ -354,17 +354,17 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
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)
vv = view(ws.v, 1:ndata, t)
vvH = view(vH, pattern, pattern)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
viFv = view(ws.iFv, 1:ndata)
viFv = view(ws.iFv, 1:ndata, t)
vFinf = view(ws.F, 1:ndata, 1:ndata)
vFstar = view(ws.Fstar, 1:ndata, 1:ndata)
vZPinf = view(ws.ZP, 1:ndata, :)
vZPstar = view(ws.ZPstar, 1:ndata, :)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, 1)
vcholH = view(ws.cholH, 1:ndata, 1:ndata, 1)
viFv = view(ws.iFv, 1:ndata)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
vcholH = changeH ? view(ws.cholH, 1:ndata, 1:ndata, t) : view(ws.cholH, 1:ndata, 1:ndata)
viFv = view(ws.iFv, 1:ndata, t)
vKinf = view(ws.Kinf, 1:ndata, :, 1)
vKstar = view(ws.K, 1:ndata, :, 1)
......@@ -394,8 +394,9 @@ 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)
# att = a + Kinf*v (5.13) DK(2012)
get_updated_a!(vatt, va, vKinf, vv)
# a = d + T*att
# a1 = d + T*att (5.13) DK(2012)
update_a!(va1, vd, vT, vatt)
# Pinf_tt = Pinf - Kinf'*Z*Pinf %(5.14) DK(2012)
get_updated_Ptt!(vPinftt, vPinf, vKinf, vZPinf)
......@@ -409,7 +410,7 @@ function diffuse_kalman_filter_init!(Y::AbstractArray{U},
end
t += 1
end
t
return t
end
function diffuse_kalman_filter!(Y::AbstractArray{U},
......
src/kalman_likelihood.jl
View file @ 5d975ceb
......@@ -285,7 +285,7 @@ function kalman_likelihood_monitored(Y::AbstractArray{U},
if !steady
# F = Z*P*Z' + H
get_F!(vF, vZP, Z, P, vH)
info = get_cholF!(ws.cholF, ws.F)
info = get_cholF!(vcholF, ws.F)
if info != 0
# F is near singular
if !cholHset
......@@ -656,6 +656,7 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U},
t = start
iy = (start - 1)*ny + 1
diffuse_kalman_tol = 1e-8
l2pi = log(2*pi)
kalman_tol = 1e-8
while t <= last
pattern = data_pattern[t]
......@@ -685,7 +686,7 @@ function diffuse_kalman_likelihood_init!(Y::Matrix{U},
ws.lik[t] += ndata*l2pi + univariate_step(Y, t, vZsmall, H, T, QQ, a, Pinf, Pstar, diffuse_kalman_tol, kalman_tol, pattern, ws, pattern)
end
else
ws.lik[t] = log(det_from_cholesky(ws.cholF))
ws.lik[t] = ndata*l2pi + log(det_from_cholesky(vcholF))
# Kinf = iFinf*Z*Pinf %define Kinf'=T^{-1}*K_0 with M_{\infty}=Pinf*Z'
copy!(vK, vZP)
LAPACK.potrs!('U', vcholF, vK)
......
src/kalman_smoother.jl
View file @ 5d975ceb
......@@ -8,7 +8,7 @@ struct KalmanSmootherWs{T, U} <: KalmanWs{T, U}
QQ::Matrix{T}
v::Matrix{T}
F::Matrix{T}
cholF::Matrix{T}
cholF::Array{T}
cholH::Matrix{T}
iF::Array{T}
iFv::Array{T}
......@@ -49,7 +49,7 @@ struct KalmanSmootherWs{T, U} <: KalmanWs{T, U}
RQ = Matrix{T}(undef, ns, np)
QQ = Matrix{T}(undef, ns, ns)
F = Matrix{T}(undef, ny, ny)
cholF = Matrix{T}(undef, ny, ny)
cholF = Array{T}(undef, ny, ny, nobs)
cholH = Matrix{T}(undef, ny, ny)
iF = Array{T}(undef, ny, ny, nobs)
a1 = Vector{T}(undef, ns)
......@@ -137,7 +137,7 @@ function kalman_smoother!(Y::AbstractArray{U},
fill!(ws.N_1,0.0)
ny = size(Y, 1)
for t = last: -1: 1
for t = last: -1: start
#inputs
pattern = data_pattern[t]
ndata = length(pattern)
......@@ -154,8 +154,10 @@ function kalman_smoother!(Y::AbstractArray{U},
vR = changeR ? view(R, :, :, t) : view(R, :, :)
vK = view(ws.K, 1:ndata, :, t)
vKDK = view(ws.KDK, :, 1:ndata, 1) # amounts to K_t (4.22): here KDK = T*K'
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
viF = view(ws.iF, 1:ndata, 1:ndata, t)
viFv = view(ws.iFv, 1:ndata, t)
viFZ = view(ws.iFZ, 1:ndata, :)
mul!(vKDK, vT, transpose(vK))
......@@ -167,8 +169,8 @@ function kalman_smoother!(Y::AbstractArray{U},
length(epsilonh) > 0 ||
length(etah) > 0)
# N_{t-1} = Z_t'iF_t*Z_t + L_t'N_t*L_t (DK 4.44)
get_iFZ!(ws.iFZ, viF, vZsmall)
update_N!(ws.N, vZsmall, ws.iFZ, ws.L, ws.N_1, ws.PTmp)
get_iFZ!(viFZ, vcholF, vZsmall)
update_N!(ws.N, vZsmall, viFZ, ws.L, ws.N_1, ws.PTmp)
end
if length(epsilonh) > 0
vepsilonh = view(epsilonh, :, t)
......@@ -176,6 +178,7 @@ function kalman_smoother!(Y::AbstractArray{U},
get_epsilonh!(vepsilonh, vH, viFv, vKDK, ws.r_1, ws.tmp_ny, ws.tmp_ns)
if length(Vepsilon) > 0
vVepsilon = view(Vepsilon,:,:,t)
get_iF!(viF, vcholF)
# D_t = inv(F_t) + KDK_t'*N_t*KDK_t (DK 4.69)
get_D!(ws.D, viF, vKDK, ws.N_1, ws.tmp_ny_ns)
# Vepsilon_t = H - H*D_t*H (DK 4.69)
......@@ -389,7 +392,6 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
N2 = ws.N2
N2_1 = ws.N2_1
fill!(r0_1, 0.0)
fill!(r1_1, 0.0)
ny = size(Y, 1)
......@@ -418,7 +420,7 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
vFstar =view(ws.Fstar, 1:ndata, 1:ndata, t)
viF = view(ws.iF, 1:ndata, 1:ndata, t)
vcholF = view(ws.cholF, 1:ndata, 1:ndata, t)
viFv = view(ws.iFv, 1:ndata)
viFv = view(ws.iFv, 1:ndata, t)
mul!(vKDKinf, vT, transpose(vKinf))
mul!(vKDK, vT, transpose(vK))
......@@ -429,9 +431,9 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
get_L!(L0, vT, vKDKinf, vZsmall)
# L1_t = - KDK*Z (DK 5.12)
get_L1!(L1, vKDK, vZsmall)
# r_{t-1} = Z_t'*iFinf_t*v_t + L0_t'r_t + L1_t'*r0_t (DK 5.21)
update_r!(r1, vZsmall, viFv, L1, r1_1, L0, r0_1)
# rinf_{t-1} = L0_t'rinf_t (DK 5.21)
# r1_{t-1} = Z_t'*iFinf_t*v_t + L0_t'r1_t + L1_t'*r0_t (DK 5.21)
update_r!(r1, vZsmall, viFv, L0, r1_1, L1, r0_1)
# r0_{t-1} = L0_t'r0_t (DK 5.21)
mul!(r0, L0, r0_1)
if (length(alphah) > 0 ||
length(epsilonh) > 0 ||
......@@ -478,7 +480,7 @@ function diffuse_kalman_smoother_coda!(Y::AbstractArray{U},
if length(alphah) > 0
valphah = view(alphah, :, t)
# alphah_t = a_t + Pstar_t*r0_{t-1} + Pinf_t*r1_{t-1} (DK 5.24)
get_alphah!(valphah, va, vPstar, vPinf, r0, r1,)
get_alphah!(valphah, va, vPstar, vPinf, r0, r1)
end
if length(Valpha) > 0
......
test/test_KalmanFilterTools.jl
View file @ 5d975ceb
......@@ -333,7 +333,7 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
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]))
@show llk_3
# Dynare returns minus log likelihood
@test llk_3 ≈ -vars["dLIK"]
@test a ≈ vars["a"]
......@@ -381,7 +381,7 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
Pstar = copy(Pstar_0)
t = KalmanFilterTools.diffuse_kalman_likelihood_init!(Y, Z, H, T, ws4.QQ, a, Pinf, Pstar, 1, nobs, 1e-8, ws4, full_data_pattern)
llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
llk_3 = -0.5*sum(ws4.lik[1:t])
# Dynare returns minus log likelihood
@test llk_3 ≈ -vars["dLIK"]
......@@ -392,7 +392,7 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
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, full_data_pattern)
llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
llk_3 = -0.5*sum(ws4.lik[1:t])
# Dynare returns minus log likelihood
@test llk_3 ≈ -vars["dLIK"]
......@@ -412,8 +412,9 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
ws6 = DiffuseKalmanSmootherWs{Float64, Int64}(ny, ns, np, nobs)
a = copy(a_0)
Pinf = zeros(ns, ns, nobs + 1)
aa = zeros(ns, nobs + 1)
aa[:, 1] .= a_0
Pinf = zeros(ns, ns, nobs + 1)
Pinftt = zeros(ns, ns, nobs + 1)
Pstar = zeros(ns, ns, nobs + 1)
Pstartt = zeros(ns, ns, nobs + 1)
......@@ -427,13 +428,37 @@ full_data_pattern = [collect(1:ny) for o = 1:nobs]
Valphah = zeros(ns, ns, nobs)
Vepsilonh = zeros(ny, ny, nobs)
Vetah = zeros(np, np, nobs)
llk_6 = diffuse_kalman_smoother!(Y, c, z, H, d, T, R, Q, aa, att,
Pinf, Pinftt, Pstar, Pstartt,
alphah, epsilonh, etah, Valphah,
Vepsilonh, Vetah, 1, nobs, 0,
1e-8, ws6)
@test llk_6 ≈ llk_4
llk_6a = diffuse_kalman_smoother!(Y, c, Z, H, d, T, R, Q, aa, att,
Pinf, Pinftt, Pstar, Pstartt,
alphah, epsilonh, etah, Valphah,
Vepsilonh, Vetah, 1, nobs, 0,
1e-8, ws6)
@test llk_6a ≈ llk_4
@test Y ≈ alphah[z, :]
aa = zeros(ns, nobs + 1)
aa[:, 1] .= a_0
Pinf = zeros(ns, ns, nobs + 1)
Pinftt = zeros(ns, ns, nobs + 1)
Pstar = zeros(ns, ns, nobs + 1)
Pstartt = zeros(ns, ns, nobs + 1)
Pinf[:, :, 1] = Pinf_0
Pinftt[:, :, 1] = Pinf_0
Pstar[:, :, 1] = Pstar_0
Pstartt[:, :, 1] = Pstar_0
alphah = zeros(ns, nobs)
epsilonh = zeros(ny, nobs)
etah = zeros(np, nobs)
Valphah = zeros(ns, ns, nobs)
Vepsilonh = zeros(ny, ny, nobs)
Vetah = zeros(np, np, nobs)
llk_6b = diffuse_kalman_smoother!(Y, c, z, H, d, T, R, Q, aa, att,
Pinf, Pinftt, Pstar, Pstartt,
alphah, epsilonh, etah, Valphah,
Vepsilonh, Vetah, 1, nobs, 0,
1e-8, ws6)
@test llk_6b ≈ llk_4
@test Y ≈ alphah[z, :]
end
@testset "start and last" begin
......@@ -441,7 +466,7 @@ end
ws1 = KalmanLikelihoodWs{Float64, Int64}(ny, ns, np, nobs)
P_0 = randn(ns, ns)
P_0 = P_0'P_0
P_0 = P_0'*P_0
a = copy(a_0)
P = copy(P_0)
......@@ -494,5 +519,6 @@ end
@testset "smoother" begin
end
nothing
test/test_kalman_base.jl
View file @ 5d975ceb
......@@ -182,6 +182,20 @@ LAPACK.potrf!('U', cholF)
KalmanFilterTools.get_Kstar!(Kstar, z, Pstar, Fstar, K, cholF)
@test Kstar ≈ F\(Pstar[z,:] - Fstar*K)
# L_t = T - K(DK)_t*Z (DK 4.29)
z = [2, 3, 1]
Z = zeros(ny, ns)
for i=1:length(z)
Z[i, z[i]] = 1.0
end
T = randn(ns, ns)
KDK = randn(ns, ny)
L = Matrix{Float64}(undef, ns, ns)
KalmanFilterTools.get_L!(L, T, KDK, Z)
@test L ≈ T - KDK*Z
KalmanFilterTools.get_L!(L, T, KDK, z)
@test L ≈ T - KDK*Z
Z = randn(ny, ns)
T = randn(ns, ns)
K = randn(ny, ns)
......@@ -197,6 +211,7 @@ K2[z,:] .= K1
KalmanFilterTools.get_L!(L, T, K1, z, L1)
@test L ≈ T - T*K2'
M = rand(ny, ny)
ZW = rand(ny, ny)
KalmanFilterTools.get_M!(M, cholF, ZW)
......
Supports Markdown
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