Procedimientos recursivos y de ventana móvil para evaluar el rendimiento del pronóstico ARIMA (auto.arima en R)

En la siguiente entrada se muestra cómo utilizar bucles en R para rotar una ventana de estimaciones de distintas formas, con el fin de realizar un pronóstico y evaluar dicho pronóstico a un horizonte determinado utilizando el promedio de los errores de pronóstico obtenidos y la representación histórica de los errores de pronósticos. Las siguientes líneas generan una serie aleatoriamente (export) con la que ilustraremos el ejemplo, de la que posteriormente se crea una tasa de crecimiento discreta (tcexport) y se verifica su evolución gráfica.

library(forecast)
export <- cumsum(rnorm(100))
plot(export, type="l")
tcexport <- diff(export)/export[-length(export)]
plot(tcexport, type="l")
Ventana recursiva

Cargados los datos, la idea del procedimiento recursivo es sencilla, se parte de una muestra de tamaño determinado, elegimos una proporción de los datos totales para iniciar nuestras estimación, sobre esta proporción se estima un modelo ARIMA, usando la función auto.arima (aunque note que pudiese cambiar la metodología de estimación por alternativas como holt, ets u otras disponibles en el paquete forecast), se realiza el pronóstico a los próximo h pasos (6 en este ejemplo) y se comparan los valores que se proyectan con los valores observados durante dicho horizonte, de esta comparación se obtienen las medidas de rendimiento del pronóstico (se usa la función accuracy). Luego, se agrega una observación a la muestra usada para la estimación, es decir, si en la primera estimación la muestra iba de la observación 1 a la 50, en esta segunda ira de la 1 a la 51, y así sucesivamente se va agregando una observación en cada paso hasta n-h (n-h porque se debe contar siempre con h valores observados para poder comparar nuestros resultados).   

Para realizar el procedimiento anterior, primero se crean tres elementos: t es el total de observaciones; ventana es el intervalo de la ventana de entrenamiento, es decir, la que se usara para estimar el modelo con el cual se hará el pronóstico (en el ejemplo se usara el 70% de la muestra para realizar la primar estimación, posteriormente se agregará 1 observación en cada corrida); finalmente, n.valid representa el horizonte relevante.

t<-length(tcexport)
ventana<-round(t*0.7)
n.valid <- 6

El siguiente bucle utiliza un procedimiento recursivo donde obs es un contador que va desde ventana hasta t-n.valid (de 69 a 93 en este ejemplo), para luego indexado en el vector que se desea estimar tcexport [1:obs], ir agregando una observación (en el vector usado para estimar el modleo) en cada iteración del bucle. Posteriormente, a la serie elegida (tcexport [1:obs]) se le realiza un autoarima (del paquete forecast) y se realiza el pronóstico (xh). Finalmente se evalúa la precisión de este (xh) al comparar con los valores observados en ese periodo (xhobs). Nótese que esta comparación se realiza utilizando la función accuracy, también del paquete forecast. La matriz matrixRend la usaremos para guardar las series de rendimientos de los modelos, agregando los errores de medidas estimados como una fila adicional en cada estimación.

for (obs in ventana:(t-n.valid)){
  x<-tcexport[1:obs]
  model <- auto.arima(x)
  xh<-forecast::forecast(model, h=n.valid)$mean
 
  xhobs <- tcexport[(obs+1):(obs+n.valid)]
  matrixRend <- rbind(matrixRend, accuracy(xh,xhobs))
}

Ahora, la matrixRend muestra donde se han guardado los errores de pronósticos a 6 pasos, obtenido en cada una de las iteraciones. Una descripción del significado de cada una de las medidas de error, puede encontrarse fácilmente en la red.

head(matrixRend)

                   ME       RMSE        MAE MPE MAPE
Test set -0.007480647 0.11702889 0.10167476 100  100
Test set  0.035803108 0.10833884 0.09413361 100  100
Test set  0.012349578 0.09964723 0.08128739 100  100
Test set  0.036271693 0.11525404 0.10465435 100  100
Test set  0.019662793 0.09575623 0.08804545 100  100
Test set  0.024933024 0.09351676 0.08277522 100  100

Usando la función vectorizada apply podemos obtener el promedio global de cada medida de error:

apply(matrixRend, 2, mean)
          ME         RMSE          MAE          MPE         MAPE
  0.02475966   0.09155146   0.07858338 100.00000000 100.00000000

O la evolución histórica del error de pronóstico:

plot(matrixRend[,"RMSE"], type="l")
Ventana móvil

El procedimiento es el mismo del anterior, la única diferencia es que la ventana de estimación de los modelos mantiene fijo el punto de inicio, siempre agregándose una observación en cada paso. Es decir, ahora la ventana de estimación mantiene el mismo tamaño (en nuestro ejemplo 70 observaciones). Por ejemplo, en la primera estimación se usa una muestra que va desde 1:70, en la segunda 2:71, 3:72, y así sucesivamente hasta recorrer todo el rango.

Metodología
Paso 1
Paso 2
Paso 3
Recursiva
1:x
1:x+1
1:x+2
Móvil
1:x
2:x+1
3:x+2

Para realizar la estimación mediante ventana móvil solo se modifica la forma de recorrer el bucle. En tal sentido, i es ahora el contador que indica en cual observación inicia la muestra, que, al ser siempre del mismo tamaño, solo necesitamos sumar las siguientes 70 observaciones (ventana) observaciones posteriores [i:(i+ventana-1)].

matrixRend1 <- c()

for (i in 1:(t-ventana-n.valid+1)){
     x<-tcexport[i:(i+ventana-1)]
     model <- auto.arima(x)
     xh1<-forecast::forecast(model, h=n.valid)$mean
 
    xhobs <- tcexport[(obs+1):(obs+n.valid)]
     matrixRend1 <- rbind(matrixRend1, accuracy(xh1,xhobs))
  }

Una vez obtenido los resultados, podemos comparar el error promedio usando ambas estimaciones, notándose menor RMSE derivado de la segunda estimación:

rbind(recursivo = apply(matrixRend, 2, mean),
      movil = apply(matrixRend1, 2, mean))

                   ME       RMSE        MAE      MPE     MAPE
recursivo  0.024759664 0.09155146 0.07858338 100.0000 100.0000
movil     -0.005826735 0.08222809 0.06966551 124.4853 138.9471

También podemos comparar la evolución histórica de ambos error de pronóstico:

plot(matrixRend[,"RMSE"], type="b", pch=19, col = "blue")
lines(matrixRend1[,"RMSE"], col="red", type="b")
legend("topright", legend=c("Recursiva", "Móvil"),
       col=c("blue", "red"), lty=1:2, cex=0.8)

Entradas populares de este blog

Valores perdidos (NA) en R: identificación y tratamiento (I)

Ejemplos de tablas en Stata

Métodos y técnicas de Investigación Económica