danielkv - 1 year ago 198

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!

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

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)
```

Recommended from our users: **Dynamic Network Monitoring from WhatsUp Gold from IPSwitch**. ** Free Download**