DLM results looking wonky

Question

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.

Answer 1

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.

Answer 2

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.