helske / KFAS

KFAS: R Package for Exponential Family State Space Models
64 stars 17 forks source link

Fitting a Multivariate Local Linear Trend model #46

Closed Vera1996 closed 3 years ago

Vera1996 commented 4 years ago

Dear Helske,

I'm trying to fit a Local Linear Trend model but unfortunately I do not understand the update function and the initial values in the fitSSM line. Can you please help me?

covid_new <- ts(data = covid_new, frequency = 1, start = c(1)) ts.plot(window(covid_new[, 1:4]), col = 1:4, ylab = "Log first differences confirmed cases", xlab = "Time") legend("topright",col = 1:4, lty = 1, legend = colnames(covid_new)[1:4])

covid_new <- window(covid_new, start = 1, end = 145) model <- SSModel(covid_new[, 1:4] ~ SSMtrend(2, Q = list(matrix(NA, 4, 4), matrix(NA, 4, 4))), H = matrix(NA, 4, 4))

updatefn <- function(pars, model, ...) { Q <- diag(exp(pars[1:4])) Q[upper.tri(Q)] <- pars[5:8] model["Q", etas = "level"] <- crossprod(Q) Q <- diag(exp(pars[9:12])) Q[upper.tri(Q)] <- pars[13:16] model["Q", etas = "slope"] <- crossprod(Q) model }

init <- chol(cov(covid_new[, 1:4]) / 10) fitinit <- fitSSM(model, updatefn = updatefn, inits = rep(c(log(diag(init)), init[upper.tri(init)], init[upper.tri(init)]), 3), method = "BFGS") -fitinit$optim.out$val

fit <- fitSSM(model, updatefn = updatefn, inits = fitinit$optim.out$par, method = "BFGS", nsim = 250) -fit$optim.out$val

varcor <- fit$model["Q", etas = "level"] varcor[upper.tri(varcor)] <- cov2cor(varcor)[upper.tri(varcor)] print(varcor, digits = 2)

varcor <- fit$model["Q", etas = "slope"] varcor[upper.tri(varcor)] <- cov2cor(varcor)[upper.tri(varcor)] print(varcor, digits = 2)

(out <- KFS(fit$model, nsim = 1000))

Kind regards,

Vera Vial

helske commented 4 years ago

What is it exactly that you do not understand? updatefn is a function that updates your model object according to the argument pars, and the initial values for that argument are given to the fitSSM using the argument inits.

The second run of fitSSM is not needed above as you seem to have fully gaussian model so argument nsim is ignored

Vera1996 commented 4 years ago

The update function above, isn't right. I need to write an update function for a Multivariate Local Linear Trend model but I don't know how because of the two Q matrices. Because of the two Q matrices, I don't know what input I need for the initial values. Can you help me with the update function of a Multivariate Local Linear Trend model and the initial values?

helske commented 4 years ago

There is only one Q even though you define the model blockwise. See model$Q .

Vera1996 commented 4 years ago

Oké yes. I see. But how do you define the update function with the covariances? I really can’t figure it out how to formulate the update function..

helske commented 4 years ago

You can define the update function in any way you want, the only requirements for it are that it takes the same arguments as in the above example and that it returns an SSModel object. The structure of matrix Q of your model should match your model definition, but if it feels difficult to create it from the scratch you can also just call the model building function SSModel (or SSMtrend-part) inside updatefn. The parameterization in the example is based on the common cholesky parameterization of Q, which takes care of the positive-definite requirement of covariance matrices, and thus you only need to define the diagonal and the upper diagonal parts of the cholesky factor of Q. See for example http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.31.494&rep=rep1&type=pdf