adding tests

test/runtests.jl

+include("test_kalman_base.jl")
+include("test_KalmanFilterTools.jl")

test/test_KalmanFilterTools.jl

+using HDF5
+using LinearAlgebra
+using LinearAlgebra.BLAS
+using KalmanFilterTools
+using MAT
+using Test
+
+path = dirname(@__FILE__)
+
+# Using test case from NYFED package StateSpaceRoutines
+# Initialize arguments to function
+file = h5open("$path/reference/kalman_filter_args.h5", "r")
+y = read(file, "data")
+T, R, C    = read(file, "TTT"), read(file, "RRR"), read(file, "CCC")
+Q, Z, D, H = read(file, "QQ"), read(file, "ZZ"), read(file, "DD"), read(file, "EE")
+s_0, P_0   = read(file, "z0"), read(file, "P0")
+close(file)
+
+ny, ns = size(Z)
+nobs = size(y, 2)
+np = size(R, 2)
+
+# Removing measure equation constant from observations
+y .-= D
+
+# Create data_pattern for all observations are available
+
+full_data_pattern = [collect(1:ny) for o = 1:nobs]
+
+# Simple Kalman Filter
+P = copy(P_0)
+s = copy(s_0)
+@testset "Basic Kalman Filter" begin
+    ws1 = KalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+
+    copy!(P, P_0)
+    llk_1 = kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1)
+
+    h5open("$path/reference/kalman_filter_out.h5", "r") do h5
+        @test read(h5, "log_likelihood") ≈ llk_1
+    end
+    
+    copy!(s, s_0)
+    copy!(P, P_0)
+    llk_2 = kalman_likelihood_monitored(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1)
+    @test llk_2 ≈ llk_1
+
+    copy!(s, s_0)
+    copy!(P, P_0)
+    llk_2 = kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1, full_data_pattern)
+    @test llk_2 ≈ llk_1
+
+    copy!(s, s_0)
+    copy!(P, P_0)
+    llk_2 = kalman_likelihood_monitored(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1, full_data_pattern)
+    @test llk_2 ≈ llk_1
+
+end
+
+# Fast Kalman Filter
+@testset "Fast Kalman Filter" begin
+    ws1 = KalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+    ws2 = FastKalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+    P = copy(P_0)
+    s = copy(s_0)
+    
+
+    copy!(s, s_0)
+    llk_1 = kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1)
+
+    copy!(P, P_0)
+    copy!(s, s_0)
+    llk_2 = fast_kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws2)
+    @test llk_2 ≈ llk_1
+
+    copy!(P, P_0)
+    copy!(s, s_0)
+    llk_3 = fast_kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws2, full_data_pattern)
+    @test llk_3 ≈ llk_1
+
+end
+
+
+@testset "Z as selection matrix" begin
+    ws1 = KalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+    ws2 = FastKalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+
+    fill!(Z, 0.0)
+    Z[1, 4] = 1
+    Z[2, 3] = 1
+    Z[3, 2] = 1
+    z = [4, 3, 2]
+
+    s = copy(s_0)
+    P = copy(P_0)
+    llk_1 = kalman_likelihood(y, Z, H, T, R, Q, s, P, 1, nobs, 0, ws1)
+    
+    s = copy(s_0)
+    P = copy(P_0)
+    llk_2 = kalman_likelihood(y, z, H, T, R, Q, s, P, 1, nobs, 0, ws1)
+    @test llk_1 ≈ llk_2
+
+    s = copy(s_0)
+    P = copy(P_0)
+    llk_2 = kalman_likelihood(y, z, H, T, R, Q, s, P, 1, nobs, 0, ws1, full_data_pattern)
+    @test llk_1 ≈ llk_2
+
+    s = copy(s_0)
+    P = copy(P_0)
+    llk_3 = fast_kalman_likelihood(y, z, H, T, R, Q, s, P, 1, nobs, 0, ws2)
+    @test llk_1 ≈ llk_3
+
+    s = copy(s_0)
+    P = copy(P_0)
+    llk_3 = fast_kalman_likelihood(y, z, H, T, R, Q, s, P, 1, nobs, 0, ws2, full_data_pattern)
+    @test llk_1 ≈ llk_3
+
+end
+
+@testset "Kalman Filter" begin
+    c = zeros(ny)
+    d = zeros(ns)
+    s = copy(s_0)
+    P = copy(P_0)
+    nobs1 = 1
+    ws1 = KalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs1)
+
+    kalman_filter!(y, c, Z, H, d, T, R, Q, s, P, 1, nobs1, 0, ws1, full_data_pattern)
+    
+    cs = zeros(ny, nobs)
+    Zs = zeros(ny, ns, nobs)
+    Hs = zeros(ny, ny, nobs)
+    ds = zeros(ns, nobs)
+    Ts = zeros(ns, ns, nobs)
+    Rs = zeros(ns, np, nobs)
+    Qs = zeros(np, np, nobs)
+    ss = zeros(ns, nobs+1)
+    Ps = zeros(ns, ns, nobs)
+
+    for i = 1:nobs
+        cs[:, i] = c
+        Zs[:, :, i] = Z
+        Hs[:, :, i] = H
+        ds[:, i] = d
+        Ts[:, :, i] = T
+        Rs[:, :, i] = R
+        Qs[:, :, i] = Q
+    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)
+    @test ss[:, nobs1+1] ≈ s
+    @test Ps[:, : , nobs1+1] ≈ P
+
+end
+
+# Replication data computed with Dynare
+vars = matread("$path/reference/test_data.mat")
+
+Y = vars["Y"]
+Z = vars["Z"]
+H = vars["H"]
+T = vars["T"]
+R = vars["R"]
+Q = vars["Q"]
+Pinf_0 = vars["Pinf"]
+Pstar_0 = vars["Pstar"]
+
+ny, nobs = size(Y)
+ns, np = size(R)
+
+a_0 = zeros(ns)
+if H == 0
+    H = zeros(ny, ny)
+end
+
+full_data_pattern = [collect(1:ny) for o = 1:nobs]
+
+@testset "Diffuse Kalman Filter" begin
+    ws4 = DiffuseKalmanLikelihoodWs{Float64, Integer}(ny, ns, np, nobs)
+    
+    a = copy(a_0)
+    Pinf = copy(Pinf_0)
+    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)
+    llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
+
+    # Dynare returns minus log likelihood
+    @test llk_3 ≈ -vars["dLIK"]
+    @test a ≈ vars["a"]
+    @test Pstar ≈ vars["Pstar1"]
+
+    z = [4, 3]
+    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)
+    llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
+
+    # Dynare returns minus log likelihood
+    @test llk_3 ≈ -vars["dLIK"]
+    @test a ≈ vars["a"]
+    @test Pstar ≈ vars["Pstar1"]
+
+    a = copy(a_0)
+    Pinf = copy(Pinf_0)
+    Pstar = copy(Pstar_0)
+    llk_4 = diffuse_kalman_likelihood(Y, Z, H, T, R, Q, a, Pinf, Pstar, 1, nobs, 0, 1e-8, ws4)
+
+    a = copy(a_0)
+    Pinf = copy(Pinf_0)
+    Pstar = copy(Pstar_0)
+    llk_5 = diffuse_kalman_likelihood(Y, z, H, T, R, Q, a, Pinf, Pstar, 1, nobs, 0, 1e-8, ws4)
+    @test llk_5 ≈ llk_4 
+
+    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, full_data_pattern)
+    llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
+
+    # Dynare returns minus log likelihood
+    @test llk_3 ≈ -vars["dLIK"]
+    @test a ≈ vars["a"]
+    @test Pstar ≈ vars["Pstar1"]
+
+    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, full_data_pattern)
+    llk_3 = -0.5*(t*ny*log(2*pi) + sum(ws4.lik[1:t]))
+
+    # Dynare returns minus log likelihood
+    @test llk_3 ≈ -vars["dLIK"]
+    @test a ≈ vars["a"]
+    @test Pstar ≈ vars["Pstar1"]
+
+    a = copy(a_0)
+    Pinf = copy(Pinf_0)
+    Pstar = copy(Pstar_0)
+    llk_4 = diffuse_kalman_likelihood(Y, Z, H, T, R, Q, a, Pinf, Pstar, 1, nobs, 0, 1e-8, ws4, full_data_pattern)
+
+    a = copy(a_0)
+    Pinf = copy(Pinf_0)
+    Pstar = copy(Pstar_0)
+    llk_5 = diffuse_kalman_likelihood(Y, z, H, T, R, Q, a, Pinf, Pstar, 1, nobs, 0, 1e-8, ws4, full_data_pattern)
+    @test llk_5 ≈ llk_4 
+end
+
+nothing

