Pronostico de Ventas de Cerveza

0.1 Introducción

En el ámbito del análisis de series de tiempo, un aspecto fundamental es la capacidad de prever y entender patrones en datos históricos. En este ejercicio, nos enfocamos en la base de datos drinks, que detalla las ventas de bebidas alcohólicas, específicamente bajo la serie denominada “Merchant Wholesalers, Except Manufacturers’ Sales Branches and Offices Sales: Nondurable Goods: Beer, Wine, and Distilled Alcoholic Beverages Sales”, de la FRED. El objetivo principal es aplicar técnicas de series de tiempo para realizar estimaciones a futuro, identificar los meses con mayores ventas y determinar el valor promedio en esos periodos destacados. Además, se pretende estimar las ventas futuras en los meses de mayor demanda identificados. Este análisis no solo permitirá proyectar tendencias futuras, sino también comprender los factores estacionales que influyen en la serie.

Acá el link: https://fred.stlouisfed.org/series/S4248SM144NCEN

0.1.1 Datos

Los datos de este artículo son ventas de bebidas alcohólicas originalmente del sitio web del Banco de la Reserva Federal de St. Louis data(“drinks”)

0.1.2 Librerias

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

0.1.2.1 API FREDR

#library(fredr)
key <- "914908a8840e2583e5e560d95d78c353"
fredr_set_key(key)
dta <- fredr(
  series_id = "S4248SM144NCEN",
  observation_start = as.Date("1992-01-01"),
  observation_end = Sys.Date()
)
glimpse(dta)

## Rows: 389
## Columns: 5
## $ date           <date> 1992-01-01, 1992-02-01, 1992-03-01, 1992-04-01, 1992-0…
## $ series_id      <chr> "S4248SM144NCEN", "S4248SM144NCEN", "S4248SM144NCEN", "…
## $ value          <dbl> 3459, 3458, 4002, 4564, 4221, 4529, 4466, 4137, 4126, 4…
## $ realtime_start <date> 2024-07-10, 2024-07-10, 2024-07-10, 2024-07-10, 2024-0…
## $ realtime_end   <date> 2024-07-10, 2024-07-10, 2024-07-10, 2024-07-10, 2024-0…

0.1.3 Objetivo

Estimar las ventas proyectadas

0.1.4 Proceso

Descomposición de la serie
Análisis de la serie
Pruebas
Ajuste
Construcción de modelo de estimación
Comparación de modelos
Proyección a futuro
Visualización
Guardar datos

0.1.5 Notas técnicas

0.1.5.1 a) Composición de una serie de tiempo

Tendencia secular: Es el movimiento general a largo plazo de la serie. Puede ser ascendente, descendente o horizontal. La tendencia secular refleja los cambios estructurales o de fondo que afectan a la variable a lo largo del tiempo, como el crecimiento económico, el cambio tecnológico o los cambios demográficos.
Variación estacional: Son las fluctuaciones regulares que se repiten en un período fijo, como un año, un trimestre o un mes. La estacionalidad se debe a factores como las estaciones del año, las vacaciones o los días festivos.
Variación cíclica: Son oscilaciones de mayor duración que la estacionalidad, típicamente de dos a diez años. Los ciclos económicos son un ejemplo de variación cíclica.
Variación irregular: Son fluctuaciones aleatorias o impredecibles en la serie de tiempo. La variación irregular puede ser causada por eventos aleatorios,

0.1.5.2 Modelos ARIMA

AR: Regresión sobre el mismo
I: Raices unitarias. Diferencias
MA: Errores. Tiene que ver con los errores
S: Estacionalidad. Fenomenos que se repiten

\[ Parametros (p,d,q) \]

\(p\) orden
\(i\) diferencia
\(q\) nivel de media movil

0.2 A) Serie Tradicional

#Convertir serie de tiempo
ts_data <- ts(dta$value, start = c(1992,1), frequency = 12)

0.2.1 1.Visualizaación de la Serie

