TEMA2_Timeseries+GEFA - Ventas Hershey México (Pronóstico 6 meses)

EJERCICIO 2 “Hershey´s”

Instrucciones: Modela la sig. Serie de tiempo y elige la mejor opcion. Pronostica los siguientes 6 meses

Se cuenta con las ventas históricas mensuales de leche saborizada Hershey México (en miles de dólares), de enero 2017 a diciembre 2019 (36 observaciones). El objetivo es probar 5 modelos de pronóstico, compararlos y elegir el que mejor se ajuste a los datos históricos, para proyectar las siguientes 6 meses (enero–junio 2020).

INSTALACIÓN DE LIBRERÍAS

#install.packages("forecast")
#install.packages("readxl")
library(forecast)
library(readxl)

CARGA Y ACOMODO DE LOS DATOS

El archivo Ventas_Históricas_Lechitas.xlsx viene organizado en 3 bloques de columnas (uno por año: 2017, 2018 y 2019), cada uno con un par de columnas “Mes” / “Ventas”. Se elige el archivo, se leen los 3 bloques y se concatenan en una sola serie continua de 36 meses.

# En vez de file.choose() (que requiere selección manual e interactiva),
# se indica la ruta directa del archivo ya subido:
ruta_archivo <- "Ventas_Históricas_Lechitas.xlsx"

# Se brincan las primeras 2 filas (título y subtítulo "Miles de dólares")
# y se toma la fila 3 como encabezado (Mes, Ventas, Mes, Ventas, Mes, Ventas)
datos_crudos <- read_excel(ruta_archivo, sheet = "Hoja1", skip = 2)
colnames(datos_crudos)  # readxl renombra los duplicados como ...1, ...2, etc.

## [1] "Mes...1"    "Ventas...2" "Mes...3"    "Ventas...4" "Mes...5"   
## [6] "Ventas...6"

# Bloque 2017 = columnas 1 y 2 ; Bloque 2018 = columnas 3 y 4 ; Bloque 2019 = columnas 5 y 6
ventas_2017 <- na.omit(datos_crudos[[2]])
ventas_2018 <- na.omit(datos_crudos[[4]])
ventas_2019 <- na.omit(datos_crudos[[6]])

ventas <- c(ventas_2017, ventas_2018, ventas_2019)
mes <- 1:length(ventas)

dataframe2 <- data.frame(mes, ventas)
head(dataframe2, 12)

##    mes   ventas
## 1    1 25520.51
## 2    2 23740.11
## 3    3 26253.58
## 4    4 25868.43
## 5    5 27072.87
## 6    6 27150.50
## 7    7 27067.10
## 8    8 28145.25
## 9    9 27546.29
## 10  10 28400.37
## 11  11 27441.98
## 12  12 27852.47

nrow(dataframe2)  # deben ser 36 meses (12 x 3 años)

## [1] 36

# Serie de tiempo MENSUAL (frequency = 12), iniciando en Enero 2017
tsDOS <- ts(ventas, start = c(2017,1), frequency = 12)
tsDOS

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2017 25520.51 23740.11 26253.58 25868.43 27072.87 27150.50 27067.10 28145.25
## 2018 28463.69 26996.11 29768.20 29292.51 29950.68 30099.17 30851.26 32271.76
## 2019 32496.44 31287.28 33376.02 32949.77 34004.11 33757.89 32927.30 34324.12
##           Sep      Oct      Nov      Dec
## 2017 27546.29 28400.37 27441.98 27852.47
## 2018 31940.74 32995.93 32197.12 31984.82
## 2019 35151.28 36133.07 34799.91 34846.17

plot(tsDOS, main = "Ventas históricas de leche saborizada Hershey México",
     ylab = "Ventas (miles de dólares)", xlab = "Año", type = "o", col = "blue")

A diferencia del ejercicio anterior, aquí sí se observa una tendencia clara y sostenida: las ventas pasan de aproximadamente $23,700 a $36,100 miles de dólares entre 2017 y 2019, con ligeras caídas en enero/febrero de cada año (posible efecto estacional post-diciembre).

