Pronóstico directo vs. iterado: metodología y validación cruzada

En la siguiente entrada se desarrolla una función de R para proyectar paso a paso en un horizonte determinado. Puntualmente, se comparan dos estrategias: una donde se usa la función auto.arima (ver) para proyectar todos los valores en un horizonte determinado a partir del modelo identificado; y una segunda estrategia que consiste en proyectar a un solo paso, tomar este valor proyectado como observado y posteriormente volver a identificar el modelo ARIMA correspondiente. Posteriormente, se repite este procedimiento tantas veces como sea especificado en la ventana de pronósticos.

 Para el ejemplo utilizaremos la serie temporal AirPassengers:

 library(dplyr)
library(readxl)
library(tibble)
library(ggplot2)
library(forecast)
library(tseries)
require(lubridate)
library(tidyr)
 
autoplot(AirPassengers) +
  theme_minimal()

Ahora definimos el horizonte de pronóstico y usamos la funciones auto.arima y forecast (como hay varios paquetes que utilizan esta función, se indexa dentro del paquete para hacer referencia a esta y evitar ambigüedades en otros programas) para obtener las proyecciones del modelo ARIMA seleccionado por el algoritmo de automatización:

h.valid <- 12
 
# alternativa 1: direct
for1 <- forecast::forecast(auto.arima(AirPassengers), h=h.valid)
for1
         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1961       445.6349 430.8903 460.3795 423.0851 468.1847
Feb 1961       420.3950 403.0907 437.6993 393.9304 446.8596
Mar 1961       449.1983 429.7726 468.6240 419.4892 478.9074
Apr 1961       491.8399 471.0270 512.6529 460.0092 523.6707
May 1961       503.3945 481.5559 525.2330 469.9953 536.7937
Jun 1961       566.8624 544.2637 589.4612 532.3007 601.4242
Jul 1961       654.2602 631.0820 677.4383 618.8122 689.7081
Aug 1961       638.5975 614.9704 662.2246 602.4630 674.7320
Sep 1961       540.8837 516.9028 564.8647 504.2081 577.5594
Oct 1961       494.1266 469.8624 518.3909 457.0177 531.2356
Nov 1961       423.3327 398.8381 447.8273 385.8715 460.7939
Dec 1961       465.5076 440.8229 490.1923 427.7556 503.2596

Ahora utilizamos la función propia AutoArima_iterated para realizar el procedimiento descrito anteriormente. En este sentido se hace una proyección a un paso [forecast::forecast(for_x, h=1)], luego esta proyección se toma como observada [x<-c(x[-1],fhat)], aunque eliminamos la primera observación del vector x [x[-1]] esto para evitar la ventana se mueva de forma recursiva (creciendo una observación en cada paso), y se mueva entonces como una ventana móvil, manteniendo el tamaño de la ventana durante todo el horizonte. Una vez se modifica este valor observado, obtenemos un nuevo vector de datos observados para alimentar el modelo, pero teniendo el dato que proyectamos anteriormente como valor observado. Este proceso se repite n veces, como cuantos horizontes deseamos obtener.

# alternativa 2: iterated
source("AutoArima_iterated.R")
 
function(x,hn=h.valid) {
  for (hj in 1:hn){
    for_x <- auto.arima(x)
    fhat <- forecast::forecast(for_x, h=1)$mean
    x<-c(x[-1],fhat)  # -1 mantiene el tamano de la ventana
  }
  return(tail(x,h.valid))
}
 
for2 <-AutoArima_iterated(x=AirPassengers,hn=h.valid)
for2
 
       Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov
1961 445.6349 460.7414 470.6938 497.8221 502.7595 495.1462 484.4825 495.1744 496.0780 491.9802 487.1089
          Dec
1961 483.7122

Ahora, para representar ambas proyecciones en un gráfico usando la función autoplot, copiamos las propiedades guardada en la serie temporal del primer pronostico (for1) en este segundo vector de datos (for2). En este caso, recordemos que la función original de pronóstico es forecast y esta guarda el pronóstico en formato de serie temporal, mientras que el segundo lo hace como un vector.

 for2 <- ts(for2, star=start(for1$mean), frequency = frequency(for1$mean))
 
tail(AirPassengers,24) %>%
autoplot() +
  autolayer(for1$mean, color = "blue") +
  autolayer(for2, color = "red") +
  theme_minimal()

