Pronostico de Ventas para diferentes tiendas

Store Item Demand Forecasting Challenge

0.1 Introducción

La intención de este ejercicio es poder estudiar las librerias de trabajo existentes para pronósticar series de tiempo; De forma sencilla y detallada, se busca desarrollar un código para el proposito ya mencionado, por lo que se plantea a continuación la estrategia desarrollada para el ejercicio.

0.1.1 Descripción del reto

Esta competencia se ofrece como una forma de explorar diferentes técnicas de series de tiempo en un conjunto de datos relativamente simple y limpio.

Se le proporcionan 5 años de datos de ventas de artículos de la tienda y se le pide que prediga 3 meses de ventas para 50 artículos diferentes en 10 tiendas diferentes.

¿Cuál es la mejor manera de lidiar con la estacionalidad? ¿Las tiendas se deben modelar por separado o se pueden agrupar? ¿El aprendizaje profundo funciona mejor que ARIMA? ¿Puede alguno superar a xgboost?

Esta es una gran competencia para explorar diferentes modelos y mejorar tus habilidades de pronóstico.

0.1.2 Descripción del conjunto de datos

El objetivo de esta competencia es predecir 3 meses de datos de ventas a nivel de artículo en diferentes ubicaciones de tiendas.

0.1.3 Descripciones de archivos

• train.csv: Datos de entrenamiento

• test.csv: Datos de prueba (Nota: la división pública/privada se basa en el tiempo)

• sample_submission.csv : un archivo de envío de muestra en el formato correcto

0.1.4 Campos de información

• Fecha : Fecha de los datos de venta. No hay efectos de vacaciones ni cierres de tiendas.

• Tienda: ID de tienda

• Artículo: ID del artículo

• Ventas : Número de artículos vendidos en una tienda determinada en una fecha determinada.

0.2 Librerias de trabajo

if(!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,#Manipulación
               tidymodels,#Modelos ML
               janitor,#Limpieza
               lubridate,#Fechas
               parsnip,#modelo
               timetk,#split
               tserie,#Time series,
               forecast,#pronostico
               lmtest,#Pruebas
               tseries,#prueba adf
               seastests,#prueba estacionalidad
               rugarch,fGarch,#Modelos Garch
               FinTS,#prueba
               stats,
               zoo,#grafico
               modeltime,#Algoritmos
               openxlsx,#Excel
               fs)#Directorio

0.3 1. Datos iniciales

# Cargar bases de datos
data_train <- read_csv('data/raw/train.csv')

#Estructura datos de entrenamiento
glimpse(data_train)

## Rows: 913,000
## Columns: 4
## $ date  <date> 2013-01-01, 2013-01-02, 2013-01-03, 2013-01-04, 2013-01-05, 201…
## $ store <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ item  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ sales <dbl> 13, 11, 14, 13, 10, 12, 10, 9, 12, 9, 9, 7, 10, 12, 5, 7, 16, 7,…

sprintf("Los datos de entrenamiento tienen %d filas y %d columnas", nrow(data_train), ncol(data_train))

## [1] "Los datos de entrenamiento tienen 913000 filas y 4 columnas"

# Convertimos variables a factor
data_train$store <- factor(data_train$store)
data_train$item <- factor(data_train$item)

0.3.1 Valores perdidos

data_train %>% 
  summarise(across(everything(), ~ sum(is.na(.)))) %>% 
  pivot_longer(everything(),
               names_to = "Variables",
               values_to = "NA's")

No hay valores perdidos detectados en la serie

0.4 2. Análisis de la serie de tiempo

# Total de Ventas por tienda
data_train %>%
  # Crear variable (fecha) de agrupación
  mutate(anio_mes = floor_date(date, "month")) %>%
  # Agrupar por tienda y por mes ~ año
  group_by(store, anio_mes) %>%
  # Calcular total de ventas mensuales por tienda
  summarise(ventas_totales = sum(sales)) %>%
  ungroup() %>%
  # Crear gráfico de líneas con facetas
  ggplot(aes(x = anio_mes, y = ventas_totales)) +
  geom_line() +
  facet_wrap(~ store, nrow = 2, ncol = 5, labeller = labeller(store = function(x) paste("Store", x)))+
  labs(title = "Ventas Totales Mensuales por Tienda",
       x = " ",
       y = " ") +
  scale_y_continuous(labels = comma) +
  theme_bw()

# Total de ventas por articulo en cada tienda
g1 <- data_train %>%
  # Crear variable (fecha) de agrupación
  mutate(anio_mes = floor_date(date, "month")) %>%
  # Agrupar por tienda, artículo y mes ~ año
  group_by(store, item, anio_mes) %>%
  # Calcular total de ventas mensuales por tienda y artículo
  summarise(ventas_totales = sum(sales)) %>%
  ungroup() %>%
  # Crear gráfico de líneas con facetas
  ggplot(aes(x = anio_mes, y = ventas_totales, color = item, group = item)) +
  geom_line() +
  facet_wrap(~ store, nrow = 2, ncol = 5, labeller = labeller(store = function(x) paste("Store", x))) +
  labs(title = "Ventas Totales Mensuales por Tienda y Artículo",
       x = " ",
       y = " ") +
  scale_y_continuous(labels = comma) +
  theme_minimal() +
  theme(legend.position = "bottom", # Coloca la leyenda en la parte inferior
        legend.box = "horizontal") + # A lo largo del gráfico
  guides(color = guide_legend(nrow = 5, ncol = 10, title = "Item")) # leyenda en 5 filas y 10 columnas
