Las predicciones del modelo BSTS (en R) fallan completamente

15

Después de leer esta publicación de blog sobre los modelos de series de tiempo estructurales bayesianas, quería analizar su implementación en el contexto de un problema para el que había usado ARIMA anteriormente.

Tengo algunos datos con algunos componentes estacionales conocidos (pero ruidosos): definitivamente hay componentes anuales, mensuales y semanales para esto, y también algunos efectos debido a días especiales (como feriados federales o religiosos).

He utilizado el bstspaquete para implementar esto y, por lo que puedo decir, no he hecho nada malo, aunque los componentes y la predicción simplemente no se ven como esperaba. No me queda claro si mi implementación es incorrecta, está incompleta o tiene algún otro problema.

La serie a tiempo completo se ve así:

Datos completos

Puedo entrenar el modelo en algún subconjunto de datos, y el modelo generalmente se ve bien en términos de ajuste (la gráfica está abajo). El código que estoy usando para hacer esto está aquí:

library(bsts)

predict_length = 90
training_cut_date <- '2015-05-01'
test_cut_date <- as.Date(training_cut_date) + predict_length

df = read.csv('input.tsv', sep ='\t')

df$date <- as.Date(as.character(df$date),format="%Y-%m-%d")
df_train = df[df$date < training_cut_date,]

yts <- xts(log10(df_train$count), order.by=df_train$date)

ss <- AddLocalLinearTrend(list(), yts)
ss <- AddSeasonal(ss, yts, nseasons = 7)
ss <- AddSeasonal(ss, yts, nseasons = 12)
ss <- AddNamedHolidays(ss, named.holidays = NamedHolidays(), yts)

model <- bsts(yts, state.specification = ss, niter = 500, seed=2016)

El modelo se ve razonable:

Parcela modelo

Pero si trazo la predicción, en primer lugar, la tendencia es completamente incorrecta, y en segundo lugar, la incertidumbre crece MUY rápidamente, hasta el punto en que no puedo mostrar la banda de incertidumbre en el mismo diagrama que las predicciones sin hacer que el eje y esté en un registro. escala. El código para esta parte está aquí:

burn <- SuggestBurn(0.1, model)
pred <- predict(model, horizon = predict_length, burn = burn, quantiles = c(.025, .975))

La predicción pura se ve así:

predicción pura

Y luego, cuando se reduce a la distribución inicial (con la línea de puntos que muestra la transición del entrenamiento a la predicción, los problemas son obvios:

distribución completa

Intenté agregar más tendencias estacionales, eliminar tendencias estacionales, agregar un término AR, cambiar AddLocalLinearModel a AddGeneralizedLocalLinearTrend y varias otras cosas relacionadas con ajustar el modelo, pero nada ha resuelto los problemas y ha hecho que las predicciones sean más significativas. En algunos casos, la dirección cambia, por lo que, en lugar de caer a 0, la predicción continúa aumentando en función del tiempo. Definitivamente no entiendo por qué el modelo se está rompiendo de esta manera. Cualquier sugerencia sería muy bienvenida.

anthr
fuente
2
¿Por qué no publica sus datos? Intentaré ayudarlo ... No podré responder por qué el modelo se está rompiendo, ya que no uso este enfoque, ya que tiene demasiados supuestos incorporados. precisa sobre cuántos valores se retuvieron, la fecha de inicio y el país de origen.
IrishStat
Muchas gracias por su comentario. He subido los datos en bruto aquí en caso de que tenga tiempo para echar un vistazo. Los datos van desde principios de 2013 hasta finales de este año. También intenté pronosticar con un modelo ARIMA, pero las predicciones a partir de eso tampoco coincidían con los datos de reserva. Los datos retenidos son básicamente solo una fracción de 2015 o 2016, dependiendo de la cantidad de datos de entrenamiento que quiera usar.
anthr
Tengo problemas para descargarlo ... envíe un archivo csv a mi dirección de correo electrónico
IrishStat

Respuestas:

26

Steve Scott aquí. Escribí el paquete bsts. Tengo algunas sugerencias para ti. Primero, tus componentes estacionales no están haciendo lo que crees que están haciendo. Creo que tiene datos diarios, porque está tratando de agregar un componente de 7 temporadas, que debería funcionar correctamente. Pero le ha dicho a su componente estacional anual que se repita cada 12 días. Obtener un componente estacional mensual con datos diarios es un poco difícil de hacer, pero puede hacer una temporada de 52 semanas AddSeasonal(..., nseasons = 52, season.duration = 7).

El seasonal.durationargumento le dice al modelo cuántos puntos de tiempo debería durar cada temporada. El nseasonsargumento le dice cuántas estaciones hay en un ciclo. El número total de puntos de tiempo en un ciclo es season.duration * nseasons.

La segunda sugerencia es que quizás desee pensar en un modelo diferente para la tendencia. El LocalLinearTrendmodelo es muy flexible, pero esta flexibilidad puede aparecer como una variación no deseada en los pronósticos a largo plazo. Hay algunos otros modelos de tendencia que contienen un poco más de estructura. GeneralizedLocalLinearTrend(perdón por el nombre no descriptivo) supone que el componente "pendiente" de la tendencia es un proceso AR1 en lugar de una caminata aleatoria. Es mi opción predeterminada si quiero pronosticar en el futuro. La mayor parte de la variación de su serie temporal parece provenir de la estacionalidad, por lo que puede intentarlo AddLocalLevelo incluso en AddArlugar de hacerlo AddLocalLinearTrend.

Finalmente, en general, si está obteniendo pronósticos extraños y desea averiguar qué parte del modelo tiene la culpa, intente plot(model, "components")ver el modelo descompuesto en las piezas individuales que ha solicitado.

Steve Scott
fuente
FYI: Tengo problemas muy similares con mis datos, que también son diarios. He implementado todas sus sugerencias enumeradas aquí y ninguna parece ayudar.
ZakJ
1
@ Steve Scott Perdón por molestarte Steve, quiero preguntarte esto: si estoy tratando de modelar varias series de tiempo y estoy en un marco de modelo mixto jerárquico, ¿puedo modelar esto usando tu paquete? Por cierto: ¡muchas gracias por tu paquete!
Tommaso Guerrini
4

Creo que también puedes cambiar la grabación predeterminada. Como he usado bsts, creé una cuadrícula de valores de quema y niter con MAPE como mi estadística sobre el período de espera. También intente usar AddStudentLocalLinearTrend en su lugar si sus datos tienen una gran variación para que el modelo espere tal variación

Elderkm2012
fuente
1
Fue útil en mi caso cuando tenía pocos puntos de datos (20)
SCallan