#Con ggplot
dta %>% 
  ggplot(aes(y=value, x = date)) + 
  geom_line(color = "#581845", alpha = 0.7)+
  scale_y_continuous(labels = comma)+
  labs(title = "Ventas de cerveza",
       x = "",
       subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx")

#Con ts
#plot(ts_data)

0.2.2 2.Descomposición de la serie

descomposicion <- stl(ts_data, s.window = "periodic")
autoplot(descomposicion)

Se observa: - Tendencia creciente - Irregularidad de los residuos - Patrones estacionales

0.2.3 3. Pruebas

0.2.3.1 i. Autocorrelación

Detectar autorcorrelación nos puede indicar la existencia de no estacionariedad.

par(mfrow=c(1,2), mar=c(4,4,4,1)+.1)
acf(ts_data, main="Serie no Estacionaria")
pacf(ts_data)

ACF - Como ve, la inercia temporal es muy grande, el decaimiento lento sugiere ser un proceso no estacionario en media ya, por lo que puede existir relación entre los valores de la serie en el tiempo actual y los valores rezagados. (AR)

PACF - Correlación entre los valores de la serie en el tiempo actual y los valores rezagados, controlando por la autocorrelación en rezagos anteriores.(MA)

Podemos deducidr que es en los modelos MA(q), ocurre que la ACF tiene palos distintos de cero, mientras que la PACF muestra un decaimiento rápido.

Es muy común que la desviación típica de una serie temporal crezca con su nivel (o tendencia). Si pasa eso, es que hay una relación multplicativa entre la tendencia de la serie y el error.

0.2.3.2 ii. Estacionariedad

# 5. Prueba de estacionariedad (Test de Dickey-Fuller)
#library(tseries)
adf.test(ts_data)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_data
## Dickey-Fuller = -5.3346, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

#Prueba estacionariedad
tseries::pp.test(ts_data)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  ts_data
## Dickey-Fuller Z(alpha) = -298.63, Truncation lag parameter = 5, p-value
## = 0.01
## alternative hypothesis: stationary

Por ahora, podemos decir que los resultados anteriores sugieren que la serie de tiempo es estacionaria, se tiene que la prueba de hipótesis con una H0 = No Estacionaria, los valores al ser menor al valor critico de 0.05 rechazan la HO por lo tanto la serie es estacionaria, sin embargo, al contrastar con los gráficos ACF Y APCF la serie aún podría tener algunos componentes no estacionarios o estacionales.

0.2.3.3 iv. Estacionalidad

# Dado que los datos son mensuales
month <- factor(cycle(ts_data))
kruskal.test(ts_data ~ month)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  ts_data by month
## Kruskal-Wallis chi-squared = 23.28, df = 11, p-value = 0.01613

El valor bajo de la prueba Kruskal-Wallis sugiere estacionalidad.

Diferencias de la Serie para convertirla estacionaria*

#Library(forecast)
#Diferenciar la serie
serie.dif <- diff(ts_data)
#Serie estacionaria
plot(serie.dif, main="Ventas de cerveza", ylab="Total ventas", xlab="")

#Necesita más diferencias?
print(sprintf("Diferencias faltantes de la nueva serie: %d", ndiffs(serie.dif)))

## [1] "Diferencias faltantes de la nueva serie: 0"

par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot(ts_data, ylab = "ventas")
acf(ts_data, main="Serie no Estacionaria")
plot(serie.dif, ylab = "ventas")
acf(serie.dif, main="Serie Estacionaria")

A pesar de haber realizado una diferencia a la serie, se sigue teniendo presenta de retardos en el ACF.

0.2.4 4.Modelo Sin diferencia

# Ajustar un modelo SARIMA
sarima_model <- auto.arima(ts_data, seasonal=TRUE)
summary(sarima_model)

## Series: ts_data 
## ARIMA(4,1,2)(0,1,2)[12] 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2     sma1     sma2
##       0.3038  -0.1918  0.0890  -0.4902  -1.1484  0.4944  -0.3328  -0.3017
## s.e.  0.1424   0.0671  0.0633   0.0631   0.1543  0.1563   0.0627   0.0548
## 
## sigma^2 = 164212:  log likelihood = -2790.79
## AIC=5599.57   AICc=5600.07   BIC=5634.94
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 11.83396 394.1408 290.6193 0.03316235 3.213084 0.6029777
##                     ACF1
## Training set 0.007996391

# Pronóstico con SARIMA
sarima_forecast <- forecast(sarima_model, h=24)
plot(sarima_forecast)

0.2.4.1 1. Residuos

# Para la prueba de autocorrelación
#library(lmtest)
# Obtener los residuos del modelo
residuals <- residuals(sarima_model)

# Aucorrelación residual
Box.test(residuals, type="Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 0.025066, df = 1, p-value = 0.8742

El Box-Ljung test sugiere que no hay autocorrelación significativa en los residuos.

forecast::checkresiduals(sarima_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,1,2)(0,1,2)[12]
## Q* = 120.11, df = 16, p-value < 0.00000000000000022
## 
## Model df: 8.   Total lags used: 24

Hay un mayor control en los residuos anque aun existen perturbaciones.

0.2.4.2 2. Patrones estacionales

monthplot(ts_data, main="Gráfico Estacional", ylab="Valor")

líneas que tienden a estar consistentemente altas o bajas para ciertos meses, esto indica un patrón estacional, por ejemplo, puede deberse a que las ventas siempre son altas en diciembre por las fiestas de sembrinas y bajas en enero, posterior a las fiestas.

#library(seastests)
isSeasonal(ts_data, test = "qs")

## [1] TRUE

Este resultado proporciona evidencia estadística que respalda la presencia de patrones estacionales en los datos, confirmando lo que observado en el monthplot.

acf(residuals)

pacf(residuals)

qqnorm(residuals)
qqline(residuals)

0.2.5 5. GARCH

0.2.5.1 1. Especificar el Modelo ARIMA-GARCH

Supuestos Modelo 1: - Los coeficientes ARMA y GARCH están interrelacionados. - Capturar interacciones entre la dinámica de la media y la varianza. - Más complejo, puede ser más difícil de estimar.

Se sospecha que la dinámica de la media y la varianza están fuertemente relacionadas.

#library(rugarch)
#library(fGarch)
#library(FinTS)

# Especificar el modelo 1 
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                   mean.model = list(armaOrder = c(3, 1), include.mean = TRUE),
                   distribution.model = "norm")

# Esta parte define la especificación del modelo GARCH
# - variance.model: Especifica un modelo GARCH estándar (sGARCH) de orden (1,1)
# - mean.model: Especifica un modelo ARMA(3,1) para la media condicional
# - distribution.model: Asume una distribución normal para los errores

Supuestos Modelo 2: - Los coeficientes GARCH son independientes del modelo SARIMA inicial. - Separar la modelación de la media y la varianza. - Potencialmente más simple, puede ser más estable en la estimación. - La estructura de la media está bien capturada por SARIMA y solo la varianza residual necesita modelarse.

# Especificar el modelo 2
garch_spec <- ugarchspec(
  variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
  mean.model = list(armaOrder = c(0, 0), include.mean = TRUE),
  distribution.model = "norm"
)
# - variance.model: Especifica un modelo GARCH estándar (sGARCH) de orden (1,1)
# - mean.model: Especifica un modelo ARMA(0,0) con la media condicional
# - distribution.model: Asume una distribución normal para los errores

0.2.5.2 2. Ajuste de Modelo

##____________________ Modelo 1 _________________
#Residuos SARIMA
residuals_sarima <- residuals(sarima_model)


#### Ajustar el modelo
sarima_garch_model <- ugarchfit(spec = spec, data = diff(ts_data, lag=12))


# Esta línea ajusta el modelo GARCH a los datos
# - spec: Usa la especificación definida arriba
# - data: Aplica el modelo a la diferencia estacional de orden 12 de ts_data

#Residuos del Modelo
residuals_std <- residuals(sarima_garch_model, standardize=TRUE)

##____________________ Modelo 2 _________________

#Residuos SARIMA
residuals_sarima <- residuals(sarima_model)

# Ajustar el modelo GARCH a los residuos del SARIMA
garch_model <- ugarchfit(spec = garch_spec, data = residuals_sarima)

#Residuos del Modelo
residuals_std_2 <- residuals(garch_model, standardize=TRUE)

# Test Modelo 1
ArchTest(residuals_std)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuals_std
## Chi-squared = 8.229, df = 12, p-value = 0.767

# Test Modelo 2
ArchTest(residuals_std_2)

## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  residuals_std_2
## Chi-squared = 5.6168, df = 12, p-value = 0.9342

La ausencia de efectos ARCH sugiere que la serie no contiene períodos de volatilidad. Siendo la condición de que hay uno o más puntos de datos en una serie para los cuales la varianza del término de error actual o innovación es una función de los tamaños reales de los términos de error de los períodos de tiempo anteriores.

0.2.5.3 3. Mejor Modelo

### Por criterios

# Para el Modelo 1
info_criteria_1 <- infocriteria(sarima_garch_model)

# Para el Modelo 2
info_criteria_2 <- infocriteria(garch_model)

# Comparar
print(info_criteria_1)

##                      
## Akaike       14.75455
## Bayes        14.83799
## Shibata      14.75367
## Hannan-Quinn 14.78767

print(info_criteria_2)

##                      
## Akaike       14.51496
## Bayes        14.55571
## Shibata      14.51475
## Hannan-Quinn 14.53112

En conclusión:

Seleccionamos el modelo 2, el cual obtiene valores más bajos para cada criterio evaluado.

#plot(garch_model)
ugarchforecast(garch_model, n.ahead=10)

## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 10
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=may. 2024]:
##      Series Sigma
## T+1   24.26 446.6
## T+2   24.26 447.9
## T+3   24.26 449.3
## T+4   24.26 450.6
## T+5   24.26 452.0
## T+6   24.26 453.3
## T+7   24.26 454.6
## T+8   24.26 456.0
## T+9   24.26 457.3
## T+10  24.26 458.6