plotly::ggplotly(g1)

0.4.1 Metricas de ventas

#Tabla con metricas relevantes solo para el ITEM 1
metricas <- data_train %>%
  group_by(item, store) %>%
  summarise(
    numero_ventas = n(),
    precio_promedio = mean(sales, na.rm = TRUE),
    ventas_totales = sum(sales, na.rm = TRUE)
  ) %>%
  ungroup()

# Visualizar las primeras filas del resultado
head(metricas, 10)

#Tabla con metricas relevantes solo para 5 tiendas para todos los item
metricas %>%
  group_by(item) %>%
  slice_head(n = 5)

# Alternativamente, si quieres ver todos los registros:
# print(resultado, n = Inf)

0.5 3. Preparar datos

0.5.1 i. Datos semanales

# Preparar los datos semanales
datos_preparados <- data_train %>%
  # Filtrar para incluir solo datos desde el 1 de enero de 2013
  filter(date >= as.Date("2013-01-01")) %>%
  # Añadir una columna para la semana, comenzando desde el 1 de enero de 2013
  mutate(week = floor_date(data_train$date - days(1), "week", week_start = 1) + days(1)) %>%
  # Agrupar por item, store y week (semana)
  group_by(store, item, week) %>%
  # Sumar las ventas semanales
  summarise(ventas_semanales = sum(sales, na.rm = TRUE), .groups = "drop") %>%
  # Ordenar los datos
  arrange(store, item, week)

0.5.2 ii. Datos-lista

datos_lista <- datos_preparados %>%
  group_by(store, item) %>%
  summarise(data_list = list(tibble(date = week, ventas_semanales = ventas_semanales)))

0.5.2.1 ii.i Cinco tiendas

##Seleccionar solo cinco tiendas
data_list_five <- datos_lista[c(5,56,103),]

# Gráfico de series

data_list_five %>%
  unnest(data_list) %>%
  mutate(date = as.Date(date)) %>%
  arrange(store, item, date) %>% 
  ggplot(aes(x = date, y = ventas_semanales,
             color = as.factor(store),
             group = interaction(store, item))) +
  geom_line() +
  theme_minimal() +
  labs(title = "Ventas Semanales Originales", color = "Tienda") +
  facet_wrap(~item, scales = "free_y")

#print(data_list_five$data_list[[1]])

0.5.3 iii. Convertir a serie de tiempo

# Definir la función para convertir a series de tiempo
my_ts <- function(df) {
  #datos
  ts(df$ventas_semanales,
     #Periodo de inicio
     start = c(year(min(df$date)), week(min(df$date))),
     #frecuencia semanal
     frequency = 52)
}

# Aplicar la función a los datos
data_list_five_ts <- data_list_five %>% 
  mutate(
    #aplicar función a los elementos dentro lista dentro de mi columna
    value_ts = map(data_list, my_ts)
  )

Verificamos que se haya realizado la transformación de mis datos a series de tiempo

print(data_list_five_ts$value_ts[[1]])

## Time Series:
## Start = c(2013, 1) 
## End = c(2018, 1) 
## Frequency = 52 
##   [1]  66  59  54  66  74  70  68  68  71  77  89  84  98 100  93 102  98 100
##  [19] 104 123 125 108 113 103 104 137 107 114 127 128  94  94 144 125 108  87
##  [37] 106 103 107 113  83  93  97  90  97  88 113  97  79  86  89  73  58  73
##  [55]  65  93  81  70  79  82  94  90 120 103 114 119 133 117 121 106 143 136
##  [73] 144 133 124 133 136 134 144 145 149 130 137 129 137 120 140 119 105 115
##  [91] 111 114 104 109  94 112 109 128 105 115  92  67  75  84  76  82  73  77
## [109] 101  72  89  87  99  97 107 114 102 132 114 107 115 121 129 125 126 113
## [127] 158 148 129 149 131 145 129 165 128 145 135 143 141 126 123 102 129  97
## [145] 118 102 115 111 116 100 119 148  91  80  87  92  75  91  88 100 105  98
## [163] 101 106  96 127 118 103 122 110 132 132 157 130 115 161 145 137 168 161
## [181] 135 130 177 166 164 169 162 158 143 151 154 113 146 138 144 117 134 119
## [199] 127 141 148 115 117 138 128  97  84 101 110 105  82  95  78  93 103 110
## [217] 119 110 117 104 125 147 130 153 147 139 141 149 163 141 177 166 170 145
## [235] 145 158 160 151 160 128 152 145 154 149 125 134 147 127 131 115 133 134
## [253] 125 128 134 158 110  90 110  86  87

0.5.4 iv. Desestacionalizar serie

Para poder trabajar con la serie, es necesario desestacionalizarla, como se observo en los graficos iniciales, existen una tendencia y creciente y patrones marcados a lo largo de la serie, por lo que se considera la existencia de estacionalidad.

