Evaluación de pronóstico en R: validación cruzada de modelos de series temporales

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

Descomposición de una serie temporal en R: componente transitorio vs permanente

La siguiente entrada utiliza la inflación mensual de la República Dominicana para obtener una descomposición histórica de la evolución de la inflación trimestral en un componente transitorio vs. componente subyacente. Primero se activan las librearas requeridas y se importa la data (dataset), que incluye la fecha mensual y el IPC con datos usando como base la canasta 2019-2020.

# Load the data
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
 
library(readxl)
library(forecast)
library(lubridate)
library(tidyverse)
 
dataset <- read_excel("data.xlsx")
 
# A tibble: 472 x 2
   fecha                 ipc
   <dttm>              <dbl>
 1 1984-01-01 00:00:00  1.38
 2 1984-02-01 00:00:00  1.42
 3 1984-03-01 00:00:00  1.44
 4 1984-04-01 00:00:00  1.46
 5 1984-05-01 00:00:00  1.48
 6 1984-06-01 00:00:00  1.54
 7 1984-07-01 00:00:00  1.56
 8 1984-08-01 00:00:00  1.57
 9 1984-09-01 00:00:00  1.64
10 1984-10-01 00:00:00  1.68
# ... with 462 more rows

Posteriormente, se agrega la serie trimestralmente usando el promedio del IPC. Se identifica el trimestre de cada mes (quarter), posteriormente se crea un agregado del promedio trimestral del IPC, y sobre este IPC trimestral, calculamos la inflación trimestral (q_inf). En una entrada anterior se explicó en detalle la agregación temporal y las transformaciones de series temporales. 

 q_dataset <- dataset |>
  mutate(quarter = zoo::as.yearqtr(fecha, "%Q")) |> # create a new variable for the quarter
  group_by(quarter)  |>  # group by the quarter variable
  summarise(q_ipc = mean(ipc)) |>
  mutate(q_inf = ((q_ipc/dplyr::lag(q_ipc))-1)*100) |>
  slice(-1)
 
# A tibble: 157 x 3
   quarter   q_ipc q_inf
   <yearqtr> <dbl> <dbl>
 1 1984 Q2    1.49  5.66
 2 1984 Q3    1.59  6.53
 3 1984 Q4    1.76 11.0
 4 1985 Q1    2.11 19.6
 5 1985 Q2    2.24  6.43
 6 1985 Q3    2.32  3.53
 7 1985 Q4    2.41  3.67
 8 1986 Q1    2.44  1.23
 9 1986 Q2    2.39 -2.01
10 1986 Q3    2.42  1.40
# ... with 147 more rows

Ahora, definimos una inflación trimestral como una serie temporal usando la función ts.

 ts_data <- ts(q_dataset$q_inf, start = c(1984,2), frequency = 4)
 
            Qtr1        Qtr2        Qtr3        Qtr4
1984              5.66392672  6.52527285 10.99767476
1985 19.64583811  6.43449238  3.53425010  3.66683217
1986  1.23093134 -2.00846916  1.39780793  4.35942001
1987  1.45783902  4.75590834  5.30846177  6.77347683

Ahora, usamos la función stl (ver referencia: https://otexts.com/fpp2/stl.html) para obtener una descomposición de la serie, recuperando los componentes estacionales, la tendencia y el componente aleatorio de las series.

# Decompose the time series using the STL function
decomp <- stl(ts_data, s.window = "periodic")
 
# Extract the seasonal, trend, and remainder components
seasonal <- decomp$time.series[, "seasonal"]
trend <- decomp$time.series[, "trend"]
remainder <- decomp$time.series[, "remainder"]

Ahora, agregamos los componentes (remainder) transitorios vs. permanentes (trend + seasonal).

q_dataset$permanent <- trend + seasonal
q_dataset$transient <- remainder
 
# A tibble: 157 x 5
   quarter   q_ipc q_inf permanent transient
   <yearqtr> <dbl> <dbl>     <dbl>     <dbl>
 1 1984 Q2    1.49  5.66    5.67    -0.00924
 2 1984 Q3    1.59  6.53    8.72    -2.19  
 3 1984 Q4    1.76 11.0    10.7      0.325 
 4 1985 Q1    2.11 19.6    11.4      8.23  
 5 1985 Q2    2.24  6.43    8.27    -1.83  
 6 1985 Q3    2.32  3.53    5.81    -2.28  
 7 1985 Q4    2.41  3.67    3.29     0.376 
 8 1986 Q1    2.44  1.23    1.57    -0.340 
 9 1986 Q2    2.39 -2.01   -0.0908  -1.92  
10 1986 Q3    2.42  1.40    1.39     0.00589
# ... with 147 more rows

Finalmente, se crea un gráfico apilado asumiendo la incidencia de cada componente.  

q_dataset |>
  dplyr::filter(quarter >= "2000 Q1") |>
  gather(id, value, -c(quarter,q_ipc,q_inf)) |>
  ggplot(aes(x = quarter, y = value, fill = id)) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c("#56B4E9", "#E69F00"), name = element_blank(), labels = c("Permantente", "Transitorio")) +
  geom_line(aes(x = quarter, y = q_inf), size=0.8) +
  theme_classic() +
  theme(legend.position = "bottom")