0.2.5.4 3. Pronostico

# Paso 1: Generar pronósticos

#Meses a pronosticar
forecast_length <- 24

# Pronósticos del modelo SARIMA
sarima_forecast <- forecast(sarima_model, h = forecast_length)

# Predicciones del modelo GARCH
garch_forecast <- ugarchforecast(garch_model, n.ahead = forecast_length)
garch_pred <- as.numeric(sigma(garch_forecast))

# Combinar las predicciones de SARIMA y GARCH
combined_forecast <- sarima_forecast$mean + garch_pred


# Generar las fechas de las predicciones
prediction_dates <- seq(ceiling_date(max(dta$date, na.rm = TRUE), "month"),
                        by = "month",
                        length.out = forecast_length)

# Crear data frames para gráficos
#Predicciones
predicted_df <- data.frame(
  time = prediction_dates,
  predicted = combined_forecast
)

# Paso 2: Crear el gráfico
#plot(garch_forecast)

#Gráfico de predicciones
gprediccion <- 
  ggplot() +
  geom_line(data = dta, aes(x = date, y = value, color = "Orginal")) +
  geom_line(data = predicted_df, aes(x = time, y = predicted, color = "Predicción")) +
  labs(title = "Predicciones del Modelo SARIMA-GARCH",
       x = " ",
       y = "Valor",
       subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx") +
  theme_minimal()+
  scale_y_continuous(labels = comma)+
  scale_color_manual(values = c("Original" = "#1A5276",
                                "Predicción" = "#D80845")) +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.title = element_text(hjust = 0.5),
        legend.box = "horizontal") +
  guides(color = guide_legend(title = NULL, title.position = "top", title.hjust = 0.5))