En el siguiente código se introduce el uso de tryCatch, el cual ayuda a que si ocurre un error durante la descomposición, dentro de tryCatch se crea la función safe_stl para capturar el error y devuelver una lista con NULL en el resultado y la descomposición, e incluye el error capturado.

#library(stats)
# Función para aplicar STL de manera segura
safe_stl <- function(x) {
  tryCatch(
    {
      #Descomposición STL de la serie de tiempo
      # El componente estacional es periodico, y el algóritmo debe ser robusto a valores atípicos
      decomp <- stl(x, s.window = "periodic", robust = TRUE)
      # Descomponer componente estacional
      seasadj <- x - decomp$time.series[, "seasonal"]
      list(result = seasadj, decomposition = decomp, error = NULL)
    },
    # Código que puede generar errores
    error = function(e) {
      # Código para manejar errores,  NULL como indicador de que no hubo error.
      list(result = NULL, decomposition = NULL, error = e)
    }
  )
}

# Aplicar STL a tus series de tiempo
data_list_five_ts_stl <- data_list_five_ts %>%
  mutate(
    #aplicamos función a cada serie de tiempo 
    value_sa = map(value_ts, safe_stl)
  )

Cada elemento de value_sa es una lista que contiene:

• result: La serie ajustada estacionalmente.
• decomposition: El objeto de descomposición stl.
• error: NULL si no hubo error, o el error capturado si hubo uno.

# Tabla con valores de tendencia y residual de la serie descompuesta
#data_list_five_ts_stl$value_sa[[1]]

0.5.4.1 Gráfica Original vs Desestacionalizada

# Función para graficar la serie original y la desestacionalizada
plot_original_and_adjusted <- 
  # *data_ts*, *data_desestacionalzida*, *titulo grafico*
  function(original, adjusted, title) {
    plot(original, main = title, ylab = "Ventas", col = "blue", lty = 1)
    lines(adjusted, col = "red", lty = 2)
    legend("topleft", legend = c("Original", "Desestacionalizada"),
           col = c("blue", "red"), lty = c(1, 2))
}


#_____________ Verificación y graficación de la primera serie de tiempo __________

#Extraer la primera serie
first_series <- data_list_five_ts_stl$value_sa[[1]]

# Si existe error imprime el mensaje de error
if (!is.null(first_series$error)) {
  print(paste("Error en la primera serie:", first_series$error))

  # Si no hay procede a graficar
} else {
  plot_original_and_adjusted(
    data_list_five_ts_stl$value_ts[[1]],
    first_series$result,
    paste("Serie para Store", data_list_five_ts_stl$store[1], 
          "Item", data_list_five_ts_stl$item[1])
  )
}

# Verificar cuántas series se procesaron exitosamente

success_count <- 
  #sumar el número de TRUE(s) encontrados por la función
  sum(sapply(data_list_five_ts_stl$value_sa,
             #función anonima para identificar si es "NULL"
             function(x) is.null(x$error)))
print(paste("Series procesadas exitosamente:", success_count, "de", nrow(data_list_five_ts_stl)))

## [1] "Series procesadas exitosamente: 3 de 3"

0.5.5 v. Autoarima

Para este caso aplicamos la función auto.arima() a nuestra serie de tiempo para obtener el mejor modelo para nuestros datos, de la misma manera a como se ha realizado anteriormente, se utiliza la función tryCatch() para poder realizar la función a todas las observaciones.

#library(forecast)

# Función para aplicar auto.arima *de manera segura*

# list(result = model, error = NULL):
# Si el ajuste del modelo tiene éxito, se devuelve una lista
# con el modelo ajustado en *result*
# y *NULL* en *error* para indicar que no hubo error.

safe_auto_arima <- function(x) {
  tryCatch(
    {
    # *try*: Contiene el código que intenta ajustar el modelo ARIMA.
      model <- auto.arima(x)
      list(result = model, error = NULL)
    },
    
    #catch: Contiene el código que se ejecuta si ocurre un error
  
    #Esta función anónima captura el error "e" que ocurre
    error = function(e) {
      list(result = NULL, error = e)
    }
  )
}

# Aplicar auto.arima a tus series de tiempo
data_list_five_ts_arima <- data_list_five_ts_stl %>%
  mutate(
    arima_model = map(value_ts, safe_auto_arima)
  )

# Función para imprimir un resumen del modelo ARIMA
print_arima_summary <- 
  # model lista resultante de la función safe_auto_arima
  # store: Identificador de la tienda.
  # item: Identificador del ítem.
  function(model, store, item) {
    #Verificación de errores: (no Null)
    if (!is.null(model$error)) {
    cat("Error en Store", store, "Item", item, ":", model$error, "\n")
      #Si hay un error, imprime un mensaje con cat:
      } else {
    cat("Resumen para Store", store, "Item", item, ":\n")
    print(summary(model$result))
    cat("\n")
  }
}

#Imprimir resúmenes de los modelos ARIMA
# for (i in 1:nrow(data_list_five_ts_arima)) {
#   print_arima_summary(
#     data_list_five_ts_arima$arima_model[[i]],
#     data_list_five_ts_arima$store[i],
#     data_list_five_ts_arima$item[i]
#   )
# }