MODELO 1. NAIVE

naive2 <- naive(tsDOS, h = 6)
summary(naive2)

## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = tsDOS, h = 6) 
## 
## Residual sd: 1112.6435 
## 
## Error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 266.4474 1112.644 902.8514 0.8162353 3.013136 0.2579763 -0.5648044
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       34846.17 33420.26 36272.08 32665.43 37026.91
## Feb 2020       34846.17 32829.63 36862.71 31762.14 37930.20
## Mar 2020       34846.17 32376.42 37315.92 31069.02 38623.32
## Apr 2020       34846.17 31994.35 37697.99 30484.69 39207.65
## May 2020       34846.17 31657.74 38034.60 29969.88 39722.46
## Jun 2020       34846.17 31353.42 38338.92 29504.47 40187.87

plot(naive2)

MODELO 2. Promedio Móvil (k=3)

# Promedio móvil de orden 3 sobre los datos históricos (referencia / suavizado)
ma3 <- ma(tsDOS, order = 3)

# Pronóstico recursivo con promedio móvil simple k=3
ma_forecast <- function(serie, k = 3, h = 6){
  serie <- as.numeric(serie)
  pronosticos <- numeric(h)
  for (i in 1:h){
    ultimos <- tail(serie, k)
    siguiente <- mean(ultimos)
    pronosticos[i] <- siguiente
    serie <- c(serie, siguiente)
  }
  return(pronosticos)
}

pronostico_ma3 <- ma_forecast(tsDOS, k = 3, h = 6)
pronostico_ma3

## [1] 35259.72 34968.60 35024.83 35084.38 35025.94 35045.05

MODELO 3. Promedio Móvil Ponderado (k=3), Pesos: 3/6, 2/6, 1/6

pesos <- c(3/6, 2/6, 1/6)  # del más reciente al más antiguo

wma_forecast <- function(serie, pesos = c(3/6, 2/6, 1/6), h = 6){
  serie <- as.numeric(serie)
  k <- length(pesos)
  pronosticos <- numeric(h)
  for (i in 1:h){
    ultimos <- tail(serie, k)
    ultimos_ordenado <- rev(ultimos)
    siguiente <- sum(ultimos_ordenado * pesos)
    pronosticos[i] <- siguiente
    serie <- c(serie, siguiente)
  }
  return(pronosticos)
}

pronostico_wma3 <- wma_forecast(tsDOS, pesos = pesos, h = 6)
pronostico_wma3

## [1] 35045.23 34937.99 34958.44 34966.09 34958.85 34961.20

Tabla de datos (histórico + promedio móvil de referencia)

tabla_referencia <- data.frame(
  mes = mes,
  ventas = ventas,
  MA3 = as.numeric(ma3)
)
tabla_referencia

##    mes   ventas      MA3
## 1    1 25520.51       NA
## 2    2 23740.11 25171.40
## 3    3 26253.58 25287.37
## 4    4 25868.43 26398.29
## 5    5 27072.87 26697.26
## 6    6 27150.50 27096.82
## 7    7 27067.10 27454.28
## 8    8 28145.25 27586.21
## 9    9 27546.29 28030.64
## 10  10 28400.37 27796.21
## 11  11 27441.98 27898.27
## 12  12 27852.47 27919.38
## 13  13 28463.69 27770.76
## 14  14 26996.11 28409.33
## 15  15 29768.20 28685.60
## 16  16 29292.51 29670.46
## 17  17 29950.68 29780.79
## 18  18 30099.17 30300.37
## 19  19 30851.26 31074.06
## 20  20 32271.76 31687.92
## 21  21 31940.74 32402.81
## 22  22 32995.93 32377.93
## 23  23 32197.12 32392.62
## 24  24 31984.82 32226.13
## 25  25 32496.44 31922.84
## 26  26 31287.28 32386.58
## 27  27 33376.02 32537.69
## 28  28 32949.77 33443.30
## 29  29 34004.11 33570.59
## 30  30 33757.89 33563.10
## 31  31 32927.30 33669.77
## 32  32 34324.12 34134.23
## 33  33 35151.28 35202.82
## 34  34 36133.07 35361.42
## 35  35 34799.91 35259.72
## 36  36 34846.17       NA

