EDIT: FIXED - Computational Instability in R Prediction Pack?

Question

The original question:

I have daily daily data:

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 put it using tbats from the R forecasts package. I also want to model it with a weekly correlation:

  library(forecast) x.msts = msts(series,seasonal.periods = 7) model <- tbats(x.msts) # shows "--- loading profile ---" 

Learning / building a model with str shows a huge variance of 4.9e+17 .

And, making a forecast ahead, we observe massive fluctuations:

 > 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 “right” solution for the TBATS model selection procedure or is there an error in the forecast package? If not a mistake, can someone help me understand mathematically why this normal time series gives these estimates?

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

Update after reply:

I posted a github bug report

Also, some people noticed that I do not use several factors of seasonality, so I want to show here that the error is still a problem:

 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 

Answer 1

This is probably the same problem as described here , so the reason seems to be a bug in the forecast package. I'm not sure that the following option will give you the desired result, but you can leave series as it is and put seasonal.periods=7 in the tbats call:

 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 deviation of model_2 much better than the variance 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" >