# Contar cuántos modelos se ajustaron exitosamente
success_count <- sum(sapply(data_list_five_ts_arima$arima_model, function(x) is.null(x$error)))
cat("Modelos ARIMA ajustados exitosamente:", success_count, "de", nrow(data_list_five_ts_arima), "\n")

## Modelos ARIMA ajustados exitosamente: 3 de 3

0.5.6 vi. Tabla con datos pronosticados

# Función para crear una tibble con los datos pronosticados
create_forecast_tibble <- function(ts_data, arima_model) {
  if (!is.null(arima_model$result)) {
    forecast_result <- forecast(arima_model$result, h = 52)  # Pronóstico para 52 semanas (1 año)
    
    # Obtener la fecha de inicio de la serie temporal original
    start_date <- as.Date(paste(start(ts_data)[1], 1, 1), format="%Y %U %u") + 7 * (start(ts_data)[2] - 1)
    
    # Crear fechas para el pronóstico
    last_date <- start_date + 7 * (length(ts_data) - 1)
    forecast_dates <- seq(last_date + 7, by = "week", length.out = 52)
    
    # Crear tibble con fechas y valores pronosticados
    forecast_tibble <- tibble(
      fecha = forecast_dates,
      pronostico = as.numeric(forecast_result$mean)
    )
    
    return(forecast_tibble)
  } else {
    return(NULL)
  }
}

# Aplicar función para tibbles de pronóstico
data_list_five_ts_arima <- data_list_five_ts_arima %>% 
  mutate(
    forecast_data = map2(value_ts, arima_model, create_forecast_tibble)
  )

# Imprimir resúmenes de los modelos ARIMA
# for (i in 1:nrow(data_list_five_ts_arima)) {
#   cat("Resumen para Store", data_list_five_ts_arima$store[i], 
#       "Item", data_list_five_ts_arima$item[i], ":\n")
#   if (!is.null(data_list_five_ts_arima$arima_model[[i]]$result)) {
#     print(summary(data_list_five_ts_arima$arima_model[[i]]$result))
#   } else {
#     cat("Error en el ajuste del modelo ARIMA\n")
#   }
#   cat("\n")
# }

0.5.7 vi. Gráfico-Forecast

# Prueba: gráficar la primera serie con su pronóstico
if (!is.null(data_list_five_ts_arima$arima_model[[1]]$result)) {
  first_model <- data_list_five_ts_arima$arima_model[[1]]$result
  forecast_result <- forecast(first_model, h = 52)  # Pronóstico para 52 semanas (1 año)
  
  plot(forecast_result, 
       main = paste("Pronóstico ARIMA para Store", data_list_five_ts_arima$store[1], 
                    "Item", data_list_five_ts_arima$item[1]),
       xlab = "Tiempo", ylab = "Ventas semanales")
}

# Función para crear un gráfico del modelo ARIMA
create_arima_plot <- function(ts_data, arima_model, store, item) {
  if (is.null(arima_model$error)) {
    forecast_result <- forecast(arima_model$result, h = 52)  # Pronóstico para 52 semanas (1 año)
    
    # Crear el gráfico usando ggplot2
    p <- autoplot(forecast_result) +
      ggtitle(paste("Pronóstico ARIMA para Store", store, "Item", item)) +
      xlab("Tiempo") +
      ylab("Ventas semanales") +
      theme_minimal()
    
    return(p)
  } else {
    return(NULL)
  }
}

#library(fs): Para manejar rutas de archivos de manera más eficiente


# Crear la estructura de carpetas
dir_create("img/arima_models",
           #asegurar que se creen ambas carpetas si no existen
           recurse = TRUE)

# Guardar todos los gráficos como archivos PNG
for (i in 1:nrow(data_list_five_ts_arima)) {
  #Donde es diferentes de null salta
  if (!is.null(data_list_five_ts_arima$arima_plot[[i]])) {
    #Nombre del archivo imagen
    filename <- paste0("img/arima_models/arima_plot_store_", 
                       data_list_five_ts_arima$store[i], 
                       "_item_", 
                       data_list_five_ts_arima$item[i], 
                       ".png")
    #guardar
    ggsave(
      filename = filename,
      plot = data_list_five_ts_arima$arima_plot[[i]],
      width = 10, height = 6
    )
    
    cat("Gráfico guardado:", filename, "\n") 
  }
}

0.5.8 Excel datos de los modelos

0.5.8.1 Datos pronosticados

# Crear un nuevo archivo Excel
wb <- createWorkbook()

# Función para agregar hojas al libro de trabajo con los datos de pronóstico
add_forecast_sheet <- function(forecast_data, item, store) {
  # Generar el nombre de la hoja concatenando el nombre del artículo (item) y la tienda (store)
  sheet_name <- paste("Item", item, "Store", store, sep = "_")
  addWorksheet(wb, sheet_name)
  
  # Escribir datos de pronóstico en la hoja
  writeData(wb, sheet_name, forecast_data)
}

# Aplicar la función a cada observación
data_list_five_ts_arima %>%
  select(item, store, forecast_data) %>%
  pwalk(~ add_forecast_sheet(..3, ..1, ..2))

