model0 or model In file na_kalman.R?

[happyzhichengzhao]

Line 200 of file na_kalman.R:
mod <- stats::StructTS(data, ...)$model0
where model0 is the initial state of the filter, not the final state fitted by maximum likelihood .
Should we change model0 to model?

[SteffenMoritz]

Hello @happyzhichengzhao , thanks for opening an issue.
I don't think this has to be changed.

Later in the na_kalman function KalmanSmooth or KalmanRun is called with model0 (and thus log-likelihood is computed)

# 2.2 Selection if KalmanSmooth or KalmanRun
if (smooth == TRUE) {
  kal <- stats::KalmanSmooth(data, mod, nit)
  erg <- kal$smooth # for kalmanSmooth
}
else {
  kal <- stats::KalmanRun(data, mod, nit)
  erg <- kal$states # for kalmanrun
}

The tsSmooth function from the stats package does this in a similar way. Input is a StructTS object and time series data and model0 is given to KalmanSmooth : KalmanSmooth(object$data, object$model0, -1).

Here the whole function:

tsSmooth <- function(object, ...) UseMethod("tsSmooth")

tsSmooth.StructTS <- function(object, ...)
{
    res <- KalmanSmooth(object$data, object$model0, -1)$smooth
    dn <- dim(fitted(object))
    res <- res[, 1L:dn[2L], drop = FALSE]
    dimnames(res) <- dimnames(fitted(object))
    ts(res, start = object$xtsp[1L], frequency = object$xtsp[3L])
}

But you are correct, there is some redundancy in the StructTS source code for our purpose.
We don't need the final state (model) at this point. Yet it is calculated.

When looking at StructTS source code, it is this line, which seems reduntant.
z <- KalmanRun(y, Z, -1, update = TRUE)
But this should not matter too much I think - I just profiled the function and it seems this line does not take too much time. It seems to be the call to optim

res <- optim(init[mask], getLike, method = "L-BFGS-B", lower = rep(0, 
        np + 1L), upper = rep(Inf, np + 1L), control = optim.control)

, which consumes nearly all computing time. But this part is anyway still needed. But I have to profile once again, was not fully satisfied by the results.

Correct me, if you think otherwise or if you just have questions - better to discuss this one time too often, than having an error in the function. Thanks for your help and support for the package ! :)