ggplotly(gprediccion)

0.3 B) Tidymodels

# 1. Transformar data
# No es necesario transformar la data a serie de tiempo para tidymodels.
# Serie de tiempo para algunas visualizaciones y desestacionalizar
drinks_ts <- ts(dta$value, start = c(1992,1), frequency = 12)

# 2. Partición de datos
set.seed(123)
splits <- time_series_split(dta,
                            #tamaño de datos prueba
                            assess = "12 months", 
                            #establecer un unico segmento de datos
                            cumulative = TRUE)

training(splits) %>% mutate(type="train") %>% 
  bind_rows(testing(splits) %>% mutate(type = "test")) %>% 
  ggplot(aes(x=date, y=value, color = type))+
  geom_line()

autoplot(drinks_ts)+
  labs(title = "Ventas Mensuales de Bebidas Alcohólicas", y = "Ventas", x = "Fecha")

De nuevo, se puede notar en los picos largos un comportamiento sistemático que se presenta al interior de cada año y se repite.

0.3.0.0.1 Box Plot la estacionalidad

Recordemos que los resultados obtenidos para analizar la serie, nos arroja que esta no es estacinaria, presenta tendencia y estacionalidad.

0.3.0.0.2 Prueba KPSS

\[H0: \text{El proceso es estacionario en tendencia} \\Ha: \text{La serie tiene una raíz unitaria (la serie no es estacionaria) }\]

kpss.test(drinks_ts)

## 
##  KPSS Test for Level Stationarity
## 
## data:  drinks_ts
## KPSS Level = 6.3482, Truncation lag parameter = 5, p-value = 0.01

0.3.1 Modelaje

# Modelo Prophet
model_prophet <- prophet_reg(
  seasonality_yearly = TRUE,
  seasonality_weekly = TRUE,
  seasonality_daily = TRUE
) %>%
  set_engine("prophet", 
             yearly_seasonality = 'auto',
             weekly_seasonality = 'auto',
             daily_seasonality = 'auto') %>%
  fit(value ~ date, data = training(splits))

# Modelo ARIMA
model_arima <- arima_reg(
  seasonal_period = 12  # Para datos mensuales
) %>%
  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
  fit(value ~ date, data = training(splits))

# Modelo ETS (Exponential Smoothing)
model_ets <- exp_smoothing(
  seasonal_period = 12,  # Para datos mensuales
  error = "multiplicative", 
  trend = "multiplicative", 
  season = "multiplicative"
) %>%
  set_engine("ets") %>%
  fit(value ~ date, data = training(splits))

# ARIMA Boost
model_boost <- arima_boost(
  seasonal_period = 12,  # Para datos mensuales
  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) %>% 
  fit(value ~ date + as.numeric(date) + factor(month(date, label = TRUE)), 
      data = training(splits))

# Crear una tabla de modelos
models_tbl <- modeltime_table(
  model_prophet,
  model_arima,
  model_ets,
  model_boost
)

# Paso 10: Calibrar los modelos
calibration_table <- models_tbl %>%
  modeltime_calibrate(new_data = testing(splits))

# Paso 11: Evaluar el rendimiento de los modelos
calibration_table %>%
  modeltime_accuracy()%>% 
  table_modeltime_accuracy(.interactive = FALSE)

.model_id	.model_desc	.type	mae	mape	mase	smape	rmse	rsq
Accuracy Table
1	PROPHET	Test	1066.94	7.47	0.60	7.02	1340.79	0.89
2	ARIMA(0,1,1)(2,1,1)[12]	Test	475.15	3.16	0.27	3.10	572.66	0.91
3	ETS(M,M,M)	Test	895.82	5.84	0.51	5.60	1078.97	0.86
4	ARIMA(0,1,1)(2,1,1)[12] W/ XGBOOST ERRORS	Test	461.52	3.05	0.26	3.01	553.95	0.90

# Paso 12: Visualizar las predicciones
calibration_table %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = dta
  ) %>%
  plot_modeltime_forecast(.interactive = TRUE)

# Paso 13: Seleccionar el mejor modelo (basado en RMSE)
best_model <- calibration_table %>%
  modeltime_accuracy() %>%
  arrange(rmse) %>%
  slice(1)

# Paso 14: Reajustar el mejor modelo a todos los datos
refit_table <- calibration_table %>%
  modeltime_refit(data = dta)

# Paso 15: Generar pronósticos para los próximos 20 meses
future_data <- dta %>%
  future_frame(
    .date_var   = date,
    .length_out = "24 months"
  )

forecast_24m <- refit_table %>%
  modeltime_forecast(
    new_data    = future_data,
    actual_data = dta
  )