# Guardar el archivo Excel
saveWorkbook(wb, "pronosticos_arima.xlsx", overwrite = TRUE)

0.5.8.2 Descomposición serie

El siguiente código genera una tabla por cada serie u observación en sus componentes en descomposición de la serie.

# Función para extraer la descomposición
extract_decomposition <- function(stl_result) {
  if (!is.null(stl_result$decomposition)) {
    decomposition <- stl_result$decomposition$time.series
    tibble(
      date = unique(floor_date(data_train$date - days(1), "week", week_start = 1) + days(1)),
      seasonal = decomposition[, "seasonal"],
      trend = decomposition[, "trend"],
      remainder = decomposition[, "remainder"]
    )
  } else {
    tibble(date = NA, seasonal = NA, trend = NA, remainder = NA)
  }
}

# Frame de las series stl
data_list_five_ts_arima <- data_list_five_ts_arima %>%
  mutate(
    decomposition_df = map(value_sa, extract_decomposition)
  )

# Crear un nuevo archivo Excel
wb <- createWorkbook()

# Función para agregar hojas al libro de trabajo
add_sheet_to_workbook <- function(data, item, store) {
  # Genera el nombre de la hoja concatenando el nombre del artículo (item) y la tienda (store)
  sheet_name <- paste("Item", item, "Store", store, sep = "_")
  # Agrega una nueva hoja de trabajo al libro de Excel (wb) con el nombre generado
  addWorksheet(wb, sheet_name)
  # Escribe los datos (data) en la hoja recién creada del libro de Excel (wb)
  writeData(wb, sheet_name, data)
}

# Aplicar la función add_sheet_to_workbook a cada observación en data_list_five_ts_arima
data_list_five_ts_arima %>%
  select(item, store, decomposition_df) %>%
  # 3 = decomposition_df; 1 = item; 2 = store
  pwalk(~add_sheet_to_workbook(..3, ..1, ..2))

# Guardar el archivo Excel
saveWorkbook(wb, "descomposicion_stl.xlsx", overwrite = TRUE)

0.6 i. Tidymodels

0.6.1 Familia de funciones:

0.6.1.1 Para el flujo de trabajo

Este funciona para solo una observación, es decir la store 1 para el producto 1.

# Función para crear un flujo de trabajo de modelo seguro
create_workflow_safe <- function(model_spec, model_name) {
  tryCatch({
    workflow() %>%
      add_model(model_spec) %>%
      add_recipe(recipe(ventas_semanales ~ date, data = train_data))
  }, error = function(e) {
    message(paste("Error in", model_name, "workflow:", e$message))
    NULL
  })
}

0.6.1.2 Función para ajustar el modelo

# Función para ajustar el modelo de manera segura
fit_model_safe <- function(workflow, train_data, model_name) {
  tryCatch({
    workflow %>%
      fit(train_data)
  }, error = function(e) {
    message(paste("Error in", model_name, "model fit:", e$message))
    NULL
  })
}

0.6.2 División de datos: Entrenamiento y prueba

# Datos de entrenamiento y prueba
train_data <- data_list_five$data_list[[1]] %>%
  slice(1:(n() - 52))
test_data <- data_list_five$data_list[[1]] %>%
  slice((n() - 51):n())

0.6.3 Especificación de los modelos

# Especificaciones de los modelos para datos semanales
model_specs <- list(
  model_boost = arima_boost(
    seasonal_period = 52,  # Para datos semanales
    learn_rate = 0.015, 
    min_n = 2
  ) %>% 
    set_engine("auto_arima_xgboost",
               seasonal = TRUE,  # Forzar la consideración de componentes estacionales
               stepwise = FALSE,  # Búsqueda más exhaustiva
               approximation = FALSE) %>% 
    set_mode("regression"),
  
  model_arima = arima_reg(
    seasonal_period = 52  # Para datos semanales
  ) %>%
    set_engine("auto_arima",
               stepwise = FALSE,  # Realiza una búsqueda más exhaustiva
               approximation = FALSE,  # Usa cálculos exactos de verosimilitud
               seasonal = TRUE) %>%  # Fuerza la inclusión de términos estacionales
    set_mode("regression"),
  
  model_ets = exp_smoothing(
    error = "multiplicative", 
    trend = "multiplicative", 
    season = "multiplicative"
  ) %>%
    set_engine("ets") %>%
    set_mode("regression"),
  
  model_prophet = prophet_reg(
    seasonality_yearly = TRUE,
    seasonality_weekly = TRUE,
    seasonality_daily = FALSE  # Los datos son semanales, así que la estacionalidad diaria no es relevante
  ) %>%
    set_engine("prophet",
               yearly_seasonality = 'auto',
               weekly_seasonality = 'auto',
               daily_seasonality = FALSE) %>%
    set_mode("regression")
  
)

0.6.4 Aplicación de funciones: Flujo de trabajo y Ajuste de Modelos

Aplicamos las especificaciones de los modelos a los flujos de trabajo.

# Crear workflows de manera segura
workflows <- map2(model_specs, names(model_specs), create_workflow_safe) %>%
  compact()

Posteriormente ajustamos los modelos de nuestro flujo de trabajo.

