Boostime: Combinando ARIMA con Catboost

Introducción a Boostime con un ejemplo práctico mediante ARIMA y Catboost.

¿Qué es Boostime?

Boostime es un paquete de R que surge para aplicar ciertos métodos sobre los residuos de un modelo para tratar de aumentar la precisión del mismo. Podríamos resumirlo de la siguiente manera: predicciones = modelo1(serie original) + modelo2(residuos)

Para comprender este enfoque, son necesarias algunas nociones básicas sobre series temporales. En primer lugar, debemos ver la serie como una estructura aditiva o multiplicativa de componentes:

$ Y_{t} = T_{t} + S_{t} + C_{t} + I_{t}$

o bien:

$ Y_{t} = T_{t} S_{t} C_{t} I_{t}$

donde el valor de la serie en el tiempo t vendrá dado por la adición (o multiplicación) de las componentes de tendencia, estacionalidad, ciclicidad e irregularidad. Aunque no es objeto de esta publicación entrar en detalle a explicar cada una de estas componentes, sí es necesario remarcar que los tres primeros son patrones estructurales de los datos, mientras que el último son patrones no estructurales.

Antes de continuar, debemos explicar brevemente también el concepto de ruido blanco.

¿Qué es el ruido blanco?

El ruido blanco es una serie de variables aleatorias no correladas entre sí, es decir, la relación entre las observaciones es aleatoria. Normalmente, se asume que el ruido blanco son variables aleatorias independientes que provienen de una distribución aleatoria identicamente distribuida (i.i.d) con media 0 y varianza

$\sigma^{2}$. Veamos un ejemplo de ruido blanco:

1
2
3
set.seed(123)
ruido_blanco <- rnorm(500, 0, 1)
plot.ts(ruido_blanco, col = "red", ylab = "", main = "Ruido Blanco")

Además del método visual, una forma de corroborar si una serie es ruido blanco es a través del gráfico de autocorrelación, donde deberemos comprobar que no existe correlación entre la serie y sus lags:

1
TSstudio::ts_cor(ts(ruido_blanco), lag.max = 24, type = "acf")

Como vemos en el gráfico superior, excepto el lag 12 están contenidos dentro del nivel de significancia por lo que se ve claramente que estamos ante una serie de ruido blanco. Otro método para ver que no existe correlación entre la serie y sus lags podría ser utilizar el test de Ljung-Box. La hipótesis nula de dicho test es que no existe correlación entre la serie y los lags, por lo que p-valor menor a 0.05 (para un nivel de significancia del 95%) indicaría la existencia de correlación entre las series.

1
2
3
4
5
6
7
8
p_values <- c()

for(i in 1:24){
  p_values[i] <- Box.test(ruido_blanco, lag = i, type = "Ljung-Box")$p.value
}

plot(p_values, ylim = c(0, 1))
abline(h = 0.05, col = "red")

En ambas gráficas podemos ver claramente como se confirma lo esperado, que el ruido blanco efectivamente lo es!

¿Qué pasa con los residuos?

En boostime contamos con dos algoritmos principales (por el momento) para modelar la serie, Arima y Prophet. Estos modelos serán los encargados en primera instancia de intentar capturar los patrones estructurales de la serie (tendencia, estacionalidad y ciclicidad). En un escenario ideal, los residuos de estos modelos (recuerda, la diferencia entre el valor observado en la serie y el valor predicho por el modelo) deberían ser ruido blanco. Sin embargo, si observamos correlación en los residuos es un indicador de que el modelo no ha sido capaz de capturar todos los patrones estructurales, y es aquí donde entra en juego boostime. La idea es extraer esos patrones estructurales que el primer algoritmo no ha sido capaz de capturar mediante un segundo algoritmo, en este caso, (y por el momento), mediante Catboost o LightGBM. Los residuos de este segundo algoritmo deberían ser ruido blanco.

Ejemplo: AirPassengers

La serie que vamos a intentar modelizar se llama “AirPassengers” y contiene información del número de pasajeros de avión desde el año 1949 hasta el año 1960. Veamos el gráfico:

1
plot.ts(AirPassengers, col = "red")

Claramente podemos observar que nos encontramos ante una serie de tipo multiplicativo, como puede apreciarse al aumentar la amplitud año tras año en la componente estacional.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
library(tidyverse)

df <- tibble::tibble(date = zoo::as.Date(time(AirPassengers)),
                     value = as.numeric(AirPassengers),
                     trend = decompose(AirPassengers, type = "multiplicative")$trend %>% as.numeric(),
                     seasonal = decompose(AirPassengers, type = "multiplicative")$seasonal %>% as.numeric())

df <- df %>% 
      timetk::future_frame(.date_var = date, .length_out = "2 years", .bind_data = TRUE) %>%
      dplyr::mutate(month = lubridate::month(date))

train <- df %>% dplyr::filter(!is.na(value))
      
test <- df %>% dplyr::filter(is.na(value)) 

splits <- timetk::time_series_cv(train, assess = "1 year", cumulative = TRUE, slice_limit = 1)

Vamos a ver un gráfico sobre la distribución de los datos para entrenamiento y test:

1
timetk::plot_time_series_cv_plan(splits, .date_var = date, .value = value)

Y vamos a ver una descomposición de la serie:

1
TSstudio::ts_decompose(AirPassengers, type = "multiplicative")

Modelos