Sin embargo, vemos que la metodología anterior asume los componentes estacionales como transitorios. Por lo que, la política monetaria estaría reaccionando a elementos estacionales al considerarlos como elementos permanentes, sin embargo, entendemos que este elemento debe agregarse como elemento transitorio.

Gráfico de recesión (Reccesion plot) en R y ggplot2

En la siguiente entrada se muestra un ejemplo de gráficos con bandas de recesión en ggplot2 de R. La idea se toma de “rstudio-pubs”[1]. Desde Excel, se cargan los datos conteniendo los precios mensuales del precio del petróleo.

library(readxl)
library(dplyr)
library(ggplot2)
library(tidyr) #gather
library(lubridate)
 
series_mes <- read_excel("series_examples.xlsx", sheet = "mes")
 
series_mes  %>%
  select(fecha,p_wti) %>%
  na.omit()
 
# A tibble: 446 x 2
   fecha               p_wti
   <dttm>              <dbl>
 1 1986-01-01 00:00:00  22.9
 2 1986-02-01 00:00:00  15.5
 3 1986-03-01 00:00:00  12.6
 4 1986-04-01 00:00:00  12.8
 5 1986-05-01 00:00:00  15.4
 6 1986-06-01 00:00:00  13.4
 7 1986-07-01 00:00:00  11.6
 8 1986-08-01 00:00:00  15.1
 9 1986-09-01 00:00:00  14.9
10 1986-10-01 00:00:00  14.9
# ... with 436 more rows

Posteriormente, se crea una base de referencia llamada reccess, donde se coloca una variable con las fechas de inicio y final de cada crisis, estas son llamadas begin y end, respectivamente. Respecto al código original, aquí se agrega la crisis financiera de 2003 dada su relevancia para la economía dominicana, además de la crisis del COVID-19. Adicionalmente, en la variable evento coloca el nombre del evento asociado a cada crisis. Esto le permitirá colocar etiquetas sobre el gráfico posteriormente.

recess <- data.frame(
  begin = c("1969-12-01","1973-11-01","1980-01-01","1981-07-01",
             "1990-07-01","2001-03-01","2003-01-01","2007-12-01",
             "2020-03-01"),
   end = c("1970-11-01","1975-03-01","1980-07-01","1982-11-01",
           "1991-03-01","2001-11-01","2004-01-01","2009-07-30",
           "2021-12-01"),
 event = c( "Fiscal & Monetary ", "1973 Oil crisis","Double dip I",
            "Double dip II", "Oil price shock", "Dot-com bubble",
             "Crisis Bancaria RD","Sub-prime crisis", "COVID-19"),
  y =  c(.014, 0.020, 0.029,  0.0341,  0.027, 0.02,0.02, 0.025, 0.2),
  stringsAsFactors = F
)

La variable y es la referencia sobre el valor de ese eje en donde entendemos debe aparecer las etiquetas asociadas a cada crisis. Es decir, que debemos cambiar estos valores para cada serie que quisiéramos representar. Sin embargo, dado que tenemos un vector de fechas en nuestra serie temporal este proceso es necesario automatizar. La opción utilizada en este ejemplo utiliza la función %in% combinada con which para identificar las posiciones de fechas en nuestra base de datos que coincide con las crisis. Estas posiciones se guardan en el vector recesion_index. Finalmente, se usa la indexación numérica para recuperar los calores de mi serie que coinciden con estas fechas. Por ejemplo, imaginemos que una recesión inició en enero de 2003, la idea es buscar el valor de mi serie en ese mes.  Estos valores se guardan en el vector position_labels.

#ubicando las posiciones de las fechas
recesion_index <- which(format(series_mes$fecha, "%y-%b") %in% format(recess$begin, "%y-%b"))
recesion_index
[1]  78 206 228 287 434
 
position_labels <- series_mes$p_wti[recesion_index]
position_labels
[1] 16.70 29.61 29.46 94.77 50.54

Finalmente, se usan diferentes geom de ggplot2 para obtener la representación deseada.

·         geom_line, gráfico de línea con el precio del petróleo.

·         geom_rect, franja roja asociada a cada crisis.

·         geom_label, posición de las etiquetas.

ggplot(series_mes, aes(x = fecha, y = p_wti)) +
  geom_line(size=1.1) +
  geom_rect(data = dplyr::filter(recess, begin >= as.Date("1986-01-01")),
            aes(xmin = begin, xmax = end, ymin = -Inf, ymax = +Inf, fill = "Recession"),
            inherit.aes = FALSE, alpha = 0.2) +
  geom_label(data = dplyr::filter(recess, begin >= as.Date("1986-01-01")),
             aes(x = end, y = position_labels, label=event), size = 3) +
  scale_fill_manual(name = "", values="red", label="Recessions")  +
  theme_minimal() +
  ggtitle(c("Precio del petróleo (WTI) \n Datos mensuales 1986-2023")) +
  theme(legend.position = "none") +
  ylab(NULL) +  xlab(NULL)
 
ggsave("plot0.png")

---

[1] Introduction to Visualization with ggplot2.