# Ajustar los modelos de manera segura
fitted_models <- map2(workflows, names(workflows), fit_model_safe, train_data = train_data) %>%
  compact()

0.6.5 Calibrar modelos

Evaluamos nuestos modelos ajustados.

# Calibrar los modelos con datos de prueba
calibration_table <- map_dfr(fitted_models, function(model) {
  tryCatch({
    modeltime_table(model) %>%
      modeltime_calibrate(test_data)
  }, error = function(e) {
    message(paste("Error in calibration:", e$message))
    NULL
  })
}, .id = "model_name")

0.6.5.1 Mejor modelo ajustado

# Evaluar el rendimiento de los modelos
calibration_table %>%
  modeltime_accuracy() %>% 
  arrange(rmse)

# Seleciona solo el mejor modelo
best_model <- calibration_table %>%
  modeltime_accuracy() %>% 
  arrange(rmse) %>%
  slice(1)

# Obtener el mejor modelo ajustado desde calibration_table
best_model_fit <- calibration_table %>%
  filter(.model_desc == best_model$.model_desc)

# Reajustar solo el mejor modelo a todos los datos
refit_table <- best_model_fit %>%
  modeltime_refit(data = train_data)

0.6.5.2 Pronóstico

# Realizar el pronóstico a 52 semanas con el mejor modelo
future_data <- data_list_five$data_list[[1]] %>%
  slice((n() - 51):n()) %>%
  mutate(date = seq(max(date) + 7, by = "week", length.out = 52))

forecast <- refit_table %>%
  modeltime_forecast(new_data = future_data, actual_data = data_list_five$data_list[[1]])

# Mostrar el pronóstico
forecast %>%
  plot_modeltime_forecast()

0.7 ii. Tidymodels

0.7.1 3. Mejor modelo estimado para cada serie

# 1. Preparación de los datos
expanded_data <- data_list_five %>%
  unnest(data_list) %>%
  mutate(date = as.Date(date)) %>%
  arrange(store, item, date)

# 2. Función para procesar cada grupo (store-item)
process_group <- function(data) {
  # División de los datos
  split_data <- time_series_split(data, date, assess = "52 weeks", cumulative = TRUE)
  
  # Definición de modelos
  model_spec_arima <- arima_reg(seasonal_period = 52) %>%
    set_engine("auto_arima", stepwise = FALSE, approximation = FALSE, seasonal = TRUE) %>%
    set_mode("regression")
  
  model_spec_xgboost <- arima_boost(
    seasonal_period = 52,
    learn_rate = 0.015,
    min_n = 2) %>%
    set_engine("auto_arima_xgboost", seasonal = TRUE, stepwise = FALSE, approximation = FALSE) %>%
    set_mode("regression")
  
  model_spec_prophet <- prophet_reg(
    seasonality_yearly = TRUE, seasonality_weekly = TRUE, seasonality_daily = FALSE) %>%
    set_engine("prophet", yearly_seasonality = 'auto', weekly_seasonality = 'auto', daily_seasonality = FALSE) %>%
    set_mode("regression")
  
  # Entrenamiento de modelos
  workflow_arima <- workflow() %>%
    add_model(model_spec_arima) %>%
    add_recipe(recipe(ventas_semanales ~ date, data = training(split_data))) %>%
    fit(training(split_data))
  
  workflow_xgboost <- workflow() %>%
    add_model(model_spec_xgboost) %>%
    add_recipe(recipe(ventas_semanales ~ date, data = training(split_data))) %>%
    fit(training(split_data))
  
  workflow_prophet <- workflow() %>%
    add_model(model_spec_prophet) %>%
    add_recipe(recipe(ventas_semanales ~ date, data = training(split_data))) %>%
    fit(training(split_data))
  
  # Calibración y evaluación
  calibration_table <- modeltime_table(workflow_arima, workflow_xgboost, workflow_prophet) %>%
    modeltime_calibrate(new_data = testing(split_data))
  
  # Selección del mejor modelo
  best_model <- calibration_table %>%
    modeltime_accuracy() %>%
    arrange(mape) %>%
    slice(1)
  
  # Reajuste del mejor modelo
  refit_table <- calibration_table %>%
    filter(.model_id == best_model$.model_id) %>%
    modeltime_refit(data = data)
  
  # Generación de pronósticos
  future_data <- data %>%
    future_frame(.date_var = date, .length_out = "52 weeks")
  
  forecast_table <- refit_table %>%
    modeltime_forecast(new_data = future_data, actual_data = data)
  
  # Diagnóstico de residuos
  residuals <- refit_table %>%
    modeltime_residuals()
  
  return(list(best_model = best_model, forecast = forecast_table, residuals = residuals))
}

# 3. Aplicar el proceso a cada grupo
results <- expanded_data %>%
  group_by(store, item) %>%
  nest() %>%
  mutate(model_results = map(data, process_group)) %>%
  ungroup()

# 4. Extraer los mejores modelos, pronósticos y residuos
best_models <- results %>%
  mutate(best_model = map(model_results, ~.x$best_model)) %>%
  select(store, item, best_model) %>%
  unnest(best_model)

