En el
siguiente ejemplo se muestra una validación cruzada de series temporales en R. Necesita
evaluar la capacidad de pronóstico de un modelo para elegir alguna alternativa
en función de su rendimiento. Es decir, necesito verificar cual hubiese sido el
pronóstico histórico de una serie temporal, por ejemplo con un modelo AR(1),
para comparar con la precisión de otros modelos y decidir cual estrategia seguir
en lo adelante. En este sentido, necesitaríamos realizar el ejercicio de proyección
en distintos puntos temporales, para posteriormente recuperar los errores y
comparar estas medidas de error entre distintos modelos.
1.
Simulamos
la serie temporal
(en el caso de una base de datos, deberías importarlas en R).
library(writexl)
# Exportar a Excel
library(readxl)
library(lubridate)
# dates
library(tidyverse)
library(forecast)
library(timeSeries)
library(zoo)
# Importando datos ---------------
#
********************************************************************
#
********************************************************************
# Configuración inicial
set.seed(4)
# Asegura reproducibilidad
n
<- 250 # Número de días (simulación de aproximadamente un año laboral)
mu
<- 0 # Promedio del cambio diario (sin
tendencia)
sigma
<- 0.01 # Volatilidad diaria
# Generación de la caminata aleatoria
tipo_cambio
<- cumsum(c(1, rnorm(n, mean = mu, sd = sigma))) # Inicia en 100 y suma los cambios diarios
# Creación de un data.frame con fechas
fechas <-
seq.Date(from = as.Date("2023-01-01"), by = "day",
length.out = n + 1)
df
<- data.frame(fecha = fechas, tipo_cambio = tipo_cambio)
# Gráfico de series temporales
library(ggplot2)
ggplot(df,
aes(x = fecha, y = tipo_cambio)) +
geom_line(color = "blue") +
labs(title = "Simulación de Caminata
Aleatoria del Tipo de Cambio",
x = "Fecha",
y = "Tipo de Cambio") +
theme_minimal()
library(xts)
tipo_cambio_xts
<- xts(tipo_cambio, order.by = fechas)
2.
Validación
cruzada
En
resumen, lo que se hace es definir una ventana de estimación y otra de validación
(pronóstico). En este ejemplo son de 60 y 5 observaciones respectivamente. 60 observaciones
serán usadas para estimar el modelo, mientras que 5 será la longitud de la
ventana del pronóstico. De esta forma, en la primera iteración se va estimar los
modelos con observaciones del 1:60 y otra ventana de pronóstico de la
61-65, luego 2:61
y 62:66, … y así respectivamente hasta recorrer
toda la muestra. El documento de trabajo “Pronóstico del consumo privado en la
República Dominicana” de Alejandro J. Balcácer y Nerys F. Ramírez, muestra en detalle
esta estrategia.
# Análisis de precisión. Error histórico ---------------h.valid
<- 5
horizonte
<- paste0("h", 1:h.valid)
vent
<- 60
Posteriormente,
se genera un for para realizar este recorrido usando
una serie. Primero se define el valor del ultimo dato observado en cada iteración
del bucle (seq_for), este dato sirve de referencia
para el resto de las ventanas. Además, se define un contador para el punto de
inicio de la ventana de estimación (cont). Los puntos clave de esta
secuencia son:
-
Indicar
la secuencia con la posición de las observaciones que vamos a usar en cada punto
de la serie temporal [obser
in seq_for].
-
La
serie temporal se segmenta en dos partes: la estimación (tipo_cambio_xts[cont:obser]) y la otra parte de pronósticos (tipo_cambio_xts[(obser+1):(obser+h.valid)]). Fíjese que la primera va desde
el contador hasta obser, recuerde que obser es el valor final hasta donde irá
la ventana de estimación y el valor que en el presente ejemplo sirve de
referencia. Luego, la parte para evaluar las proyecciones se definen a partir
de la última observación de la parte de estimación (obser+1):(obser+h.valid).
-
La
tercera parte importante del bucle es incluir el modelo de pronóstico. Aquí debemos
usar alguna estrategia para proyectar los periodos que nos resulten de interés.
En este ejemplo son 5 observaciones, y se usa un modelo ar(1).
-
Finalmente,
se crea un data.frame para guardar el pronóstico, el
dato observado, los horizontes y el nombre del model.
Además, se guarda cada data.frame en una lista creada anteriormente
(data_Rbind). Primero se llaman todos los objetos
que inician con “df_” (tener pendiente guardar cada
nuevo modelo con ese nombre).
# Esto lo que hace, es que permite recorrer las ventanas
seq_for
<- vent:(nrow(tipo_cambio_xts)-h.valid)
# OJO CORRER EL BULCE DESDE AQUI
cont
<- 1
data_Rbind
<- list() # Lista que guarda resultados
for
(obser in seq_for){
tc_train <- tipo_cambio_xts[(obser-vent+1):obser]
tc_test <-
tipo_cambio_xts[(obser+1):(obser+h.valid)]
fecha_ciclo <- fechas[obser]
fecha_hat <- fechas[(obser+1):(obser+h.valid)]
ordena_data
<- function(objetos,...){
# objetos = lista de nombre de vectore
# Iterar sobre los nombres de los objetos y
cambiar los nombres de las columnas
for (i in seq_along(objetos_df)) {
xdf_name <- get(objetos_df[i])
new_column_names <-
c("fechaCiclo","fechaHat","horizonte","obs","pro","model")
names(xdf_name) <- new_column_names
assign(objetos_df[i], xdf_name)
}
lista_objetos <- lapply(objetos_df, get)
df1_final <-
dplyr::bind_rows(lista_objetos)
return(df1_final)
}
### 1. ESTRATEGIA
AR(1) . Ventana móvil ************************
model_ar1 <- Arima(tipo_cambio_xts,
order=c(1,0,0))
ht_f_ar1 <-
forecast::forecast(model_ar1, h=h.valid)$mean
df_ar1 <- data.frame(fecha_ciclo,
fecha_hat, horizonte, #fechas
obsT=tc_test, # dato
observado
proT=ht_f_ar1,
model="1.Ar(1)_móvil")
### 1. ESTRATEGIA
AR(2,1) . Ventana móvil ************************
model_ar21 <- Arima(tipo_cambio_xts,
order=c(2,0,1))
ht_f_ar21 <- forecast::forecast(model_ar21,
h=h.valid)$mean
df_ar21 <- data.frame(fecha_ciclo,
fecha_hat, horizonte, #fechas
obsT=tc_test, # dato
observado
proT=ht_f_ar21,
model="2.Ar(21)_móvil")
###
******************************************************************
###
******************************************************************
objetos
<- ls()
objetos_df
<- objetos[grep("^df_", objetos)]
data_Rbind[[(obser-vent+1)]]
<- ordena_data(objetos_df)
#
Esto al final siempre
cont
<- cont+1 # cont y esto es equivalente
(obser-vent+1):
}
Luego
de culminar, se utiliza una función propia para ordenar la data (ordena_data) y convertirla en un único data.frame con toda la data:
data_forecast_final
#
A tibble: 1,870 × 6
fechaCiclo fechaHat horizonte
obs proy model
<date> <date> <chr> <dbl> <dbl> <chr>
1 2023-03-01 2023-03-02 h1 0.969 1.01 1.Ar(1)_móvil
2 2023-03-01 2023-03-03 h2 0.981 1.01 1.Ar(1)_móvil
3 2023-03-01 2023-03-04 h3 0.972 1.01 1.Ar(1)_móvil
4 2023-03-01 2023-03-05 h4 0.978 1.00 1.Ar(1)_móvil
5 2023-03-01 2023-03-06 h5 0.987 1.00 1.Ar(1)_móvil
6 2023-03-01 2023-03-02 h1 0.969 1.01 1.Ar(21)_móvil
7 2023-03-01 2023-03-03 h2 0.981 1.01 1.Ar(21)_móvil
8 2023-03-01 2023-03-04 h3 0.972 1.01 1.Ar(21)_móvil
9 2023-03-01 2023-03-05 h4 0.978 1.00 1.Ar(21)_móvil
10
2023-03-01 2023-03-06 h5 0.987
0.999 1.Ar(21)_móvil
Con esta
data, ahora se puede realizar cualquier análisis sobre los errores promedio, su
evolución histórica o el comportamiento de estos errores condicionados a eventos
económicos, por ejemplo, los errores durante el COVID o periodos de alta
incertidumbre.
# Cuadrando la base de los errores
# Crea la data final con todos los los dtos necesarios
data_forecast_final
<- dplyr::bind_rows(data_Rbind) |> as_tibble()
names(data_forecast_final)
<-
c("fechaCiclo","fechaHat","horizonte","obs","proy","model")
# Errores total
# Total
data_forecast_final
|>
#dplyr:: filter(fechaCiclo >
as.Date("2018-01-01")) |>
mutate(error=obs-proy) |>
group_by(horizonte, model) |> #
(horizonte, model, COVID)
summarize(
mae_e_abs = mean(abs(error), na.rm = TRUE),
) |>
pivot_wider(names_from = horizonte,
values_from = mae_e_abs) |>
arrange(h1) |>
dplyr::select(model,horizonte, everything())
`summarise()` has grouped output by 'horizonte'. You can override
using the `.groups` argument.
# A tibble: 2 × 6
model h1 h2
h3 h4 h5
<chr> <dbl> <dbl>
<dbl> <dbl> <dbl>
1
1.Ar(21)_móvil 0.0479 0.0442 0.0405 0.0377 0.0350
2
1.Ar(1)_móvil 0.0483 0.0443 0.0409 0.0378
0.0353