MODELO 4. Suavizado Exponencial Simple, alpha = 0.2

ses2 <- ses(tsDOS, alpha = 0.2, h = 6, initial = "simple")
summary(ses2)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
## ses(y = tsDOS, h = 6, initial = "simple", alpha = 0.2)
## 
##   Smoothing parameters:
##     alpha = 0.2 
## 
##   Initial states:
##     l = 25520.5113 
## 
##   sigma:  1545.881
## Error measures:
##                    ME     RMSE      MAE     MPE    MAPE      MASE      ACF1
## Training set 1220.525 1545.881 1363.334 3.87354 4.44632 0.3895523 0.2043858
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       34308.29 32327.16 36289.42 31278.42 37338.16
## Feb 2020       34308.29 32287.93 36328.65 31218.42 37398.17
## Mar 2020       34308.29 32249.44 36367.14 31159.56 37457.03
## Apr 2020       34308.29 32211.66 36404.92 31101.78 37514.81
## May 2020       34308.29 32174.55 36442.03 31045.02 37571.56
## Jun 2020       34308.29 32138.08 36478.51 30989.23 37627.35

plot(ses2)

MODELO 5. ARIMA

arima2 <- auto.arima(tsDOS)
summary(arima2)

## Series: tsDOS 
## ARIMA(1,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ar1     sar1     drift
##       0.6383  -0.5517  288.8979
## s.e.  0.1551   0.2047   14.5026
## 
## sigma^2 = 202701:  log likelihood = -181.5
## AIC=371   AICc=373.11   BIC=375.72
## 
## Training set error measures:
##                    ME    RMSE    MAE        MPE      MAPE       MASE      ACF1
## Training set 25.22158 343.864 227.17 0.08059932 0.7069542 0.06491044 0.2081026

pronostico_arima <- forecast(arima2, level = c(95), h = 6)
pronostico_arima

##          Point Forecast    Lo 95    Hi 95
## Jan 2020       35498.90 34616.48 36381.32
## Feb 2020       34202.17 33155.28 35249.05
## Mar 2020       36703.01 35596.10 37809.92
## Apr 2020       36271.90 35141.44 37402.36
## May 2020       37121.98 35982.07 38261.90
## Jun 2020       37102.65 35958.90 38246.40

plot(pronostico_arima)

COMPARACIÓN DE MODELOS Y SELECCIÓN DEL MEJOR

acc_naive <- accuracy(naive2)
acc_ses   <- accuracy(ses2)
acc_arima <- accuracy(arima2)

valor_num <- as.numeric(tsDOS)
n <- length(valor_num)

# Ajuste in-sample MA(3)
ma_fitted <- sapply(4:n, function(t) mean(valor_num[(t-3):(t-1)]))
ma_actual <- valor_num[4:n]
ma_rmse <- sqrt(mean((ma_actual - ma_fitted)^2))
ma_mae  <- mean(abs(ma_actual - ma_fitted))

# Ajuste in-sample WMA(3)
wma_fitted <- sapply(4:n, function(t){
  ventana <- rev(valor_num[(t-3):(t-1)])
  sum(ventana * pesos)
})
wma_actual <- valor_num[4:n]
wma_rmse <- sqrt(mean((wma_actual - wma_fitted)^2))
wma_mae  <- mean(abs(wma_actual - wma_fitted))

tabla_comparacion <- data.frame(
  Modelo = c("Naive", "Promedio Movil (k=3)", "Promedio Movil Ponderado (k=3)",
             "Suavizado Exponencial (alpha=0.2)", "ARIMA"),
  RMSE = c(acc_naive[1,"RMSE"], ma_rmse, wma_rmse, acc_ses[1,"RMSE"], acc_arima[1,"RMSE"]),
  MAE  = c(acc_naive[1,"MAE"],  ma_mae,  wma_mae,  acc_ses[1,"MAE"],  acc_arima[1,"MAE"])
)