test/test_kalman_base.jl

+using LinearAlgebra
+using LinearAlgebra.BLAS
+using KalmanFilterTools
+using Test
+
+ny = 3
+ns = 10
+np   = 2
+nobs = 50
+
+Z = randn(ny, ns)
+T = randn(ns, ns)
+K = randn(ny, ns)
+L = Matrix{Float64}(undef, ns, ns)
+L1 = similar(L)
+
+KalmanFilterTools.get_L!(L, T, K, Z, L1)
+@test L ≈ T - T*K'*Z
+
+z = [1, 3]
+K1 = K[z, :]
+K2 = zeros(ns, ns)
+K2[z,:] .= K1
+KalmanFilterTools.get_L!(L, T, K1, z, L1)
+@test L ≈ T - T*K2'
+
+# r_{t-1} = Z_t'*iF_t*v_t + L_t'r_t
+r = randn(ns)
+r1 = similar(r)
+iFv = randn(ny)
+
+KalmanFilterTools.update_r!(r1, Z, iFv, L, r)
+@test r1 ≈ Z'*iFv + L'r
+
+# r_{t-1} = Z_t'*iF_t*v_t + L_t'r_t
+iFv1 = iFv[z]
+KalmanFilterTools.update_r!(r1, z, iFv1, L, r)
+ZiF = zeros(ns)
+ZiF[z] .= iFv1
+@test r1 ≈ ZiF + L'r
+
+# alphah_t = a_t + P_t*r_{t-1}
+alphah = Vector{Float64}(undef, ns)
+a = randn(ns)
+P = randn(ns, ns)
+P = P'*P
+KalmanFilterTools.get_alphah!(alphah, a, P, r1)
+@test alphah ≈ a + P*r1
+
+# N_{t-1} = Z_t'iF_t*Z_t + L_t'N_t*L_t
+N = randn(ns, ns)
+N = N'*N
+N1 = similar(N)
+iFZ = randn(ny, ns)
+Ptmp = similar(P)
+KalmanFilterTools.update_N!(N1, Z, iFZ, L, N, Ptmp)
+@test N1 ≈ Z'*iFZ + L'*N*L
+
+iFZ1 = iFZ[z, :]
+ZiFZ = zeros(ns, ns)
+ZiFZ[z, :] .= iFZ1
+KalmanFilterTools.update_N!(N1, z, iFZ1, L, N, Ptmp)
+@test N1 ≈ ZiFZ + L'*N*L
+
+# V_t = P_t - P_t*N_{t-1}*P_t
+V = Matrix{Float64}(undef, ns, ns)
+KalmanFilterTools.get_V!(V, P, N1, Ptmp)
+@test V ≈ P - P*N1*P
+
+R = randn(ns, np)
+Q = randn(np, np)
+Q = Q'*Q
+QQ = zeros(ns, ns)
+RQ = zeros(ns, np)
+KalmanFilterTools.get_QQ!(QQ, R, Q, RQ)
+@test QQ ≈ R*Q*R'
+
+H = randn(ny, ny)
+H = H'*H
+F = zeros(ny, ny)
+ZP = zeros(ny, ns)
+KalmanFilterTools.get_F!(F, ZP, Z, P, H)
+@test ZP == Z*P
+@test F ≈ Z*P*Z' + H
+
+cholF = zeros(ny, ny)
+KalmanFilterTools.get_cholF!(cholF, F)
+CF = cholesky(0.5*(F + F'))
+@test triu(cholF) ≈ CF.U
+
+K = zeros(ny, ns)
+KalmanFilterTools.get_K!(K, ZP, cholF)
+@test K ≈ inv(F)*Z*P
+
+P_0 = similar(P)
+PTmp = similar(P)
+copy!(P_0, P)
+KalmanFilterTools.update_P!(P, T, QQ, K, ZP, PTmp)
+@test P ≈ T*(P_0 - K'*ZP)*T' + QQ
+
+W = rand(ns, ny)
+K = Z*P
+mul!(W, T, transpose(K))
+@test W ≈ T*P*Z'
+
+M = rand(ny, ny)
+ZW = rand(ny, ny)
+KalmanFilterTools.get_M!(M, cholF, ZW) 
+@test M ≈ -inv(F)
+
+#v  = Y[:,t] - Z*a
+a = rand(ns)
+v = rand(ny)
+y = rand(ny, 1)
+KalmanFilterTools.get_v!(v, y, Z, a, 1, ny)
+@test v ≈ y[:,1] - Z*a
+
+# iFv = inv(F)*v
+iFv = similar(v)
+KalmanFilterTools.get_iFv!(iFv, cholF, v)
+@test iFv ≈ F\v
+
+# a = T(a + K'*iFv)
+a_0 = copy(a)
+a1 = similar(a)
+KalmanFilterTools.update_a!(a, K, iFv, a1, T)
+@test a ≈ T*(a_0 + K'*iFv)
+
+# M = M + M*W'*Z'iF*Z*W*M
+M_0 = copy(M)
+ZWM = similar(M)
+iFZWM = similar(M)
+KalmanFilterTools.update_M!(M, Z, W, cholF, ZW, ZWM, iFZWM)
+@test M ≈ M_0 + M_0*W'*Z'*inv(F)*Z*W*M_0
+
+# F =  F + Z*W*M*W'Z'
+F_0 = copy(F)
+gemm!('N', 'T', 1.0, ZWM, ZW, 1.0,F)
+@test F ≈ F_0 + Z*W*M_0*W'*Z'
+
+# K = K + Z*W*M*W'
+K_0 = copy(K)
+KalmanFilterTools.update_K!(K, ZWM, W)
+@test K ≈ K_0 + Z*W*M_0*W'
+
+# W = T(W - K'*iF*Z*W)
+#K2 = K1  + T*W*M_0*W'*Z'
+W_0 = copy(W)
+ZW = Z*W
+F = randn(ny, ny)
+F = F'*F
+cholF = copy(F)
+LAPACK.potrf!('U', cholF)
+iFZW = rand(ny, ny)
+copy!(iFZW, ZW)
+LAPACK.potrs!('U', cholF, iFZW)
+@test iFZW ≈ inv(F)*ZW
+
+KtiFZW = rand(ns, ny)
+KalmanFilterTools.update_W!(W, ZW, cholF, T, K, iFZW, KtiFZW)
+@test iFZW ≈ inv(F)*ZW
+@test KtiFZW ≈ W_0 - K'*inv(F)*Z*W_0
+@test W ≈ T*(W_0 - K'*inv(F)*Z*W_0)
+#@test ws1.W ≈ T*Wold - K2*inv(ws1.F)*Z*Wold
+
+# Z as selection matrix
+fill!(Z, 0.0)
+Z[1, 4] = 1
+Z[2, 3] = 1
+Z[3, 2] = 1
+z = [4, 3, 2]
+
+P = copy(P_0)
+KalmanFilterTools.get_F!(F, ZP, z, P, H)
+@test F ≈ Z*P*Z' + H
+
+# v  = Y[:,t] - Z*a
+a = copy(a_0)
+KalmanFilterTools.get_v!(v, y, z, a, 1, ny)
+@test v ≈ y[:,1] - Z*a_0
+
+# missing observations
+vv = view(v, 1:ny)
+c = randn(ny)
+vc = view(c, 1:ny)
+vZ = view(Z, 1:ny, :)
+va = view(a_0, :)
+full_data_pattern .= [collect(1:ny)]
+pattern = full_data_pattern[1]
+
+KalmanFilterTools.get_v!(vv, y, vc, vZ, va, 1, pattern)
+@test vv ≈ y[:, 1] - c - Z*a
+
+va1 = zeros(ns)
+va = randn(ns)
+vd = zeros(ns)
+vK = randn(ny,ns)
+vv = randn(ny)
+vT = view(T, :, :)
+KalmanFilterTools.update_a!(va1, va, vd, vK, vv, a1, vT)
+@test va1  ≈ vd + vT*(va + vK'*vv)
+
+# P = T*(P - K'*Z*P)*T'+ QQ
+vP = copy(P_0)
+vZP = randn(ny, ns)
+vP1 = copy(vP)
+Ptmp = similar(P)
+KalmanFilterTools.update_P!(vP1, vT, QQ, vK, vZP, PTmp)
+@test vP1 ≈ vT*(vP - vK'*vZP)*vT' + QQ
+
+AA = randn(ns, ns)
+Pinf = AA*diagm(rand(ns))*AA'
+BB = randn(ns, ns)
+Pstar = BB*diagm(rand(ns))*BB'
+
+F = zeros(ny, ny)
+ZP = zeros(ny, ns)
+KalmanFilterTools.get_F!(F, ZP, Z, Pinf)
+@test F ≈ Z*Pinf*Z'
+@test ZP ≈ Z*Pinf
+
+cholF = 0.5*(F + F')
+info = LAPACK.potrf!('U', cholF)
+@test info[2] == 0
+
+K = F\Z*Pinf
+Kstar = similar(K)
+Fstar = zeros(ny, ny)
+ZPstar = zeros(ny, ns)
+KalmanFilterTools.get_F!(Fstar, ZPstar, Z, Pstar, H)
+@test Fstar ≈ Z*Pstar*Z' + H
+
+KalmanFilterTools.get_Kstar!(Kstar, Z, Pstar, Fstar, K, cholF)
+@test Kstar ≈ F\(Z*Pstar - Fstar*K)
+
+z = [4, 3, 2]
+F = 0.5*(F + F')
+cholF = copy(F)
+LAPACK.potrf!('U', cholF)
+KalmanFilterTools.get_Kstar!(Kstar, z, Pstar, Fstar, K, cholF)
+@test Kstar ≈ F\(Z*Pstar - Fstar*K)
+
+QQ = rand(ns, ns)
+RQ = rand(ns, np)
+KalmanFilterTools.get_QQ!(QQ, R, Q, RQ)
+@test QQ ≈ R*Q*R'
+