Finalmente, es relevante realizar una validación cruzada de estas proyecciones. Esto es principalmente útil en el marco de un sistema de pronósticos, pues nos ayuda a evaluar en el tiempo cual de las estrategias ha mostrado mayor precisión (menor error de pronósticos). En tal sentido usamos la función tsCV del paquete forecast de R dato que permite fácilmente obtener una serie histórica de los errores de pronóstico a cada horizonte.

1.       Primero definimos una función para proyectar valores (far1).

2.       Luego utilizamos la función de validación cruzada tsCV para obtener un histórico de errores a cada horizonte. Como son 6, aquí obtendremos una base de datos con 6 columnas y tantas observaciones como tenga nuestra data original.

3.       Finalmente obtenemos el RMSE de estos errores. Raíz del promedio de errores al cuadrado.

Sqrt(mean(e^2)) 

far1<- function(x, h){forecast::forecast(auto.arima(x), h=h)}
eAr1 <- tsCV(AirPassengers, forecastfunction = far1, h = h.valid, window=vent)
e_Ar1 <- sqrt(apply(eAr1^2, MARGIN=2, FUN=mean, na.rm = TRUE))

El procedimiento anterior lo repetimos, pero ahora usando la función AutoArima_iterated. Note que la única diferencia respecto a la evaluación anterior, es que usamos la función tsCv2, esta es la misma función tsCV, pero modificada para aceptar trabajar con cualquier pronóstico, dado que la función original solo trabaja con pronósticos arrojado en la estructura de los pronósticos del paquete forecast. Esto constituye una limitante para evaluar otras estrategias cuyo formato no coincidan. Ahora se modifica la función para que trabaje con cualquier pronostico siempre que coincida la longitud de los mismos. Dicha función modificada la colocamos como anexo en este documento.

# Auto.Arima paso a paso
source("tsCv2.R")
 
eAr1b <- tsCv2(AirPassengers, forecastfunction = AutoArima_iterated, h = h.valid, window=vent)
e_Ar1b <- sqrt(apply(eAr1b^2, MARGIN=2, FUN=mean, na.rm = TRUE))
e_Ar1b

 Ahora comparamos ambas estrategias, verificándose que la segunda estrategia es menos adecuada al intentar anticipar la serie, especialmente porque no parece capturar muy bien el componente estacional de la serie:

 t(rbind(e_Ar1,e_Ar1b))
    e_Ar1   e_Ar1b
h=1  13.42502 13.42502
h=2  16.01273 43.62504
h=3  18.80455 68.52383
h=4  20.42555 81.19967
h=5  20.51924 86.99055
h=6  20.98445 85.87618
h=7  21.54783 82.59922
h=8  22.37709 81.59745
h=9  22.52545 77.86818
h=10 23.80513 72.09909
h=11 24.62540 62.98470
h=12 25.83688 48.98253

Anexo:

tsCv2 <- function (y, forecastfunction, h = 1, window = NULL, xreg = NULL,
          initial = 0, ...)
{
  y <- as.ts(y)
  n <- length(y)
  e <- ts(matrix(NA_real_, nrow = n, ncol = h))
  if (initial >= n)
    stop("initial period too long")
  tsp(e) <- tsp(y)
  if (!is.null(xreg)) {
    xreg <- ts(as.matrix(xreg))
    if (NROW(xreg) != length(y))
      stop("xreg must be of the same size as y")
    xreg <- ts(rbind(xreg, matrix(NA, nrow = h, ncol = NCOL(xreg))),
               start = start(y), frequency = frequency(y))
  }
  if (is.null(window))
    indx <- seq(1 + initial, n - 1L)
  else indx <- seq(window + initial, n - 1L, by = 1L)
  for (i in indx) {
    y_subset <- subset(y, start = ifelse(is.null(window),
                                         1L, ifelse(i - window >= 0L, i - window + 1L, stop("small window"))),
                       end = i)
   
      fc <- forecastfunction(y_subset, h = h)
   
      e[i, ] <- y[i + seq(h)] - fc[seq(h)]
  }
  if (h == 1) {
    return(e[, 1L])
  }
  else {
    colnames(e) <- paste("h=", 1:h, sep = "")
    return(e)
  }
}