diff --git a/src/kalman_base.jl b/src/kalman_base.jl
index 2c2a2a20957adf0413d6d28bce5a936b3d67b414..e31013f97b809c1577f40fb629ecfd7e13c2301f 100644
--- a/src/kalman_base.jl
+++ b/src/kalman_base.jl
@@ -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)
diff --git a/src/kalman_filter.jl b/src/kalman_filter.jl
index 0c30df68ba179793adcb4b473cdd46e81d1b3007..470cdb45fe15f773a0198a0d6c2b05fd4b23f274 100644
--- a/src/kalman_filter.jl
+++ b/src/kalman_filter.jl
@@ -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},
diff --git a/src/kalman_likelihood.jl b/src/kalman_likelihood.jl
index fd394f89aa3b6cc3294537c2cb3f404d9b8a62c3..ed34d7cdf74fcdaf1e80fff05288dd3b1b7c32d8 100644
--- a/src/kalman_likelihood.jl
+++ b/src/kalman_likelihood.jl
@@ -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)
diff --git a/src/kalman_smoother.jl b/src/kalman_smoother.jl
index de909881273df3b1dffb6f13c00b65185733c42c..8ff620730f7717d9bb770a5c5fff3ec7f757ef8a 100644
--- a/src/kalman_smoother.jl
+++ b/src/kalman_smoother.jl
@@ -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
diff --git a/test/test_KalmanFilterTools.jl b/test/test_KalmanFilterTools.jl
index 0388f7ad8c92b9e332ee6407c821055dae3d371b..9fdd551010548c80a1c00588420424d6075d3d64 100644
--- a/test/test_KalmanFilterTools.jl
+++ b/test/test_KalmanFilterTools.jl
@@ -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
diff --git a/test/test_kalman_base.jl b/test/test_kalman_base.jl
index 37d7c2365de2b65e1a998e8667d29babc804e953..a61d9040e67791d2fce9ca0d7d577acc9c1acea7 100644
--- a/test/test_kalman_base.jl
+++ b/test/test_kalman_base.jl
@@ -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)