# Paso 16: Visualizar el pronóstico final
forecast_24m %>%
  plot_modeltime_forecast(.interactive = TRUE)

0.3.2 Mejor modelo

# Extraer los valores predichos
predicciones <- forecast_24m %>%
  filter(.model_id == best_model$.model_id) %>%
  select(-.model_id, -.model_desc, -.key) %>%
  rename(ventas_predichas = .value) %>% 
  print()

## # A tibble: 24 × 4
##    .index     ventas_predichas .conf_lo .conf_hi
##    <date>                <dbl>    <dbl>    <dbl>
##  1 2024-06-01           17571.   16462.   18680.
##  2 2024-07-01           15372.   14263.   16481.
##  3 2024-08-01           17140.   16031.   18249.
##  4 2024-09-01           15936.   14827.   17045.
##  5 2024-10-01           16651.   15542.   17760.
##  6 2024-11-01           17349.   16240.   18458.
##  7 2024-12-01           17986.   16876.   19095.
##  8 2025-01-01           12718.   11609.   13827.
##  9 2025-02-01           14046.   12937.   15155.
## 10 2025-03-01           15979.   14870.   17089.
## # ℹ 14 more rows

#Visualizar predicciones
ggplot(predicciones, aes(x = .index)) +
  geom_line(aes(y = ventas_predichas, color = "Ventas Predichas")) +
  geom_ribbon(aes(ymin = .conf_lo, ymax = .conf_hi, fill = "Intervalo de Confianza"), alpha = 0.2) +
  geom_point(aes(y = ventas_predichas), color = '#1A5276') +
  labs(title = "Predicción de Ventas con Intervalos de Confianza",
       x = "",
       y = "Ventas",
       subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx") +
  scale_x_date(date_labels = "%Y-%m", date_breaks = "1 month") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom", 
        legend.title = element_blank()) +
  scale_y_continuous(labels = comma) +
  scale_color_manual(values = c("Ventas Predichas" = '#1A5276')) +
  scale_fill_manual(values = c("Intervalo de Confianza" = '#1A5276'))

0.4 C) Mejor Modelo

best_model

#GARCH
comparativa <- predicted_df %>% select(time, predicted) %>%
  #Tidymodels
  cbind(select(.data = predicciones, ventas_predichas))

# Graficar las líneas usando ggplot2
ggplot(comparativa, aes(x = time)) +
  geom_line(aes(y = ventas_predichas, color = "ARIMA(0,1,1)(2,1,1)[12] W/ XGBOOST ERRORS")) +
  geom_line(aes(y = predicted, color = "GARCH"), size = 1) +
  labs(title = "Comparación de Predicciones",
       x = "",
       y = "Ventas",
       subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx") +
  scale_x_date(date_labels = "%Y-%m", date_breaks = "1 month") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "bottom",
        legend.title = element_blank()) +
  scale_y_continuous(labels = comma) +
  scale_color_manual(values = c("ARIMA(0,1,1)(2,1,1)[12] W/ XGBOOST ERRORS" = '#1A5276', "GARCH" = '#FF5733'))

0.4.0.1 Serie completa

# Filtrar los valores originales
actual_values <- forecast_24m %>%
  filter(.model_desc == "ACTUAL") %>% 
  select(.index, .value)

# Filtrar los valores predichos por el best_model
# Suponemos que best_model_id es el ID del modelo seleccionado como mejor
predicted_values <- forecast_24m %>%
  filter(.model_id == best_model$.model_id) %>% 
  select(.index, .value)

# Unir ambas tablas por la columna de fecha (suponiendo que se llama 'date')
serie_completa <- rbind(actual_values, predicted_values)

0.5 D) Analitica del Negocio

¿Qué meses tienen el mayor número de ventas?

0.5.1 Descomponer serie

# Descomponer la serie temporal usando STL
last_date <- as.Date(max(dta$date, na.rm = TRUE))

drinks.ts <- ts(serie_completa$.value, start =c(1992,1), end = c(year(last_date),month(last_date)), frequency = 12)

decomposed <- stl(drinks.ts, s.window = "periodic")
seasonal_component <- decomposed$time.series[, "seasonal"]
trend_component <- decomposed$time.series[, "trend"]
remainder_component <- decomposed$time.series[, "remainder"]

# Convertir las componentes a un tibble
decomposed_tbl <- tibble(
  fecha = as.Date(time(drinks.ts)),
  original = drinks.ts,
  seasonal = seasonal_component,
  trend = trend_component,
  remainder = remainder_component
)

0.5.2 Estacionalidad de la serie

# Graficar la componente estacional
decomposed_tbl$fecha <- ymd(decomposed_tbl$fecha)

