P[t+1] = T*Ptt[t]*T’ + RQR' not true in KFAS

[jye1999]

I encountered a problem while using the "alcohol" example to verify (4.24) on page 85 of the book (Koopman and Durbin, 2012).

To be specific, I wanted to get P[t+1] = TPtt[t]T’ + RQR'. But it failed.

My R code is attached.

` library("KFAS") options(digits=5)

KFAS, P-5 ----

data("alcohol") deaths <- window(alcohol[, 2], end = 2007) population <- window(alcohol[, 6], end = 2007) Zt <- matrix(c(1, 0), 1, 2) Ht <- matrix(NA) Tt <- matrix(c(1, 0, 1, 1), 2, 2) # a = a + slope Rt <- matrix(c(1, 0), 2, 1) # selection Qt <- matrix(NA) a1 <- matrix(c(0, 0), 2, 1) P1 <- matrix(0, 2, 2) P1inf <- diag(2)

model_gaussian <- SSModel(deaths / population ~ -1 +SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1, P1inf = P1inf), H = Ht) fit_g <- fitSSM(model_gaussian, inits = c(0, 0), method = "BFGS") # for Ht, Qt logLik(fit_g$model) #[1] -108.9734

out_g <- KFS(fit_g$model)

attach(out_g) # database is attached to the R search path str(a) # Time-Series [1:40, 1:2] from 1969 to 2008: 0 23.7 20.6 27 26.3 ... str(att) # Time-Series [1:39, 1:2] from 1969 to 2007: 23.7 22.2 25.6 25.5 20.2 ... str(P) # num [1:2, 1:2, 1:40] 0 0 0 0 13.7 ... str(Ptt) # num [1:2, 1:2, 1:39] 9.49 0 0 0 9.49 ...

n=dim(att)[1]

H= c(model$H); H # [1] 9.4884 Q= c(model$Q); Q #[1] 4.257 Q=diag(c(Q,0)) # already RQR'

t=n-2;

check Koopman P-85, (4.24)

Ft = Z Pt Z' + H

F= Zt %% P[,,t] %% t(Zt) +H ; F # 18.877

att = at + PZ' F-¹ * vt,

att[t,] - a[t,] # -0.919809 -0.026336 da = P[,,t] %% t(Zt)/c(F) v[t]; c(da) # -0.919809 -0.026336

P[t+1] = TP[tt]T’ + Q

P[,,t+1] - Tt%% Ptt[,,t]%% t(Tt) - Q # already RQR'

[,1] [,2]

[1,] -0.133697 -0.066848 -- ?

[2,] -0.066848 0.000000

detach(out_g)` kp_alco_Ptt.r.txt

[helske]

Thanks, there seems to be a bug in recent contemporaneous filtering algorithm, I'll look into this next week. If you are in a hurry, you can use bssm package with as_gssm function which can convert KFAS model to bssm model and then you can use kfilter and smoother from bssm.

[jye1999]

Thank you, I will try.

[helske]

The issue was in copying Pt to Ptt where I had forgotten that only upper triangular part of Pt was updated at that point. So the strictly lower triangular part of Ptt was wrong. Fixed now in the Github version, should be in CRAN as well soon.