nelsonauner - 7 months ago 39
R Question

# Computational instability in R Forecast package? Edit: Fixed

## Original Question:

I have the following time series data observed daily:

``````series <- c(10, 25, 8, 27, 18, 21, 12, 9, 31, 18, 8, 30, 14, 13, 10, 14,
14, 14, 6, 9, 22, 21, 22, 8, 7, 6, 22, 21, 36, 16, 2, 13, 23,
40, 12, 27, 18, 10, 11, 37, 44, 30, 40, 25, 13, 11, 58, 56, 46,
39, 28, 27, 19, 20, 97, 90, 70, 73, 30, 22, 97, 34)
``````

and want to fit it using
`tbats`
from the R
`forecasts`
package. I also want to model it with weekly correlation:

`````` library(forecast)
x.msts = msts(series,seasonal.periods = 7)
model <- tbats(x.msts)
``````

Examing/plotting the model with
`str`
reveals a huge fitted variance of
`4.9e+17`
.

And, plotting the forecast going forward, we observe massive swings:

``````> forecast(model)\$mean

Multi-Seasonal Time Series:
Start: 9 7
Seasonal Periods: 7
Data:
[1]  1.483789e+44 -1.399297e+42 -2.566455e+44 -1.374316e+43 -1.527758e+38
[6]  2.036194e+42  5.639596e+42  8.231600e+40 -2.578859e+41 -1.355840e+43
``````

Are these estimates the "correct" solution to the TBATS model fitting procedure, or is there a bug in the
`forecast`
package? If not a bug, can someone help me understanding mathematically why this normal-looking time series produces these estimates?

This is my first post on CV so apologies if this should be on SO!

I have filed a bug report on github

Also some people have noticed that I'm not using multiple seasonality factors, so I want to show here that the bug is still an issue:

``````x2.msts <- msts(series,seasonal.periods = c(7,30))
model_x2_1 <- tbats(x2.msts) # high variance
model_x2_2 <- tbats( series, seasonal.periods = c(7,30) ) # also high variance
``````

This is perhaps the same problem as described here, so the reason is presumably a bug in the forecast package. I'm not sure if the following alternative will give you the desired result, but you can leave `series` as is and put `seasonal.periods=7` in the call of `tbats`:

``````library(forecast)

series <- c(10, 25, 8, 27, 18, 21, 12, 9, 31, 18, 8, 30, 14, 13, 10, 14,
14, 14, 6, 9, 22, 21, 22, 8, 7, 6, 22, 21, 36, 16, 2, 13, 23,
40, 12, 27, 18, 10, 11, 37, 44, 30, 40, 25, 13, 11, 58, 56, 46,
39, 28, 27, 19, 20, 97, 90, 70, 73, 30, 22, 97, 34)

x.msts <- msts(series,seasonal.periods = 7)
model_1 <- tbats(x.msts)

model_2 <- tbats( series, seasonal.periods = 7 )
``````

The variance of `model_2` is much better than that of `model_1`:

``````> str(model_1)
List of 19
\$ lambda           : num 0.21
\$ alpha            : num 0.374
\$ beta             : NULL
\$ damping.parameter: NULL
\$ gamma.values     : NULL
\$ ar.coefficients  : num [1:2] 1.296 -0.911
\$ ma.coefficients  : num [1:2] -1.62 0.98
\$ likelihood       : num 549
\$ optim.return.code: int 0
\$ variance         : num 4.9e+17
\$ AIC              : num 571
\$ parameters       :List of 2
..\$ vect   : num [1:6] 0.21 0.374 1.296 -0.911 -1.615 ...
..\$ control:List of 6
.. ..\$ use.beta    : logi FALSE
.. ..\$ use.box.cox : logi TRUE
.. ..\$ use.damping : logi FALSE
.. ..\$ length.gamma: num 0
.. ..\$ p           : int 2
.. ..\$ q           : int 2
\$ seed.states      : num [1:5, 1] 4.16 0 0 0 0
\$ fitted.values    : Time-Series [1:62] from 1 to 9.71: 19.97 19.28 4.53 21.83 56.15 ...
..- attr(*, "msts")= num 7
\$ errors           : Time-Series [1:62] from 1 to 9.71: -1.206 0.496 0.828 0.415 -2.354 ...
..- attr(*, "msts")= num 7
\$ x                : num [1:5, 1:62] 3.71 -1.21 0 -1.21 0 ...
\$ seasonal.periods : NULL
\$ y                : Time-Series [1:62] from 1 to 9.71: 10 25 8 27 18 21 12 9 31 18 ...
..- attr(*, "msts")= num 7
\$ call             : language tbats(y = x.msts)
- attr(*, "class")= chr "bats"
>
``````

.

``````> str(model_2)
List of 23
\$ lambda           : num 0.198
\$ alpha            : num 0.198
\$ beta             : NULL
\$ damping.parameter: NULL
\$ gamma.one.values : num -0.0157
\$ gamma.two.values : num 0.00991
\$ ar.coefficients  : NULL
\$ ma.coefficients  : NULL
\$ likelihood       : num 553
\$ optim.return.code: int 0
\$ variance         : num 0.969
\$ AIC              : num 571
\$ parameters       :List of 2
..\$ vect   : num [1:4] 0.19842 0.19782 -0.0157 0.00991
..\$ control:List of 6
.. ..\$ use.beta    : logi FALSE
.. ..\$ use.box.cox : logi TRUE
.. ..\$ use.damping : logi FALSE
.. ..\$ length.gamma: int 2
.. ..\$ p           : num 0
.. ..\$ q           : num 0
\$ seed.states      : num [1:5, 1] 4.1851 0.3176 0.0103 -0.5806 0.4447
\$ fitted.values    : Time-Series [1:62] from 1 to 62: 25.1 20 11.1 10.2 24.3 ...
\$ errors           : Time-Series [1:62] from 1 to 62: -1.594 0.41 -0.507 1.697 -0.552 ...
\$ x                : num [1:5, 1:62] 3.87 -0.231 0.456 -0.626 -0.125 ...
\$ seasonal.periods : num 7
\$ k.vector         : int 2
\$ y                : Time-Series [1:62] from 1 to 62: 10 25 8 27 18 21 12 9 31 18 ...
\$ p                : num 0
\$ q                : num 0
\$ call             : language tbats(y = series, seasonal.periods = 7)
- attr(*, "class")= chr [1:2] "tbats" "bats"
>
``````