When does initialisation change the results in the (R) KFAS package?

I'm trying to get an idea about the importance of the initialisation for the filtering and smoothing results of KFAS. The code below first replicates the random walk with drift example in the KFAS article with the "exact diffuse initialisation". Then I do the same with the more naive explicit initialisation P1 <- matrix(c(1e7, 0, 0, 1e7), 2, 2). The results coincide almost perfectly (as I would expect in a simple example like this). Everything fine so far.

Now I change the dataset. ("tdx" is a series of tax rates that I want to analyse using state space models.) The version with P1 <- matrix(0, 2, 2); P1inf <- diag(2) gives plausible filtering and smoothing results. But the version with P1 <- matrix(c(1e7, 0, 0, 1e7), 2, 2) looks like complete nonsense. It does not coincide at all with the reference version and shows a diverging trend.

Can someone explain to me what is going on here?

(Apparently, I'm not able to attach the data file here. It's not large (approximately 100 data points). I would be happy to send it to anyone wanting to give it a try.)

rm(list = ls())
library(KFAS)
data("alcohol", package = "KFAS")

deaths <- window(alcohol[, 2], end = 2007)
population <- window(alcohol[, 6], end = 2007)
data <- deaths / population

# load("data_tdx.RData")
# data <- tdx

# this is the model from the KFAS article (random walk with drift)
Zt <- matrix(c(1, 0), 1, 2)
Ht <- matrix(NA)
Tt <- matrix(c(1, 0, 1, 1), 2, 2)
Rt <- matrix(c(1, 0), 2, 1)
Qt <- matrix(NA)
a1 <- matrix(c(1, 0), 2, 1)
P1 <- matrix(0, 2, 2)
P1inf <- diag(2)
model_gaussian <- SSModel(data ~ -1 +
  SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1,
  P1inf = P1inf), H = Ht)

fit_gaussian <- fitSSM(model_gaussian, inits = c(0, 0), method = "BFGS")
out_gaussian <- KFS(fit_gaussian$model)

# this reproduces figure 1 from the article
temp <- cbind(data, out_gaussian$a[,1], out_gaussian$alphahat[,1])
cols <- c("black","red","blue")
lwd <- 2
lty <- c(1,1,1)
plot.ts(temp, plot.type="single", col = cols,lwd = lwd,lty=lty,xlab="",ylab="",main="")
legend("bottomright",
       legend = c("Observation Data","Trend filtered","Trend smoothed"),
       col = c("black","red","blue"), lty = c(1,1,1),
       lwd = c(2,2,2), bty="n", cex=0.9 )


# the same with explicit initialisation
P1 <- matrix(c(1e7, 0, 0, 1e7), 2, 2)
model_gaussian1 <- SSModel(data ~ -1 +
  SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1), H = Ht)

fit_gaussian1 <- fitSSM(model_gaussian1, inits = c(0, 0), method = "BFGS")
out_gaussian1 <- KFS(fit_gaussian1$model)

temp <- cbind(data, out_gaussian1$a[,1], out_gaussian1$alphahat[,1])
cols <- c("black","red","blue")
lwd <- 2
lty <- c(1,1,1)
plot.ts(temp, plot.type="single", col = cols,lwd = lwd,lty=lty,xlab="",ylab="",main="")
legend("bottomright",
       legend = c("Observation Data","Trend filtered","Trend smoothed"),
       col = c("black","red","blue"), lty = c(1,1,1),
       lwd = c(2,2,2), bty="n", cex=0.9 )

(check <- cbind(out_gaussian$a[,1], out_gaussian1$a[,1],
               out_gaussian$alphahat[,1], out_gaussian1$alphahat[,1]))

I ran the code with the alternative data set. I expected the results of the two models to be very close and deviations to be concentrated at the beginning of the series (initialisation). Actually, the two models differ fundamentally and the second one cannot be a plausible filtering outcome.

Having run the code block with the tax rate data set you provided via email, it seems the discrepancy in result comes down to the optimization routine getting stuck in some local optimum. You can see this by comparing the function values and parameters produced via fit_gaussian$optim.out and fit_gaussian1$optim.out.

This can often be an issue in state space models, as Helske comments in the KFAS vignette/JSS paper:
>Parameter estimation of a state space model is often a difficult task, as the likelihood surface contains multiple maxima, thus making the optimization problem highly dependent on the initial values. Often the unknown parameters are related to the unobserved latent states, such as the covariance matrix in this example, with little a priori knowledge. Therefore, it is challenging to guess good initial values, especially in more complex settings. Thus, multiple initial value configurations possibly with several different type of optimization routines is recommended before one can be reasonably sure that proper optimum is found.

I found several possible solutions to this issue. The first is to reduce the P1 covariance initialization to a smaller value such as 1e4 (which makes sense--you know the latent state will have low variance since the tax rate data has very low variance). Another option is to have a better initialization for the optimization, as Helske suggests. For example you could set it as the variance of your data (or rather the log of the variance since that is how the fitSSM function is designed). This will lead to the intended result:

fit_gaussian1 <- fitSSM(model_gaussian1, inits = c(-8, -8), method="BFGS")

Finally, you could follow the other suggestion of Helske and change the optimization method:

fit_gaussian1 <- fitSSM(model_gaussian1, inits = c(0, 0))

This uses the default Nelder-Mead method of the optim function, and produces more sensible results in this case (although that may not be generally true).