KFS smoothed variance bug (?)

[gragusa]

The following code suggests that the V_mu diagonal term returned by KFS is off by the variance of the observation equation. (The example is univariate, but for the multivariate case this is still the case, that is, the returned diagonal entry of the smoothed variance of the observation are off by the diagonal term of H).

## Test KFAS vs DLM

library(KFAS)
library(dlm)

## Generate fake data

set.seed(1)
theta <- array(0, 101)
y <- array(0, 101)
for (j in 1:100) {
  theta[j + 1] = 0.9 * theta[j] + rnorm(1)
  y[j + 1] <- theta[j + 1] + rnorm(1)
}
Y <- y[-1]

Z <- matrix(1)
T <- matrix(0.9)
R <- matrix(1)
Q <- matrix(1)
H <- matrix(1)
a1 <- matrix(0)
P1 <- matrix(1)

Y <- matrix(y[-1])

model_kfas = SSModel(Y ~ -1 + SSMcustom(Z  = Z,
                                    T  = T,
                                    R  = R,
                                    Q  = Q,
                                    a1 = a1,
                                    P1 = P1), 
                H  = H)


m <- attr(model_kfas, "m")
p <- attr(model_kfas, "p")
k <- attr(model_kfas, "k")

n.ahead <- 12
states <- as.integer(1:attr(model_kfas, "m"))
timespan <- 1:attr(model_kfas, "n")
endtime <- end(model_kfas$y)

timespan <- attr(model_kfas, "n") + 1:n.ahead
n <- attr(model_kfas, "n") <- attr(model_kfas, "n") + as.integer(n.ahead)
endtime <- end(model_kfas$y) + c(0, n.ahead)
model_kfas$y <- window(model_kfas$y, end = endtime, extend = TRUE)

PRED_KFAS <- KFS(model = model_kfas, smoothing = c("mean", "state"))

model_dlm <- dlm(m0 = a1, C0 = P1, FF = Z, V = H, GG = T, W = Q)
PRED_DLM <- dlmForecast(dlmFilter(Y, model_dlm), nAhead = 12)


all.equal(c(PRED_DLM$a), PRED_KFAS$muhat[101:(100+n.ahead)])

## This fail
all.equal(unlist(PRED_DLM$Q), PRED_KFAS$V_mu[,,101:(100+n.ahead)])

## This is correct
all.equal(unlist(PRED_DLM$Q), PRED_KFAS$V_mu[,,101:(100+n.ahead)] + H)

[helske]

The variable V_mu is not the variance of observation y_t but the variance of mean i.e. in Gaussian case Var(Z_t alpha_t | Y) so there is no bug. But to be honest the documentation is bit poor in KFS regarding that, I have to rewrite it better.

Instead of manipulating the model after running SSModel, you can just add NAs to observations when constructing it via SSModel:

model_kfas_pred <- SSModel(c(Y, rep(NA, 12)) ~ -1 + SSMcustom(Z  = Z,
  T  = T,  R  = R,  Q  = Q,  a1 = a1,  P1 = P1),   H  = H)

out <- KFS(model_kfas_pred)
c(out$V_mu) + H 

But based on your code (dlmForecast) it looks like you are trying to make make predictions, so you don't actually need smoothed variances. You can get the prediction error variances from F component:

out$F[101:112]

And you can also use predict method for SSModel object (I am going to add similar method to KFS object soon too), for example:

predict(model_kfas, n.ahead = 12, se.fit = TRUE, interval = "prediction")
Time Series:
Start = 101 
End = 112 
Frequency = 1 
          fit       lwr      upr   se.fit
101 0.3848579 -2.704122 3.473838 1.218154
102 0.3463721 -3.160791 3.853535 1.483900
103 0.3117349 -3.500678 4.124148 1.668408
104 0.2805614 -3.762244 4.323367 1.804080
105 0.2525053 -3.967708 4.472718 1.906911
106 0.2272547 -4.131369 4.585878 1.986306
107 0.2045293 -4.263064 4.672123 2.048361
108 0.1840763 -4.369872 4.738024 2.097280
109 0.1656687 -4.457044 4.788381 2.136084
110 0.1491018 -4.528570 4.826773 2.167005
111 0.1341916 -4.587527 4.855911 2.191732
112 0.1207725 -4.636326 4.877871 2.211558

[gragusa]

The documentation says the variance of the link function so I thought it
was the prediction variance.

I am implementing univariate filtering in Julia and I am using both KFAS
and dlm (but since are both gpl I can't look at the code to see what is
going on.

I suggestion: in many applications is useful to have the variance
covariance of the prediction. It should be easy enough to make kfs returns
F the matrix and not only the diagonal.

I am closing this.
On Fri, 11 Dec 2015 at 19:37, Jouni Helske notifications@github.com wrote:

The variable V_mu is not the variance of observation y_t but the variance
of mean i.e. in Gaussian case Var(Z_t alpha_t | Y) so there is no bug.
But to be honest the documentation is bit poor in KFS regarding that, I
have to rewrite it better.

Instead of manipulating the model after running SSModel, you can just add
NAs to observations when constructing it via SSModel:

model_kfas_pred <- SSModel(c(Y, rep(NA, 12)) ~ -1 + SSMcustom(Z = Z,
T = T, R = R, Q = Q, a1 = a1, P1 = P1), H = H)

out <- KFS(model_kfas_pred)
c(out$V_mu) + H

But based on your code (dlmForecast) it looks like you are trying to make
make predictions, so you don't actually need smoothed variances. You can
get the prediction error variances from F component:

out$F[101:112]

And you can also use predict method for SSModel object (I am going to add
similar method to KFS object soon too), for example:

predict(model_kfas, n.ahead = 12, se.fit = TRUE, interval = "prediction")
Time Series:
Start = 101
End = 112
Frequency = 1
fit lwr upr se.fit
101 0.3848579 -2.704122 3.473838 1.218154
102 0.3463721 -3.160791 3.853535 1.483900
103 0.3117349 -3.500678 4.124148 1.668408
104 0.2805614 -3.762244 4.323367 1.804080
105 0.2525053 -3.967708 4.472718 1.906911
106 0.2272547 -4.131369 4.585878 1.986306
107 0.2045293 -4.263064 4.672123 2.048361
108 0.1840763 -4.369872 4.738024 2.097280
109 0.1656687 -4.457044 4.788381 2.136084
110 0.1491018 -4.528570 4.826773 2.167005
111 0.1341916 -4.587527 4.855911 2.191732
112 0.1207725 -4.636326 4.877871 2.211558

—
Reply to this email directly or view it on GitHub
#9 (comment).

[helske]

Yes that is actually already possible with the github version of the package, there is now function mvInnovations which computes the usual multivariate prediction errors and their variances.