La
siguiente entrada muestra una alternativa para obtener una serie histórica de
errores de predicción de una serie temporal (validación cruzada) a partir de
modelos ARIMA con represores (xARIMA). La entrada se segmenta en dos
partes, primero muestra un ejemplo aplicado del uso de la función y segundo explica
brevemente la construcción de la función de validación cruzada y de sus
funciones implícitas.
a.
Ejemplo de uso de la función
El ejemplo
utiliza datos del crecimiento trimestral del PIB real de la República
Dominicana, publicados por el Banco Central de la República Dominicana, para el
periodo 1991q1-2020q1.
library(forecast)
library(readxl)
library(tidyverse)
library(zoo)
setwd("C:/Users/Nery
Ramirez/Dropbox/densidadCrecimiento/densidadpib")
data_prib <- read_excel("Data_pib_real.xlsx")
> data_prib
# A tibble: 117 x 2
t pib_cons_real_2007
<chr>
<dbl>
1
1991Q1 20.5
2
1991Q2 23.9
3
1991Q3 22.5
4
1991Q4 24.4
5
1992Q1 26.7
6
1992Q2 26.6
7
1992Q3 29.1
8
1992Q4 26.8
9
1993Q1 30.0
10 1993Q2 32.6
# ... with 107 more rows
Utilizando
el paquete dplyr del tidyverse y la función as.yearqtr de la librería zoo, se recuperan las fechas correspondientes a cada
trimestre, se calcula el crecimiento anual (utilizando cuatro rezagos en la función
lag) y posteriormente se crea una variable binaria
llamada d_cre_b1, la misma asume el valor de 1 en
los periodos de crecimiento negativo posterior al 2002, y en los 10 trimestres
donde el crecimiento fue menor antes de este año (esto porque no se registraron
periodos de caídas interanual en la producción trimestral).
data_prib1 <- data_prib %>%
mutate(q =
as.yearqtr(time(pib_cons_real_2007), format = "q"),
anio = as.numeric(substring(t, 1,4)),
creci = (pib_cons_real_2007/dplyr::lag(pib_cons_real_2007,4)) - 1,
d_cre_b1 = (creci %in% sort(creci[anio<=2002])[1:10] | creci<0)*1,
) %>%
#
Reordenar el nombre de las variables
select(t,q,everything())
data_prib1
# A tibble: 117 x 6
t q
pib_cons_real_2007 anio creci d_cre_b1
<chr> <yearqtr> <dbl> <dbl> <dbl>
<dbl>
1 1991Q1 1
Q1 20.5 1991 NA NA
2 1991Q2 2
Q1 23.9 1991 NA NA
3 1991Q3 3
Q1 22.5
1991 NA NA
4 1991Q4 4
Q1 24.4 1991 NA NA
5 1992Q1 5
Q1 26.7 1992
0.299 0
6 1992Q2 6
Q1 26.6 1992
0.114 0
7 1992Q3 7
Q1 29.1 1992
0.294 0
8 1992Q4 8
Q1 26.8 1992
0.101 0
9 1993Q1 9
Q1 30.0 1993
0.124 0
10 1993Q2 10 Q1 32.6 1993
0.226 0
# ... with 107 more rows
Ahora,
usamos la función (se explica debajo) error_historico_d1 para obtener una serie histórica de
los errores de pronósticos a diversos pasos. La serie x es la serie de crecimiento anteriormente estimada,
pero omitiendo los na, mientras que el regresor x
del modelo
corresponde a la dummy para indicar periodos de crecimiento
económico “bajo” o negativos. Los argumentos de la función son (note que solo
los dos primeros argumentos son obligatorios):
-
x: serie de tiempo sobre la que se
desea trabajar.
-
xreg_d: regresores del modelo ARIMA.
-
n.valid: número de periodos que se desean
proyectar.
-
Recursiva: indicar que queremos que la
ventana de estimación sea móvil o recursiva.
-
Ventana: tamaño de la ventana, el valor
predeterminado es 70% de las observaciones.
-
Medida_error: indicamos la serie de errores que
deseamos guardar. Las opciones disponibles se colocan debajo.
("ME",
"MSE", "RMSE", "MAE", "MPE",
"MAPE","ACF1", "Theil's U")
-
Agregar: se muestra un agregado de las
series de errores, puede ser un promedio, mediana o la función de interés,
siempre que esté disponible en R.
x<-data_prib1$creci[-c(1:4)]
xreg_d<-data_prib1$d_cre_b1[-c(1:4)]
error_historico_d1(x,
xreg_d,
n.valid=8,
recursiva=TRUE,
ventana=round(length(!is.na(x))*0.7),
medida_error="RMSE",
agregar = "mean")
El
modelo arroja por salida una lista con información sobre la estimación, el
ultimo modelo autoarima seleccionado por el algoritmo desarrollado en el
paquete forecast, la serie histórica de errores a
cada paso, y los errores agregados, que para el caso actual representa el
promedio de errores.
$Informacion
[,1]
[1,] "Medida de errro: RMSE"
[2,] "Horizonte: 8"
[3,] "Función agregar: mean"
$`Modelo Arima`
Series: x.train
Regression with ARIMA(1,0,0) errors
Coefficients:
ar1 intercept xreg
0.3090 0.1407 -0.1949
s.e.
0.1059 0.0120 0.0209
sigma^2 estimated as 0.005801: log likelihood=122.84
AIC=-237.68
AICc=-237.28 BIC=-227.06
$`Hitorico errores`
h1 h2 h3 h4 h5 h6 h7 h8
error_pronH 0.0120526909 0.01155399 0.01390647
0.01437705 0.01369995 0.03722491 0.04796603 0.04526843
error_pronH 0.0158480814 0.01798974 0.01688961
0.01571486 0.04056176 0.05165188 0.04824761 0.04778444
error_pronH 0.0138126750 0.01552025 0.01395166
0.04467589 0.05622090 0.05179590 0.05078809 0.04789154
error_pronH 0.0224250648 0.01677765 0.05223426
0.06335995 0.05719506 0.05523534 0.05154089 0.04958089
error_pronH 0.0164286782 0.06017829 0.07053814
0.06166394 0.05861779 0.05396812 0.05146070 0.05575162
error_pronH 0.0898086733 0.08977179 0.07400109
0.06776404 0.06108429 0.05734009 0.06103166 0.05711851
error_pronH 0.0597617913 0.04254667 0.04443076
0.03958809 0.03817338 0.04751249 0.04402317 0.04143310
error_pronH 0.0164400564 0.04262606 0.03739458
0.03541055 0.04663292 0.04258040 0.03962853 0.04393133
error_pronH 0.0517267059 0.03950836 0.03657621
0.04995537 0.04469983 0.04105940 0.04551537 0.05180642
error_pronH 0.0009055807 0.02340388 0.05000653
0.04339390 0.03929627 0.04430251 0.05116947 0.04800552
error_pronH 0.0334304486 0.06134581 0.05018948
0.04400640 0.04857884 0.05530050 0.05134836 0.05356872
error_pronH 0.0671412866 0.04759847 0.03959037
0.04701842 0.05568364 0.05103217 0.05349869 0.06633577
error_pronH 0.0214602097 0.01728465 0.04048314
0.05432041 0.04890278 0.05193941 0.06728899 0.06510779
error_pronH 0.0199991749 0.04827238 0.06183984
0.05388911 0.05639750 0.07222707 0.06924120 0.06718595
error_pronH 0.0729692207 0.07808759 0.06415510
0.06439812 0.08009231 0.07560880 0.07257203 0.07483408
error_pronH 0.0552389712 0.03981395 0.05033203
0.07538330 0.07084374 0.06789383 0.07156330 0.06906939
error_pronH 0.0101622099 0.04823030 0.08075583
0.07418822 0.07001982 0.07401316 0.07095847 0.07125072
error_pronH 0.0635089143 0.09744160 0.08450430
0.07734165 0.08030907 0.07596078 0.07555640 0.07396403
error_pronH 0.1471701781 0.10949240 0.09397442
0.09300588 0.08597771 0.08379617 0.08113644 0.07623440
error_pronH 0.1023847115 0.08070527 0.08495009
0.07743285 0.07689872 0.07435431 0.06939065 0.06493957
error_pronH 0.0833769891 0.08756248 0.07637803
0.07564973 0.07301224 0.06724270 0.06227617 0.06112950
error_pronH 0.1048222591 0.08473302 0.08124879
0.07688083 0.06939527 0.06336327 0.06178965 0.05841032
error_pronH 0.0266498858 0.04711935 0.05482717
0.04827173 0.04317618 0.04416090 0.04167408 0.03904759
error_pronH 0.0530020671 0.06221638 0.05169243
0.04476714 0.04557993 0.04250907 0.03942900 0.04917571
error_pronH 0.0887500842 0.06311191 0.05156113
0.05037165 0.04590345 0.04195031 0.05175833 0.05205893
error_pronH 0.0388291674 0.02796296 0.03736565
0.03424014 0.03078301 0.04705677 0.04843461 0.04533792
error_pronH 0.0049771555 0.03339035 0.02975767
0.02594643 0.04714780 0.04861765 0.04505684 0.04292588
$`Errores agregados`
h1 h2 h3 h4 h5 h6 h7 h8
0.04789196 0.05163872 0.05346425 0.05381539
0.05499571 0.05628511 0.05645721 0.05626474
a.
Funciones
La función
error_historico_d1 utiliza la función EvaluacionP, que es en realidad una modificación
de la función acurracy del paquete forcast. La misma obtiene una serie de
medidas de error (me, mse, sqrt(mse), mae, mpe, mape, r1, theil) de pronóstico para el periodo de validación
de las series (test data) (ver referencia en https://otexts.com/fpp2/accuracy.html), al comparar el valor observado
(x) con el proyectado (f).
# Solo se reescribe la original para hacerla menos
restrictiva
EvaluacionP <- function(f, x){
# versión
6/12/2019
n <-
length(x)
test <-
1:n
error
<- (x - f)[test]
pe <-
error / x[test] * 100
me <-
mean(error, na.rm = TRUE)
mse <-
mean(error ^ 2, na.rm = TRUE)
mae <-
mean(abs(error), na.rm = TRUE)
mape <-
mean(abs(pe), na.rm = TRUE)
mpe <-
mean(pe, na.rm = TRUE)
fpe <-
(c(f[2:n]) / c(x[1:(n - 1)]) - 1)[test - 1]
ape <-
(c(x[2:n]) / c(x[1:(n - 1)]) - 1)[test - 1]
theil
<- sqrt(sum((fpe - ape) ^ 2, na.rm = TRUE) / sum(ape ^ 2, na.rm = TRUE))
# Se
coloca este Condicional para poder evaluar pronosticos a 1 paso
if
(length(f)==1){
r1 <-
NA
} else {
r1 <-
acf(error, plot = FALSE, lag.max = 2, na.action = na.pass)$acf[2, 1, 1]
}
out <-
c(me, mse, sqrt(mse), mae, mpe, mape, r1, theil)
names(out)
<- c("ME", "MSE", "RMSE", "MAE",
"MPE", "MAPE","ACF1", "Theil's U")
return(out)
}
Para
correr la función para cada uno de los periodos seleccionados, se utiliza un
pequeño bucle, que modifica los vectores del set de entrenamiento (vistos en la
función anterior), para a partir de su longitud evaluar los errores a distintos
pasos. Es decir, especificando en el argumento horizonte, los pasos que
deseamos evaluar, podemos obtener la evaluación del error a dichos pasos. De
manera predeterminada están los periodos 1-6, y el 12.
EvaluacionPh <- function(f, x,
horizontes=c(1:6,12)){
#
vh <-
horizontes
Resultado<-c()
for(i in
vh){
performance<-as.vector(EvaluacionP(f[1:i],
x[1:i]))
Resultado<-rbind(Resultado, performance)
}
#
Resultados de la función
row.names(Resultado) <- as.character(horizontes)
colnames(Resultado) <- c("ME", "MSE",
"RMSE", "MAE", "MPE",
"MAPE","ACF1", "Theil's U")
return(Resultado)
}
Finalmente,
llegamos a la función de validación cruzada que usamos anteriormente a partir
de una ventada recursiva, es decir, que la ventana de estimación se va expandiendo
en la medida en que nos desplazamos en el tiempo. La función trabaja de forma
muy sencilla:
a. Segmenta la data en una sección de estimación
y otra de validación (test data). Note esta se crean en función de los pasos
que queremos evaluar. En el caso de una serie con 100 observaciones en la que
deseamos proyectar tres pasos, las primeras 97 serán del set de estimación y
las últimas tres del set de validación. Note que esto se logra mediante indexación
numérica, por tanto, data estimación: observaciones ength(data.ejercicio)
- n.valid; data
entrenamiento: test [-c(1:n.train)].
b. Luego se utiliza la data de estimación
o entrenamiento para estimar el modelo ARIMA, incluyendo la dummy, así como la proyección
correspondiente.
auto.arima(x.train,
xreg=d.train, seasonal=T)
c. Posteriormente se usa la función EvaluacionPh, para evaluar el pronóstico
obtenido en la secuencia de pasos seleccionados (n.valid)
y recuperar la medida de error indicada, mediante indexación del vector obtenido
([,medida_error]) en la función de evaluación,
llamado resultados.
EvaluacionPh(f=x.hat, x=x.test,
h=1:n.valid)[,medida_error]
d. Note que la selección del set de
entrenamiento puede ser móvil o recursiva. La clave de esto es el contador
incluido en el bucle (c<-1), dado que las observaciones del set de entrenamiento
es x[c:obs] y la evolución de c esta
condicionada a al argumento recursiva de la función if (recursiva == FALSE){c <- c+1}, en el caso verdadero no se cumple
la condición y c vale 1 en todas las vueltas que de el bucle, así la ventana de
estimación va desde 1 hasta obs (que note va creciendo en el bucle que crea los
set de entrenamiento (for (obs in
(ventana+n.valid):t_obs)).
En el caso sea falso, c iria sumando 1 en cada vuelta, dejando una ventana del
mismo tamaño en cada vuelta, en vez de ir creciendo.
error_historico_d1 <- function(x,
xreg_d,
ventana=round(length(!is.na(x))*0.7),
n.valid,
recursiva=TRUE,
medida_error="RMSE",
agregar =
"mean"){
#
********************************************************************
# x = time serie vector
# ventana
= longitud del vector de estimación. Default= 70% de la muestra
# n.valid
= periodos proyectados
#
medida_error:
ME;MSE;RMSE;MAE;MPE;MAPE;ACF1;Theil's U
# agregar:
qué función deseas agregar de los errores.
# Ejemplo:
error_historico(AirPassengers, ,ventana=70, n.valid=5, recursiva=TRUE, medida_error="RMSE" # ,agregar
= "median")
t_obs
<- length(x[!is.na(x)])
Errores_histo
<- c()
c<-1
for (obs
in (ventana+n.valid):t_obs){
# Extrae
las observaciones con que se van a trabajar
data.ejercicio <- x[c:obs]
dummy.ejercicio <- xreg_d[c:obs]
#
Prepara los datos:
n.train
<- length(data.ejercicio) - n.valid #
intervalo de estimación
x.train
<- data.ejercicio[1:n.train] #
Datos del modelo
x.test
<- data.ejercicio[-c(1:n.train)]
d.train
<- dummy.ejercicio[1:n.train]
d.test
<- dummy.ejercicio[-c(1:n.train)]
# Modelo
Autoarima Móvil
modelo
<- auto.arima(x.train, xreg=d.train, seasonal=T)
x.hat
<- forecast::forecast(modelo, xreg=d.test, h=n.valid)$mean
error_pronH
<- EvaluacionPh(f=x.hat, x=x.test, h=1:n.valid)[,medida_error]
Errores_histo <- rbind(Errores_histo,error_pronH)
if
(recursiva == FALSE){c <- c+1}
#"Movil"
} #Cierrra
el bucle
colnames(Errores_histo) <- paste("h",1:n.valid,
sep="")
errores_h
<- apply(Errores_histo, MARGIN=2, FUN=get(agregar), na.rm = TRUE)
info <-
rbind(paste("Medida de errro: ", medida_error, sep=""),
paste("Horizonte: ", n.valid,
sep=""),
paste("Función agregar: ", agregar, sep="") )
resultados<-list("Informacion"=info,
"Modelo Arima"=modelo,
"Hitorico errores"=Errores_histo,
"Errores
agregados"=errores_h)
return(resultados)
} #Cierra la función