ggplot(decomposed_tbl, aes(x = as.Date(fecha, format = "%Y-%m-%d"))) +
   geom_line(aes(y = seasonal), color = "#581845") +
   labs(title = "Componente Estacional de la Serie Temporal",
        x = " ",
        y = "Estacionalidad",
        subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx") +
   theme_minimal() +
   scale_x_date(date_breaks = "1 year", date_labels = "%Y")+
   theme(axis.text.x = element_text(angle = 45, hjust = 1))

0.5.3 Meses de mayor venta

# Agregar una columna con el mes
sales_data <- serie_completa %>%
  mutate(month = format(.index, "%m"))

# Calcular el promedio de ventas por mes
average_sales_by_month <- sales_data %>%
  group_by(month) %>%
  summarise(average_sales = mean(.value))

# Ordenar por promedio de ventas
average_sales_by_month <- average_sales_by_month %>%
  arrange(desc(average_sales))

# Graficar las ventas promedio por mes
ggplot(average_sales_by_month, aes(x = month, y = average_sales)) +
  geom_bar(stat = "identity", fill = "#581845") +
  labs(title = "Promedio de Ventas por Mes",
       x = "Mes",
       y = "Ventas Promedio",
       subtitle = "Cifras en Millones de Dólares",
       caption = "Creado por: @EspacioColectivoMx")+
  theme_minimal()+
  scale_y_continuous(labels = comma)

0.5.3.1 Predicción de meses con mayor ventas

# Mostrar los meses con mayores ventas
top_months <- head(average_sales_by_month, 3)
top_months

0.5.3.2 Ventas Futuras para los principales meses

# Filtrar las predicciones para los meses con mayores ventas
predicted_sales <- predicciones %>%
  mutate(month = format(.index, "%m")) %>%
  filter(month %in% top_months$month)

predicted_sales

# Identificar los cuatro meses de mayor venta de cada año
top_ventas <- serie_completa %>% 
  mutate(year = lubridate::year(.index)) %>%
  group_by(year) %>%
  mutate(rank = rank(-.value)) %>%
  ungroup() %>%
  mutate(highlight = ifelse(rank <= 3, "Top 3", "Other"))

# Crear el gráfico
final.graf <- ggplot(decomposed_tbl, aes(x = as.Date(fecha))) +
  geom_line(aes(y = original, color = "Original"), size = .5) +
  geom_line(aes(y = seasonal, color = "Estacional"), size = .6) +
  geom_line(aes(y = trend, color = "Tendencia"), size = .6) +
  geom_line(aes(y = remainder, color = "Residuo"), size = .6) +
  geom_line(data = predicted_df, aes(x = time, y = predicted, color = "Predicción"),size = .8) +
  geom_point(data = top_ventas %>% filter(highlight == "Top 3"), 
             aes(x = .index, y = .value, color = "Top 3"), size = .8) +
  labs(title = "Pronostico y Descomposición de la Serie Temporal",
       subtitle = "Modelo SARIMA-GARCH. Cifras en Millones de Dólares",
       x = " ",
       y = "Ventas",
       color = "Componentes",
       caption = "Creado por: @EspacioColectivoMx") +
  scale_color_manual(values = c("Original" = "#1A5276",
                                "Estacional" = "#581845",
                                "Tendencia" = "#D80845",
                                "Residuo" = "#D8D845",
                                "Top 3" = "#FF5733",
                                "Predicción" = "#D80845")) +
  theme_minimal() +
  theme(legend.position = "bottom",
        legend.title = element_text(hjust = 0.5),
        legend.box = "horizontal") +
  guides(color = guide_legend(title = NULL, title.position = "top", title.hjust = 0.5)) +
  scale_y_continuous(labels = comma) +
  scale_x_date(date_breaks = "2 years", date_labels = "%Y")

final.graf

# Crear un data frame con los datos originales
tabla_final <- dta  %>%
  select(date,value)%>%
  rename(fecha = date,
         ventas = value)

# Crear un data frame con las predicciones
forecast_data <- data.frame(fecha = prediction_dates,
                            ventas = combined_forecast)

# Combinar los datos originales y las predicciones
tabla_final <- rbind(tabla_final,forecast_data)

#Vector fechas originales + preducciones
fechas <- seq.Date(from = as.Date(min(dta$date, na.rm = TRUE)),
                   to = as.Date(max(prediction_dates, na.rm = TRUE)),
                   by = "month")

#Serie de tiempo
tabla_final <- ts(tabla_final$ventas, start = c(1992,01), frequency = 12)

# Descomponer la serie temporal usando STL
new_decomposed <- stl(tabla_final, s.window = "periodic")
new_seasonal_component <- new_decomposed$time.series[, "seasonal"]
new_trend_component <- new_decomposed$time.series[, "trend"]
new_remainder_component <- new_decomposed$time.series[, "remainder"]

# Convertir las componentes a un tibble
new_decomposed_tbl <- tibble(
  fecha = fechas,
  original = tabla_final,
  seasonal = new_seasonal_component,
  trend = new_trend_component,
  remainder = new_remainder_component
)

#Filtrar a partir de las predicciones
tabla_final <- new_decomposed_tbl %>% 
  filter(fecha > max(dta$date, na.rm = TRUE))


##_____________ Unir con valores previos
#Data original
data_previa <-  decomposed_tbl %>%
  mutate(fecha = seq.Date(from = min(dta$date, na.rm = TRUE),
                          to = max(dta$date, na.rm = TRUE), by = "month"))

#Unión 
tabla_final <- dplyr::bind_rows(data_previa, tabla_final)

tabla_final <- tabla_final %>%
  mutate(estimación = case_when(
    fecha < as.Date("2024-04-01") ~ "valor original de la serie",
    fecha >= as.Date("2024-05-01") ~ "valor estimado pronósticado",
    TRUE ~ NA_character_  # Para manejar posibles fechas que no cumplan ninguna condición
  ))

tabla_final

0.5.4 Tasa de Crecimiento Mensual

# Asegúrate de tener los datos de monthly_sales con la tasa de crecimiento mensual calculada
monthly_sales <- serie_completa %>%
  mutate(date = as.Date(.index)) %>%
  group_by(year = year(date), month = month(date)) %>%
  summarise(total_sales = sum(.value)) %>%
  ungroup() %>%
  arrange(year, month) %>%
  mutate(monthly_growth_rate = (total_sales / lag(total_sales) - 1),
         date = as.Date(paste(year, month, "01", sep = "-"))) # Crear una columna de fecha para el gráfico

# Calculamos el promedio de la tasa de crecimiento mensual
mean_mensual_growth_rate <- mean(monthly_sales$monthly_growth_rate, na.rm = TRUE)

# Graficar la tasa de crecimiento mensual
ggplot(monthly_sales, aes(x = date, y = monthly_growth_rate)) +
  geom_line(aes(color = "Tasa de Crecimiento")) +
  geom_point(color = "#581845") +
  geom_hline(aes(yintercept = mean_mensual_growth_rate, color = "Promedio"), linetype = "dashed") +
  scale_color_manual(name = " ",
                     values = c("Tasa de Crecimiento" = "#581845",
                                "Promedio" = "#1A5276")) +
  labs(title = "Tasa de Crecimiento Mensual de Ventas",
       x = "Fecha",
       y = "Tasa de Crecimiento (%)",
       caption = "Creado por: @EspacioColectivoMx") +
  theme_minimal() +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                     breaks = seq(-0.5, 0.5, by = 0.1))+theme(
    axis.text.x = element_text(angle = 45, vjust = 1),
    legend.position = "bottom"
  )