Vamos a modelar sólo en función del mes del año en el que nos encontramos para mostrar la diferencia entre usar el modelo Arima únicamente (el cuál tendrá correlación en los residuos) y usar adicionalmente Catboost para modelar los residuos. Veremos la diferencia en el ajuste.

1
library(boostime)
## Loading required package: parsnip

## == Welcome to boostime =========================================================
## If you are interested in time series, maybe you would like to check my other packages: garchmodels and bayesmodels
## </> If you find this package useful, please leave a star: https://github.com/AlbertoAlmuinha/boostime </>
1
2
library(modeltime)
library(tidymodels)
## -- Attaching packages -------------------------------------- tidymodels 0.1.3 --

## v broom        0.7.6      v rsample      0.1.0 
## v dials        0.0.9      v tune         0.1.5 
## v infer        0.5.4      v workflows    0.2.2 
## v modeldata    0.1.0      v workflowsets 0.0.2 
## v recipes      0.1.16     v yardstick    0.0.8

## -- Conflicts ----------------------------------------- tidymodels_conflicts() --
## x scales::discard() masks purrr::discard()
## x dplyr::filter()   masks stats::filter()
## x recipes::fixed()  masks stringr::fixed()
## x dplyr::lag()      masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step()   masks stats::step()
## * Use tidymodels_prefer() to resolve common conflicts.
1
2
3
4
5
splits <- splits$splits[[1]]

arima_model <- arima_reg() %>%
               set_engine("auto_arima") %>%
               fit(value ~ date + month, data = training(splits))
## frequency = 12 observations per 1 year
1
2
3
arima_catboost_model <- boost_arima() %>%
                        set_engine("auto_arima_catboost", verbose = 0) %>%
                        fit(value ~ date + month, data = training(splits))
## frequency = 12 observations per 1 year
1
arima_model$fit$models$model_1 %>% forecast::checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(1,1,1)(0,0,2)[12] errors
## Q* = 118.88, df = 19, p-value = 2.22e-16
## 
## Model df: 5.   Total lags used: 24

Claramente podemos ver que existe correlación entre los residuos del modelo Arima, como podemos ver en el gráfico de autocorrelación. Esto se ve confirmado por el p-valor del test de Ljung-Box. Por lo tanto, podemos concluir que este sería un buen caso para probar un modelo que modele también los residuos (como el segundo que hemos hecho). Vamos a comparar su desempeño con Modeltime.

1
2
3
4
modeltime_tbl <- modeltime_table(arima_model,
                                 arima_catboost_model)

modeltime_tbl
## # Modeltime Table
## # A tibble: 2 x 3
##   .model_id .model   .model_desc                                   
##       <int> <list>   <chr>                                         
## 1         1 <fit[+]> REGRESSION WITH ARIMA(1,1,1)(0,0,2)[12] ERRORS
## 2         2 <fit[+]> ARIMA(1,1,0)(0,1,0)[12] W/ CATBOOST ERRORS

Quizás en este punto te surja la pregunta al observar la tabla de porque los dos modelos arima detectados con el algoritmo auto.arima no tienen el mismo orden si a priori hemos introducido los mismos datos. La respuesta está en que en el segundo caso el modelo arima se entrena sin regresores externos para así remover sólo la tendencia y dejar la tarea entera de la estacionalidad para el segundo modelo (en este caso, Catboost, quién utilizará los regresores externos). En nuestra experiencia de esta manera se obtendrán resultados mejores en la mayoría de casuísticas.

1
2
3
4
5
modeltime_tbl %>% 
  modeltime_calibrate(new_data = testing(splits)) %>%
  modeltime_forecast(new_data = testing(splits),
                     actual_data = training(splits)) %>%
  plot_modeltime_forecast()

Claramente podemos observar como el modelo que combina Arima + Catboost ha conseguido captar la tendencia y el aumento de la amplitud en la estacionalidad mientras que el modelo Arima no ha sido capaz de modelar bien la estacionalidad. Veamos algunas métricas:

1
2
3
modeltime_tbl %>% 
  modeltime_calibrate(new_data = testing(splits)) %>%
  modeltime_accuracy()
## # A tibble: 2 x 9
##   .model_id .model_desc                .type   mae  mape  mase smape  rmse   rsq
##       <int> <chr>                      <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1         1 REGRESSION WITH ARIMA(1,1~ Test   28.9  5.81 0.598  5.85  36.6 0.922
## 2         2 ARIMA(1,1,0)(0,1,0)[12] W~ Test   18.0  4.09 0.373  3.93  24.2 0.954

Finalmente, reentrenamos sobre la serie entera y predecimos los dos siguientes años:

1
2
3
4
5
modeltime_tbl %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_refit(data = df) %>%
  modeltime_forecast(new_data = test, actual_data = train) %>%
  plot_modeltime_forecast()
## frequency = 12 observations per 1 year
## frequency = 12 observations per 1 year

Como podemos observar, la combinación de Arima y Catboost sí consigue hacer un buen trabajo y modelar correctamente la tendencia y la estacionalidad de la serie mientras que el modelo Arima inicial no es capaz de modelar la estacionalidad correctamente.

Si te ha gustado este post y el paquete boostime, puedes mostrar tu apoyo dandonos una estrella en nuestro repositorio de GitHub! :https://github.com/AlbertoAlmuinha/boostime

updatedupdated2021-08-192021-08-19
Load Comments?