R help - Nov 2014 - Function StructTS: covariance matrix of the initial state vector

Hi everyone:

My question is about the function "StructTS" of the core package
"stats", which fits by maximum likelihood the basic structural time
series model.

According to theory and the references given in "?StructTS", the
covariance matrix of the initial state vector is a diagonal matrix. However, in
"StructTS" it is defined as a matrix containing the same value in the
diagonal and off-diagonal elements of the matrix. The matrix is defined as
follows in "StructTS":

Z$P[] <- 1e+06 * vx

Does anybody know why is this matrix defined as non-diagonal? May it be some
kind of diffuse initialization? That's not the way I understand it in
textbooks though, as it still involves diagonal matrices.

I thought it may be a bug in the code. However, as "StructTS" has been
defined in this way for years and as the documentation claims that the result
for the "AirPassengers" data is superior to the results reported in
other source, I wondered whether there may be a reason to use a non-diagonal
matrix.

I used the package "stsm" [http://cran.r-project.org/package=stsm] to
compare some results based on a non-diagonal and a diagonal initial covariance
matrix.

# original result obtained with StructTS
res0 <- stats::StructTS(log(AirPassengers), type = "BSM")
res0
# Variances:
#     level      slope       seas    epsilon 
# 0.0007719  0.0000000  0.0013969  0.0000000 

# StructTS can be replicated using the interface in package "stsm" as
follows
require("stsm")
mairp <- stsm.model(model = "BSM", y = log(AirPassengers),
transPars = "StructTS")
res1 <- maxlik.td.optim(mairp, KF.version = "KFKSDS",
KF.args = list(P0cov = TRUE), method = "L-BFGS-B")
mairp1 <- set.pars(mairp, pmax(res1$par, .Machine$double.eps))
round(get.pars(mairp1)[c(4,1:3)], 7)
#      var4      var1      var2      var3
# 0.0013969 0.0000000 0.0007719 0.0000000
all.equal(get.pars(mairp1), res0$coef[c(4,1:3)],
tol = 1e-04, check.attributes = FALSE)
# [1] TRUE

# next, the likelihood function and the optimization procedure
# are defined as in StructTS except for Z$P, which is now a diagonal matrix
res2 <- maxlik.td.optim(mairp, KF.version = "KFKSDS",
KF.args = list(P0cov = FALSE), method = "L-BFGS-B")
mairp2 <- set.pars(mairp, pmax(res2$par, .Machine$double.eps))
round(get.pars(mairp2)[c(4,1:3)], 7)
#      var4      var1      var2      var3
# 0.0000642 0.0001293 0.0006995 0.0000000

The results imply a different allocation of the variance across the components,
as reflected by these ratios:

round(get.pars(mairp1)[-4]/get.pars(mairp1)[4], 3)
#  var1  var2  var3
# 0.000 0.553 0.000
round(get.pars(mairp2)[-4]/get.pars(mairp2)[4], 3)
#  var1   var2   var3
# 2.014 10.898  0.000

The estimated seasonal component for the modified "StructTS" version
(i.e., the standard approach using a diagonal matrix) looks better, because the
range is more homogeneous.

ss2 <- char2numeric(mairp2, P0cov = FALSE)
kf2 <- KFKSDS::KF(mairp1 at y, ss2)
ks2 <- KFKSDS::KS(mairp1 at y, ss2, kf2)
par(mfcol = c(2,2), mar = c(2.5,5.1,3,2.1))
plot(tsSmooth(res0)[,1], main = "StructTS", ylab = "level")
plot(tsSmooth(res0)[,3], ylab = "seasonal")
plot(ks2$ahat[,1], main = "Modified StructTS", ylab =
"level")
plot(ks2$ahat[,3], ylab = "seasonal") 

[image in attachment airp1.pdf]

Yet, I am still doubtful whether this is a bug or intentional. In some cases, I
have found it useful to use P0cov=TRUE (non-diagonal matrix) when the likelihood
function is evaluated by the optimization algorithm. Then, given the parameter
estimates, the smoothed components can be extracted setting P0cov=FALSE
(diagonal matrix) in the Kalman smoother, in order to avoid the erratic
behaviour at the beginning of the sample observed in the plot above. See for
example my answer in this post
[http://stats.stackexchange.com/questions/115506/forecasting-a-seasonal-time-series-in-r].


I wonder if this qualifies as a bug or if there is a reason that explains the
non-diagonal matrix used in "StructTS". I would be grateful for some
feedback.

Thanks.

--
javi 
http://jalobe.com
-------------- next part --------------
A non-text attachment was scrubbed...
Name: airp1.pdf
Type: application/pdf
Size: 43816 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20141105/2a470a8a/attachment.pdf>