forecasts <- results %>%
  mutate(forecast = map(model_results, ~.x$forecast)) %>%
  select(store, item, forecast) %>%
  unnest(forecast)

residuals <- results %>%
  mutate(residuals = map(model_results, ~.x$residuals)) %>%
  select(store, item, residuals) %>%
  unnest(residuals)


# 5. Visualizar los pronósticos
forecasts %>%
  ggplot(aes(x = .index, y = .value, color = .key)) +
  geom_line() +
  facet_grid(item ~ store, scales = "free_y") +
  theme_minimal() +
  labs(title = "Pronósticos por Tienda e Ítem",
       x = "Fecha", y = "Ventas Semanales",
       color = "Serie") +
  scale_color_manual(values = c("actual" = "black", "prediction" = "red"))

# 6. Diagnóstico de residuos
plot_residuals <- function(res_data) {
  res_data %>%
    plot_modeltime_residuals(.interactive = FALSE)
}
residual_plots <- residuals %>%
  group_by(store, item) %>%
  nest() %>%
  mutate(plots = map(data, ~plot_residuals(.x)))

# Mostrar los gráficos de residuos
print(residual_plots$plots)

## [[1]]

## 
## [[2]]

## 
## [[3]]

# 7. Visualización de datos originales
ggplot(expanded_data, aes(x = date, y = ventas_semanales, color = as.factor(store), group = interaction(store, item))) +
  geom_line() +
  theme_minimal() +
  labs(title = "Ventas Semanales Originales", color = "Tienda") +
  facet_wrap(~item, scales = "free_y")

# 8. Resumen de los mejores modelos
best_models %>%
  group_by(.model_desc) %>%
  summarise(count = n()) %>%
  mutate(percentage = count / sum(count) * 100)

# Primero, creamos un dataframe con la información del modelo para cada combinación de store e item
model_info <- forecasts %>%
  filter(.key == "prediction") %>%
  group_by(store, item) %>%
  summarise(model = first(.model_desc), .groups = "drop")

# Luego, unimos esta información con el dataframe original de forecasts
forecasts_with_model <- forecasts %>%
  left_join(model_info, by = c("store", "item"))

# Ahora, creamos una nueva columna que combine store, item y model
forecasts_with_model <- forecasts_with_model %>%
  mutate(facet_label = paste("Tienda", store, "Ítem", item, "\n", model))

# Finalmente, creamos el gráfico usando los datos modificados
graf_01 <- forecasts_with_model %>%
  ggplot(aes(x = .index, y = .value)) +
  geom_line(aes(color = .key)) +
  geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi), fill = "blue", alpha = 0.2) +
  facet_wrap(~ facet_label, scales = "free_y") +
  theme_minimal() +
  labs(title = "Pronósticos por Tienda e Ítem",
       x = "Fecha", y = "Ventas Semanales",
       color = "Serie") +
  scale_color_manual(values = c("actual" = "black", "prediction" = "red")) +
  theme(strip.text = element_text(size = 8),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  coord_cartesian(clip = "off")


ggsave("img/tidymodels/pronosticos_tienda_item.png", plot = graf_01,
       width = 10, height = 6,dpi = 300,bg = "white")

graf_01

##_________________________________________________________________
#Si quisiera cada uno de los gráficos creados arriba imprimo lo siguiente:

# Crea una función para generar un gráfico individual
create_individual_plot <- function(data, facet_label) {
  ggplot(data, aes(x = .index, y = .value)) +
    geom_line(aes(color = .key)) +
    geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi), fill = "blue", alpha = 0.2) +
    theme_minimal() +
    labs(title = paste("Pronóstico para", facet_label),
      x = "Fecha", y = "Ventas Semanales",
      color = "Serie",
      caption = "Elaborado por @EspacioColectivoMx") +
    scale_color_manual(values = c("actual" = "black", "prediction" = "red")) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    coord_cartesian(clip = "off")
}

# Crea una carpeta para guardar los gráficos si no existe
dir.create("img/tidymodels", recursive = TRUE, showWarnings = FALSE)

# Crea y guarda un gráfico para cada facet_label única
forecasts_with_model %>%
  group_by(facet_label) %>%
  group_walk(~ {
    p <- create_individual_plot(.x, .y$facet_label)

    # Crea un nombre de archivo válido
    filename <- paste0("pronostico_", gsub("[^a-zA-Z0-9]", "_", .y$facet_label), ".png")

    ggsave(paste0("img/tidymodels/", filename),
           plot = p,
           width = 10,
           height = 6,
           dpi = 300,
           bg = "white")
  })

# Crear una lista de gráficos
plot_list <- forecasts_with_model %>%
  group_by(facet_label) %>%
  group_map(~ create_individual_plot(.x, .y$facet_label))

# Imprimir los gráficos
for (plot in plot_list) {
  print(plot)
}

0.8 Código alternativo

0.8.1 1. Modelo de código

Este esta anidado a una función que estima todos los modelos generados para las series, por lo cual es lento y de mayor coste.

