[R] Strange behaviour of `dlm` package

[zuba]

Dear R-help!

I have encountered strange behaviour (that is, far-off filtering, smoothing
and forecast distributions under certain conditions) in the `dlm` package by
Giovanni Petris.

Here is an example:

I use the annual hotel bookings time series data, which I model using a
second order polinomial DLM.

First I perform the analysis with the data in logarithmic form and
everything seems to be ok.

I repeat the analysis with the original, non logarithmic data and the
filtering, smoothing and forecast distributions are *very* far off the
actual data. 

I again repeat the analysis, this time I also specify the expected value and
variance of the pre-sample state vector (m0 and C0). The filtering and
forecast distribution seem to be influenced by these values only and ignore
the actual data. 

I would like to know what I would have to do in order to get plausible
results from DLM with data in non-logarithmic form. Is there something that
I am missing?

Thanks in advance!

Martin Zuba


References:
dlm introductory paper
<http://cran.r-project.org/web/packages/dlm/vignettes/dlm.pdf>   by Giovanni
Petris (2009)
dlm package  help file <http://benz.nchu.edu.tw/~finmyc/dlm.pdf>  


This R code reproduces my analysis:



# DLM anomaly

rm(list=ls())

hotann <- ts(start=1973, data=c(13733531, 15445493, 16024342, 16715224,
16463957, 17896118, 18263114, 
                                18727641, 19756293, 20352440, 22489054,
23442393, 24170736, 25308044, 
                                25021560, 25724331, 25990336, 26602038,
27292250, 28891954, 29659700, 
                                31533579, 32513666, 33628559, 34494451))

plot(hotann, ylab="Annual hotel bookings")

# Analysis with log data

tsdata <- log(hotann)

buildfun <- function (x) {
  dlmModPoly(order = 2, dV = exp(x[1]), dW = c(0,exp(x[2])))
}

fit <- dlmMLE(y=tsdata, parm=c(0,0), build=buildfun)
# Warning: a numerically singular 'V' has been slightly perturbed to make it
nonsingular
fit$conv

dlmTsdata <- buildfun(fit$par)
tsdataFilter <- dlmFilter(tsdata, mod=dlmTsdata)
tsdataSmooth <- dlmSmooth(tsdata, mod=dlmTsdata)

plot(tsdata, lwd=2)
for (i in 1:10)
  lines(lty=6, col="blue", dropFirst(dlmBSample(tsdataFilter))[,1])
# looks ok!


tsdataForecast <- dlmForecast(tsdataFilter, nAhead=20)
sqrtR <- sapply(tsdataForecast$R, function(x) sqrt(x[1,1]))
pl <- tsdataForecast$a[,1] + qnorm(0.05, sd= sqrtR)
pu <- tsdataForecast$a[,1] + qnorm(0.95, sd= sqrtR)

x <- ts.union(tsdata,tsdataSmooth$s[,1],tsdataForecast$a[,1],pl,pu)
plot(x, plot.type="single", type="l",
col=c("black","blue","red","orange","orange"), ylab="TS data")
# looks ok!

# Repeat with non-log data.

tsdata <- hotann

buildfun <- function (x) {
  dlmModPoly(order = 2, dV = exp(x[1]), dW = c(0,exp(x[2])))
}

fit <- dlmMLE(y=tsdata, parm=c(0,0), build=buildfun)
fit$conv

dlmTsdata <- buildfun(fit$par)
tsdataFilter <- dlmFilter(tsdata, mod=dlmTsdata)
tsdataSmooth <- dlmSmooth(tsdata, mod=dlmTsdata)

plot(tsdata, lwd=2, ylim=c(0, 6e+07))
for (i in 1:10)
  lines(lty=6, col="blue", dropFirst(dlmBSample(tsdataFilter))[,1])
# where is the filtering distribution???


tsdataForecast <- dlmForecast(tsdataFilter, nAhead=20)
sqrtR <- sapply(tsdataForecast$R, function(x) sqrt(x[1,1]))
pl <- tsdataForecast$a[,1] + qnorm(0.05, sd= sqrtR)
pu <- tsdataForecast$a[,1] + qnorm(0.95, sd= sqrtR)

x <- ts.union(tsdata,tsdataSmooth$s[,1],tsdataForecast$a[,1],pl,pu)
plot(x, plot.type="single", type="l",
col=c("black","blue","red","orange","orange"), ylab="TS data")
# where is the forecast???

# strong sensitivity towards expected value (m0) and variance (C0) of
pre-sample state vectors
# m0 is rep(0,order) in the default specification, this might be responsible
for the results above (?)
# repeat analysis, but specify values for m0 and C0

tsdata <- hotann

buildfun <- function (x) {
  dlmModPoly(order = 2, dV = exp(x[1]), dW = c(0,exp(x[2])),
m0=c(x[3],x[4]), C0 = exp(x[5]) * diag(nrow = 2))
}

fit <- dlmMLE(y=tsdata, parm=c(0,0,2e+07,20000,20), build=buildfun)
fit$conv

dlmTsdata <- buildfun(fit$par)
tsdataFilter <- dlmFilter(tsdata, mod=dlmTsdata)
tsdataSmooth <- dlmSmooth(tsdata, mod=dlmTsdata)

plot(tsdata, lwd=2, ylim=c(0, 6e+07))
for (i in 1:10)
  lines(lty=6, col="blue", dropFirst(dlmBSample(tsdataFilter))[,1])
# the distribution only follows m0, C0 and seems ignorant of the data


tsdataForecast <- dlmForecast(tsdataFilter, nAhead=20)
sqrtR <- sapply(tsdataForecast$R, function(x) sqrt(x[1,1]))
pl <- tsdataForecast$a[,1] + qnorm(0.05, sd= sqrtR)
pu <- tsdataForecast$a[,1] + qnorm(0.95, sd= sqrtR)

x <- ts.union(tsdata,tsdataSmooth$s[,1],tsdataForecast$a[,1],pl,pu)
plot(x, plot.type="single", type="l",
col=c("black","blue","red","orange","orange"), ylab="TS data")
# the forecast only follows m0, C0 and seems ignorant of the data



--
View this message in context: http://r.789695.n4.nabble.com/Strange-behaviour-of-dlm-package-tp4683255.html
Sent from the R help mailing list archive at Nabble.com.