# DLM results looking wonky

I am teaching myself DLM's using R's dlm package and have two strange results. I am modeling a time series using three combined elements: a trend (dlmModPoly), seasonality (dlmModTrig), and moving seasonality (dlmModReg).

The first strange result is with the $f (one-step-ahead foreacast) result. Most of this forecast appears to be one month behind the actual data, which I believe I've seen in many examples of one-step-ahead forecasting online and in books. The strange thing is that the moving seasonality is NOT similarly lagged, but hits exactly where it should. Is this normal? If I use the result's $m to manually assemble the componenet, everything lines up perfectly, so it's weird, though it makes sense in a way: the moving seasonality has exogenous data to help it while the rest of the forecast does not. (Still, it'd be nice to simply lag the resulting $f and see a nice match.) More troubling is the difference I see if I change the degree of dlmModPoly's polynomial (from 1 to 2) in an attempt to get a smoother level. This introduces a huge spike in all three components at month 9. The spikes all basically cancel out in the composite, but obviously make each piece, say the level or the seasonality, look rather ridiculous there. Is this just one of those things that happens and I should be prepared to throw away the result's first year of data as "break-in"? Or is it an indication that something is wrong? (Even in the degree 1 polynomial case, the first year's moving seasonality's level is a bit unsettled, but no huge spike as when I use a degree 2 polynomial.) Here is my R code: lvl0 <- log (my.data[1]) slp0 <- mean (diff (log (my.data))) buildPTR2 <- function (x) { pm <- dlmModPoly (order=1, dV=exp (x[1]), dW=exp (x[2]), m0=lvl0) tm <- dlmModTrig (s=12, dV=exp (x[1]), q=2, dW=exp (x[3:4])) rm <- dlmModReg (moving.season, dV=exp (x[1])) ptrm <- pm + tm + rm return (ptrm) } mlptr2 <- dlmMLE (log (my.data), rep (1, 6), buildPTR2) dptr2 <- buildPTR2 (mlptr2$par)
dptrf2 <- dlmFilter (log (my.data), dptr2)

tsdiag (dptrf2)

buildPTR3 <- function (x)
{
pm <- dlmModPoly (order=2, dV=exp (x[1]), dW=c(0, exp (x[2])), m0=c(lvl0, slp0))
tm <- dlmModTrig (s=12, dV=exp (x[1]), q=2, dW=exp (x[3:4]))
rm <- dlmModReg (moving.season, dV=exp (x[1]))

ptrm <- pm + tm + rm
return (ptrm)
}

mlptr3 <- dlmMLE (log (my.data), rep (1, 8), buildPTR3)
dptr3 <- buildPTR3 (mlptr3\$par)
dptrf3 <- dlmFilter (log (my.data), dptr3)


Per the follow-on question: the data itself is monthly data for 10 years, with each month being the weekly average attendance at a theatrical production. The data definitely has seasonal and moving seasonal effects. I want to model the trend and the seasonal effects to give the management some insight, and to prepare for forecasting. (Which is not directly possible with dlm when you include a dlmModReg component, though that's the next step.)

(I am trying to use an order=2 polynomial component that I believe creates an IRW trend, which is supposed to be nicely smooth.)

If it matters, my moving seasonality is a yearly Big Bash Gala event that can fall in two different months, and I indicate it with 0 for most months and 1 for months in which the Big Bash falls.

-
As a side note, I'll say that the DLM approach to time-series modeling is nice, but it reminds me of C programming: you can do anything, but you can cut your hand off easily as well. The dlm package makes it easy to assemble your matrices, but you've got to figure out what each resulting column means and if you change something, you can easily get crazy results because the columns changed. –  Wayne Dec 14 '10 at 14:08
Could you please add what are you trying to model? –  mpiktas Dec 14 '10 at 14:35
@mpiktas: If you mean what the data is, it's the weekly average attendance (by month) at a theater. If you mean what I want to learn, it's the overall trend of attendance and the seasonal effects, to help the management understand those dynamics and in preparation for forecasting. Hope that helps. –  Wayne Dec 14 '10 at 15:44

Difficult to diagnose without looking at the results, but here is a wild guess: if in the fit the variance associated to the level component is large, that component will "follow" your data. The filtered estimates will nearly coincide with observations, and the forecasts will appear to lag them --which, as I understand, is what you observe--.

I have no explanation for what you mention in your third paragraph.

I would think that the behaviour you note in the fourth paragraph is not so abnormal: you have monthly data and the Kalman filter has to settle down. For a transient period (at least a year) you may look at a wandering state.

My understanding is that the local linear trend (dlmModPoly(order=2,...)) needs not give a smoother fit than the local level model (dlmModPoly(order=1,...)). It all depends on the values of the fitted variances of the level and slope components of the state.

-
(Sorry I could not include graphs, I don't have permission to share the data.) So there can be a Kalman filter settling time. That makes total sense, looking at the Big Bash (moving seasonality) effect: the first year has movement in the "baseline": in the one case it wanders a bit, and in the second case it has that horrible spike. After that, the baseline (Big Bash indicator 0) runs nearly flat around 0.95, while the spikes (Big Bash indicator 1) go up to 1.10-1.25 (multiplicative effect), which is what I'd expect. The 0.95 baseline I think due to me not centering the 0/1 Big Bash series. –  Wayne Dec 15 '10 at 22:46
In my previous note I used "spike" in two different ways: 1) the Bad Spike that appeared evidently due to filter settling, and 2) the Good Spikes which correspond to the once-a-year Big Bash indicator being 1. Sorry for any confusion. –  Wayne Dec 15 '10 at 22:57
@F Tusell: several good hints in one reply, thanks! –  Wayne Dec 15 '10 at 23:05

So you have monthly data with trend and seasonality and you want to both analyse the trend/seasonal components and produce forecasts. These are two separate tasks. While you can do both with dlm, there are simpler approaches if you separate the tasks.

For studying the trend and seasonality, I suggest using STL via the stl() function in R. It is robust and has nice graphics that are easy to explain to non-statisticians.

For forecasting, I would use a seasonal ARIMA model. To start with, ignore the Big Bash Gala event, and try fitting a seasonal ARIMA model using auto.arima() from the forecast package. The automatic seasonal differencing selection is not very good, so I would use

fit1 <- auto.arima(x,D=1)


where x is your data.

Then, to add in the BBG effect, set up a dummy regression variable (z) taking value 1 when BBG occurs in the month and 0 otherwise. Add that in to your model using

fit2 <- auto.arima(x,xreg=z,D=1)


Produce forecasts from each model using the forecast() function. For fit2, you will need future values of z.

Even if you want to stick with dlm, the above approaches will provide useful comparative benchmarks.

-
Yes, point well-taken. I was doing things the hard way to learn DLM's, mainly, since one of the touted benefits of the SS/DLM approach is that you model each element explicitly so your results are more open to understanding/explanation than a more black-box approach like ARIMA. In some ways true, and in some ways not. Especially if a combination of tools (stl and arima) can give you similar results with less pain and on more familiar territory. –  Wayne Dec 15 '10 at 22:55
(Not to mention that the dlm package won't forecast if you include a dlmModReg element, in the current version anyhow. As a follow-on exercise, I am going to see if I can modify their dlmForecast function (which is completely in R, thankfully) to handle this specific case.) –  Wayne Dec 15 '10 at 23:00