# # library(tidyverse)
# # library(modeltime)
# # library(timetk)
# # library(purrr)
# 
# # 1. Preparación de los datos (sin unnest)
# process_data <- function(data) {
#   data %>%
#     mutate(date = as.Date(date)) %>%
#     arrange(date)
# }
# 
# data_list_five <- data_list_five %>%
#   mutate(data_list = map(data_list, process_data))
# 
# # 2. División de los datos
# split_data <- function(data) {
#   time_series_split(data, date, assess = "52 weeks", cumulative = TRUE)
# }
# 
# # 3. Definición de modelos (igual que antes)
# model_spec_arima <- arima_reg(
#   seasonal_period = 52) %>% # Para datos semanales
#   set_engine("auto_arima",
#              stepwise = FALSE, # Realiza una búsqueda más exhaustiva
#              approximation = FALSE, # Usa cálculos exactos de verosimilitud
#              seasonal = TRUE) %>% # Fuerza la inclusión de términos estacionales
#   set_mode("regression")
# 
# model_spec_xgboost <- arima_boost(
#   seasonal_period = 52,
#   learn_rate = 0.015,
#   min_n = 2) %>%
#   set_engine("auto_arima_xgboost",
#              seasonal = TRUE,
#              stepwise = FALSE,
#              approximation = FALSE) %>%
#   set_mode("regression")
# 
# model_spec_prophet <- prophet_reg(
#   seasonality_yearly = TRUE,
#   seasonality_weekly = TRUE,
#   seasonality_daily = FALSE) %>%
#   set_engine("prophet",
#              yearly_seasonality = 'auto',
#              weekly_seasonality = 'auto',
#              daily_seasonality = FALSE) %>%
#   set_mode("regression")
# 
# # 4-11. Función para entrenar modelos, evaluar y pronosticar
# model_and_forecast <- function(data) {
#   split <- split_data(data)
#   
#   # Entrenar modelos
#   workflow_arima <- workflow() %>%
#     add_model(model_spec_arima) %>%
#     add_recipe(recipe(ventas_semanales ~ date, data = training(split))) %>%
#     fit(training(split))
#   
#   workflow_xgboost <- workflow() %>%
#     add_model(model_spec_xgboost) %>%
#     add_recipe(recipe(ventas_semanales ~ date, data = training(split))) %>%
#     fit(training(split))
#   
#   workflow_prophet <- workflow() %>%
#     add_model(model_spec_prophet) %>%
#     add_recipe(recipe(ventas_semanales ~ date, data = training(split))) %>%
#     fit(training(split))
#   
#   # Calibrar y evaluar modelos
#   calibration_table <- modeltime_table(
#     workflow_arima,
#     workflow_xgboost,
#     workflow_prophet
#   ) %>%
#     modeltime_calibrate(new_data = testing(split))
#   
#   # Seleccionar el mejor modelo
#   best_model <- calibration_table %>%
#     modeltime_accuracy() %>%
#     arrange(mape) %>%
#     slice(1)
#   
#   # Reajustar el mejor modelo
#   refit_table <- calibration_table %>%
#     modeltime_refit(data = data)
#   
#   # Generar pronósticos
#   future_data <- data %>%
#     future_frame(
#       .date_var = date,
#       .length_out = "52 weeks"
#     )
#   
#   forecast <- refit_table %>%
#     modeltime_forecast(
#       new_data = future_data,
#       actual_data = data
#     )
#   
#   return(list(calibration = calibration_table, 
#               best_model = best_model, 
#               forecast = forecast))
# }
# 
# # Aplicar el proceso a cada serie de tiempo
# results <- data_list_five %>%
#   mutate(model_results = map(data_list, model_and_forecast))
# 
# # Ahora puedes acceder a los resultados para cada serie
# # Por ejemplo, para ver el pronóstico del primer item:
# results$model_results[[1]]$forecast %>%
#   plot_modeltime_forecast(.interactive = TRUE)
# ```
# 
# ```{r}
# # Extraer las métricas de todos los modelos para todas las series
# all_metrics <- results %>%
#   mutate(metrics = map(model_results, function(x) {
#     x$calibration %>%
#       modeltime_accuracy() %>%
#       mutate(store = cur_group()$store,
#              item = cur_group()$item)
#   })) %>%
#   select(store, item, metrics) %>%
#   unnest(metrics)
# 
# # Ver las métricas
# print(all_metrics)
# 
# # Si quieres un resumen más compacto, puedes hacer algo como esto:
# summary_metrics <- all_metrics %>%
#   group_by(store, item, .model_id, .model_desc) %>%
#   summarise(across(c(mae, mape, mase, smape, rmse, rsq), mean))
# 
# print(summary_metrics)
# 
# # Para ver las métricas del mejor modelo para cada serie:
# best_models <- results %>%
#   mutate(best_model_metrics = map(model_results, ~.x$best_model)) %>%
#   select(store, item, best_model_metrics) %>%
#   unnest(best_model_metrics)
# 
# print(best_models)
# 
# # Comparar MAPE entre modelos
# ggplot(all_metrics, aes(x = .model_desc, y = mape, color = .model_desc)) +
#   geom_boxplot() +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
#   labs(title = "Comparación de MAPE entre modelos",
#        x = "Modelo", y = "MAPE")
# 
# # Ver qué modelo es el mejor más frecuentemente
# best_models %>%
#   count(.model_desc) %>%
#   arrange(desc(n))
# ```
#