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.
library(readxl)
library(lubridate) # dates
library(tidyverse)
library(forecast)
library(timeSeries)
library(zoo)
# ********************************************************************
# ********************************************************************
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)
tipo_cambio <- cumsum(c(1, rnorm(n, mean = mu, sd = sigma))) # Inicia en 100 y suma los cambios diarios
fechas <- seq.Date(from = as.Date("2023-01-01"), by = "day", length.out = n + 1)
df <- data.frame(fecha = fechas, tipo_cambio = tipo_cambio)
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()
tipo_cambio_xts <- xts(tipo_cambio, order.by = fechas)
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.
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).
seq_for <- vent:(nrow(tipo_cambio_xts)-h.valid)
cont <- 1
data_Rbind <- list() # Lista que guarda resultados
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)]
# 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:
# 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
# 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")
# 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())
# 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
