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)