tabla_comparacion[order(tabla_comparacion$RMSE), ]

##                              Modelo      RMSE       MAE
## 5                             ARIMA  343.8640  227.1700
## 3    Promedio Movil Ponderado (k=3)  958.9819  814.1926
## 2              Promedio Movil (k=3) 1031.2592  870.4659
## 1                             Naive 1112.6435  902.8514
## 4 Suavizado Exponencial (alpha=0.2) 1545.8812 1363.3336

Como la serie tiene una tendencia creciente marcada, se espera que los modelos que solo promedian valores pasados sin adaptarse a la tendencia (Naive, Promedio Móvil, Promedio Móvil Ponderado, Suavizado Exponencial simple) se queden sistemáticamente POR DEBAJO de los valores reales, mientras que ARIMA —al poder incluir un término de diferenciación (d=1) y/o deriva (“drift”)— captura la tendencia y logra el menor RMSE/MAE de los cinco modelos.

PRONÓSTICO FINAL - MODELO GANADOR (6 MESES)

mejor_modelo <- tabla_comparacion[order(tabla_comparacion$RMSE), ][1, "Modelo"]
cat("Modelo elegido según RMSE/MAE mas bajos:", mejor_modelo, "\n")

## Modelo elegido según RMSE/MAE mas bajos: ARIMA

pronostico_final <- forecast(arima2, h = 6)

fechas_pronostico <- seq(as.Date("2020-01-01"), by = "month", length.out = 6)

tabla_pronostico_final <- data.frame(
  Mes = format(fechas_pronostico, "%Y-%m"),
  Pronostico_Ventas = round(as.numeric(pronostico_final$mean), 1),
  Limite_Inferior_95 = round(as.numeric(pronostico_final$lower[,2]), 1),
  Limite_Superior_95 = round(as.numeric(pronostico_final$upper[,2]), 1)
)
tabla_pronostico_final

##       Mes Pronostico_Ventas Limite_Inferior_95 Limite_Superior_95
## 1 2020-01           35498.9            34616.5            36381.3
## 2 2020-02           34202.2            33155.3            35249.0
## 3 2020-03           36703.0            35596.1            37809.9
## 4 2020-04           36271.9            35141.4            37402.4
## 5 2020-05           37122.0            35982.1            38261.9
## 6 2020-06           37102.6            35958.9            38246.4

plot(pronostico_final, main = "Pronóstico de ventas - siguientes 6 meses (Ene-Jun 2020)",
     ylab = "Ventas (miles de dólares)", xlab = "Año")

CONCLUSIÓN

Al comparar los cinco modelos, Naive, Promedio Móvil, Promedio Móvil Ponderado y el Suavizado Exponencial (alpha=0.2) tienden a subestimar los valores futuros, ya que estos métodos solo “miran hacia atrás” y no incorporan explícitamente una tendencia; su pronóstico converge hacia el promedio de las últimas observaciones, quedándose por debajo del nivel que la serie viene alcanzando.

El modelo ARIMA, seleccionado automáticamente con auto.arima(), es el que mejor se ajusta a los datos (menor RMSE y MAE), porque logra capturar la tendencia creciente de la serie mediante diferenciación y/o un término de deriva. Por esta razón, ARIMA se elige como el modelo final para pronosticar las ventas de los siguientes 6 meses (enero a junio de 2020), proyectando que la tendencia de crecimiento observada en 2017-2019 continuará, con ventas estimadas por encima de los $35,000 mil dólares mensuales y una banda de confianza al 95% que refleja la incertidumbre natural de proyectar hacia el futuro. Desde el punto de vista del negocio, esto sugiere que Hershey México debería planear capacidad de producción, inventario y distribución al alza para los primeros meses de 2020, dado que la demanda de este producto no muestra señales de desaceleración en el periodo analizado.