0.5.5 Tasa de Crecimiento Anual

# Agrupar por año y calcular las ventas totales por año
annual_sales <- serie_completa %>%
  group_by(year = year(.index)) %>%
  summarise(total_sales = sum(.value)) %>%
  ungroup() %>%
  arrange(year) %>%
  mutate(annual_growth_rate = (total_sales / lag(total_sales) - 1))

# Calculamos el promedio de la tasa de crecimiento anual, ignorando NA.
mean_growth_rate <- mean(annual_sales$annual_growth_rate, na.rm = TRUE)

# Graficar la tasa de crecimiento anual
ggplot(annual_sales, aes(x = year, y = annual_growth_rate)) +
  geom_line(aes(color = "Tasa de Crecimiento")) +
  geom_point(color = "#581845") +
  geom_hline(aes(yintercept = mean_growth_rate, color = "Promedio"), linetype = "dashed") +
  scale_color_manual(name = " ", 
                     values = c("Tasa de Crecimiento" = "#581845", "Promedio" = "#1A5276")) +
  labs(title = "Tasa de Crecimiento Anual de Ventas",
       x = " ",
       y = "Tasa de Crecimiento (%)",
       caption = "Creado por: @EspacioColectivoMx") +
  theme_minimal() +
  scale_x_continuous(breaks = annual_sales$year) +
  scale_y_continuous(limits = c(-0.02, 0.1), 
                     breaks = seq(-0.02, 0.1, by = 0.01), 
                     labels = scales::percent_format(accuracy = 1), 
                     expand = c(0, 0)) +
  theme(
    axis.text.x = element_text(angle = 45, vjust = 1),
    legend.position = "bottom"
  )

annual_sales %>%
  ggplot(aes(x = year)) +
  # Crear barras para las ventas totales
  geom_col(aes(y = total_sales, fill = "Total Sales"), color = "black", width = 0.7, na.rm = TRUE) +
  
  # Añadir línea para la tasa de crecimiento anual
  geom_line(aes(y = annual_growth_rate * max(total_sales, na.rm = TRUE) / .50,
                color = "Growth Rate"), size = 1) +
  
  # Añadir puntos a la línea de tasa de crecimiento
  geom_point(aes(y = annual_growth_rate * max(total_sales, na.rm = TRUE) / .50, color = "Growth Rate"), size = 2) +
  
  # Configurar ejes Y: principal para ventas, secundario para tasa de crecimiento
  scale_y_continuous(
    name = "Ventas Totales",
    labels = comma,
    sec.axis = sec_axis(~ . / max(annual_sales$total_sales, na.rm = TRUE) * .50, name = "Tasa de Crecimiento Anual (%)", labels = scales::percent_format(accuracy = 1))
  ) +
  
  # Configurar eje X para mostrar todos los años
  scale_x_continuous(breaks = seq(min(annual_sales$year), max(annual_sales$year), by = 1)) +
  
  # Añadir títulos y etiquetas
  labs(title = "Ventas Totales y Tasa de Crecimiento ",
       subtitle = "Cifras Anuales en millones de dólares",
       x = " ",
       y = "Ventas Totales",
       caption = "Creado por: @EspacioColectivoMx") +

  # Definir colores para las barras y la línea
  scale_fill_manual(name = "", values = c("Total Sales" = "#1A5276")) +
  scale_color_manual(name = "", values = c("Growth Rate" = "#D8D845")) +
  
  # Aplicar tema minimalista
  theme_minimal() +
  
  # Personalizar la apariencia del gráfico
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  
  # Ajustar la leyenda
  guides(fill = guide_legend(override.aes = list(size = 1.5)))

