danielkv - 8 months ago 79

R Question

I am using 'KFAS' package from R to estimate a state-space model with the Kalman filter. My measurement and transition equations are:

y_t = Z_t * x_t + \eps_t (measurement)

x_t = T_t * x_{t-1} + R_t * \eta_t (transition),

with \eps_t ~ N(0,H_t) and \eta_t ~ N(0,Q_t).

So, I want to estimate the variances H_t and Q_t, but also T_t, the AR(1) coefficient. My code is as follows:

`library(KFAS)`

set.seed(100)

eps <- rt(200, 4, 1)

meas <- as.matrix((arima.sim(n=200, list(ar=0.6), innov = rnorm(200)*sqrt(0.5)) + eps),

ncol=1)

Zt <- 1

Ht <- matrix(NA)

Tt <- matrix(NA)

Rt <- 1

Qt <- matrix(NA)

ss_model <- SSModel(meas ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt,

Q = Qt), H = Ht)

fit <- fitSSM(ss_model, inits = c(0,0.6,0), method = 'L-BFGS-B')

But it returns:

The NA definitions for the variances works well, as documented in the package's paper. However, it seems this cannot be done for the AR coefficients. Does anyone know how can I do this?

Note that I am aware of the SSMarima function, which eases the definition of the transition equation as ARIMA models. Although I am able to estimate the AR(1) coef. and Q_t this way, I still cannot estimate the \eps_t variance (H_t). Moreover, I am migrating my Kalman filter codes from EViews to R, so I need to learn SSMcustom for other models that are more complicated.

Thanks!

Answer

It seems that you are missing something in your example, as your error message comes from the function `fitSSM`

. If you want to use `fitSSM`

for estimating general state space models, you need to provide your own model updating function. The default behaviour can only handle NA's in covariance matrices H and Q. The main goal of `fitSSM`

is just to get started with simple stuff. For complex models and/or large data, I would recommend using your self-written objective function (with help of `logLik`

method) and your favourite numerical optimization routines manually for maximum performance. Something like this:

```
library(KFAS)
set.seed(100)
eps <- rt(200, 4, 1)
meas <- as.matrix((arima.sim(n=200, list(ar=0.6), innov = rnorm(200)*sqrt(0.5)) + eps),
ncol=1)
Zt <- 1
Ht <- matrix(NA)
Tt <- matrix(NA)
Rt <- 1
Qt <- matrix(NA)
ss_model <- SSModel(meas ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt,
Q = Qt), H = Ht)
objf <- function(pars, model, estimate = TRUE) {
model$H[1] <- pars[1]
model$T[1] <- pars[2]
model$Q[1] <- pars[3]
if (estimate) {
-logLik(model)
} else {
model
}
}
opt <- optim(c(1, 0.5, 1), objf, method = "L-BFGS-B",
lower = c(0, -0.99, 0), upper = c(100, 0.99, 100), model = ss_model)
ss_model_opt <- objf(opt$par, ss_model, estimate = FALSE)
```

Same with `fitSSM`

:

```
updatefn <- function(pars, model) {
model$H[1] <- pars[1]
model$T[1] <- pars[2]
model$Q[1] <- pars[3]
model
}
fit <- fitSSM(ss_model, c(1, 0.5, 1), updatefn, method = "L-BFGS-B",
lower = c(0, -0.99, 0), upper = c(100, 0.99, 100))
identical(ss_model_opt, fit$model)
```

Source (Stackoverflow)