# Tasa de crecimiento mensual anualizado
monthly_anualizado <- monthly_sales %>%
  filter(date > max(dta$date, na.rm = TRUE)) %>% 
  arrange(date) %>% 
  select(date,total_sales, monthly_growth_rate) %>% 
  # Calcular la tasa de crecimiento anualizada
  mutate(annualized_rate = round((1 + monthly_growth_rate)^12 - 1, 4))

0.5.6 Tasa de Crecimiento Trimestral

#fig.height=5, fig.width=11
# Agrupar por año y trimestre y calcular las ventas totales por trimestre
quarterly_sales <- serie_completa  %>%
  mutate(quarter = paste0(year(.index), " Q", quarter(.index))) %>%
  group_by(quarter) %>%
  summarise(total_sales = sum(.value)) %>%
  ungroup() %>%
  arrange(quarter) %>% 
  mutate(quarterly_growth_rate = (total_sales / lag(total_sales) - 1))

# Crear un vector de años únicos
years_unique <- unique(quarterly_sales$year)

# Promedio de crecimiento
mean_trim_growth_rate <- mean(quarterly_sales$quarterly_growth_rate, na.rm = TRUE)

# Graficar la tasa de crecimiento trimestral
ggplot(quarterly_sales, aes(x = quarter, y = quarterly_growth_rate, group = 1)) +
  geom_line(aes(color = "Tasa de Crecimiento")) +
  geom_point(color = "#581845") +
  geom_hline(aes(yintercept = mean_trim_growth_rate, color = "Promedio"), linetype = "dashed")+
  scale_color_manual(name = " ",
                     values = c("Tasa de Crecimiento" = "#581845",
                                "Promedio" = "#1A5276"))+
  labs(title = "Tasa de Crecimiento Trimestral de Ventas",
       x = " ",
       y = "Tasa de Crecimiento (%)",
       caption = "Creado por: @EspacioColectivoMx") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_discrete(labels = function(x) ifelse(substr(x, 6, 7) == "Q1", substr(x, 1, 4), ""))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1),
        legend.position = "bottom")+
  scale_y_continuous(limits = c(-0.3, 0.3), 
                     breaks = seq(-0.3, 0.3, by = 0.10), 
                     labels = scales::percent_format(accuracy = 1), 
                     expand = c(0, 0))

# Primero, creamos el dataframe con los datos trimestrales
quarterly_data <- serie_completa %>%
  group_by(year = year(.index), quarter = quarter(.index)) %>%
  summarise(total_sales = sum(.value), .groups = 'drop') %>%
  arrange(year, quarter) %>%
  mutate(date = as.Date(paste(year, quarter * 3 - 2, "01", sep = "-"))) %>%
  mutate(quarterly_growth_rate = (total_sales / lag(total_sales) - 1) * 100) %>%
  select(date, year, quarter, total_sales, quarterly_growth_rate)

ggplot(quarterly_data, aes(x = factor(year), y = total_sales, fill = as.factor(quarter))) +
  geom_col(position = "stack") +
  labs(title = "Ventas Totales por Trimestre",
       subtitle = "Cifras en millones de dólares",
       x = " ",
       y = "Ventas Totales",
       fill = "Trimestre",
       caption = "Creado por: @EspacioColectivoMx") +
  scale_y_continuous(labels = comma)+
  scale_fill_manual(values = c("#581845", "#D8D845", "#1A5276", "#D80845"), 
                    name = "", 
                    labels = c("Q1", "Q2", "Q3", "Q4")) +
  theme_minimal() +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, vjust = 1))

0.5.6.0.1 Guardar Excel

# Crear la carpeta "final" si no existe
dir.create("data/final", showWarnings = FALSE)

# Crear un nuevo libro de trabajo
wb <- createWorkbook()

# Añadir la primera hoja con la tabla
addWorksheet(wb, "serie_descomp")
writeData(wb, "serie_descomp", tabla_final)

# Guardar el gráfico como imagen
ggsave("img/final.graf.png", plot = final.graf, width = 10, height = 6, dpi = 300)

# Añadir una nueva hoja para el gráfico
addWorksheet(wb, "grafico_final")

# Insertar la imagen del gráfico en la hoja
insertImage(wb, "grafico_final", "img/final.graf.png", 
            width = 8, height = 5, units = "in",
            startRow = 2, startCol = 2)

# Guardar el libro de trabajo
saveWorkbook(wb, "data/final/serie_cerveza.xlsx", overwrite = TRUE)

# Eliminar el archivo temporal de la imagen
file.remove("data_final/grafico_temporal.png")

## [1] FALSE