Series de Tiempo

Instalar paquetes y llamar librerías

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

Ejercicio. Hershey’s

INSTRUCCIONES: Modela la sigueinte Serie de Tiempo y elige la mejor opción. Pronostica las sigueintes 6 semanas.

ruta <- "/Users/adriana/Downloads/Ventas_Históricas_Lechitas.xlsx"
datos <- read_excel(ruta, sheet = 1, skip = 2)

## New names:
## • `Mes` -> `Mes...1`
## • `Ventas` -> `Ventas...2`
## • `Mes` -> `Mes...3`
## • `Ventas` -> `Ventas...4`
## • `Mes` -> `Mes...5`
## • `Ventas` -> `Ventas...6`

ventas <- as.numeric(c(datos[[2]], datos[[4]], datos[[6]]))  # concatena 2017-2018-2019
mes <- 1:length(ventas)
df2 <- data.frame(mes, ventas)
ts2 <- ts(ventas, start = c(2017,1), frequency = 12)
plot(ts2, main = "Ventas históricas de leche saborizada Hershey (Lechitas)",
     ylab = "Ventas (miles de dólares)", xlab = "Año")

# Modelo 1. Naive
naive2 <- naive(ts2, h = 6)
summary(naive2)

## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = ts2, 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

mape_naive2 <- accuracy(naive2)[1,"MAPE"]
# Modelo con menor MAPE es el mas adecuado.

# Función genérica de Promedio Móvil: promedia los últimos k valores para
# obtener el ajuste dentro de la muestra, y de forma recursiva usa los
# pronósticos ya calculados para proyectar h periodos hacia adelante.
ma_forecast <- function(serie, k, h){
  n <- length(serie)
  fitted <- rep(NA, n)
  for (i in (k+1):n) fitted[i] <- mean(serie[(i-k):(i-1)])
  mape <- mean(abs((serie[(k+1):n] - fitted[(k+1):n]) / serie[(k+1):n])) * 100
  ext <- c(serie, rep(NA, h))
  for (i in (n+1):(n+h)) ext[i] <- mean(ext[(i-k):(i-1)])
  list(fitted = fitted, mape = mape, forecast = ext[(n+1):(n+h)])
}

# Función genérica de Promedio Móvil Ponderado: los pesos se aplican del más
# antiguo al más reciente, por lo que el valor más reciente pesa más (3/6).
wma_forecast <- function(serie, pesos, h){
  k <- length(pesos)
  n <- length(serie)
  fitted <- rep(NA, n)
  for (i in (k+1):n) fitted[i] <- sum(serie[(i-k):(i-1)] * pesos)
  mape <- mean(abs((serie[(k+1):n] - fitted[(k+1):n]) / serie[(k+1):n])) * 100
  ext <- c(serie, rep(NA, h))
  for (i in (n+1):(n+h)) ext[i] <- sum(ext[(i-k):(i-1)] * pesos)
  list(fitted = fitted, mape = mape, forecast = ext[(n+1):(n+h)])
}

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

# Modelo 2. Promedio Móvil (k = 3)
ma2 <- ma_forecast(ventas, k = 3, h = 6)
ma2$mape

## [1] 2.811654

ma2$forecast

## [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, y 1/6)
wma2 <- wma_forecast(ventas, pesos, h = 6)
wma2$mape

## [1] 2.624516

wma2$forecast

## [1] 35045.23 34937.99 34958.44 34966.09 34958.85 34961.20

# Modelo 4. Suavizador Exponencial (alfa = 0.2). Menor alfa = mejor modelo
ses2 <- ses(ts2, alpha = 0.2, h = 6)
summary(ses2)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
## ses(y = ts2, h = 6, alpha = 0.2)
## 
##   Smoothing parameters:
##     alpha = 0.2 
## 
##   Initial states:
##     l = 26303.9678 
## 
##   sigma:  1574.86
## 
##      AIC     AICc      BIC 
## 661.0074 661.3710 664.1744 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 1111.747 1530.489 1335.458 3.458378 4.357771 0.3815873 0.3247295
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       34308.55 32290.28 36326.81 31221.88 37395.22
## Feb 2020       34308.55 32250.31 36366.78 31160.75 37456.34
## Mar 2020       34308.55 32211.10 36405.99 31100.78 37516.31
## Apr 2020       34308.55 32172.61 36444.48 31041.92 37575.17
## May 2020       34308.55 32134.81 36482.28 30984.10 37632.99
## Jun 2020       34308.55 32097.65 36519.44 30927.27 37689.82

mape_ses2 <- accuracy(ses2)[1,"MAPE"]

# Modelo 5. ARIMA
arima2 <- auto.arima(ts2)
arima2

## Series: ts2 
## 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

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

mape_arima2 <- accuracy(arima2)[1,"MAPE"]

# Tabla con el Resumen de los Resultados.
resumen2 <- data.frame(
  Modelo = c("Naive", "Promedio Móvil (k=3)", "Promedio Móvil Ponderado",
             "Suavizador Exponencial (alfa=0.2)", "ARIMA"),
  MAPE = round(c(mape_naive2, ma2$mape, wma2$mape, mape_ses2, mape_arima2), 2)
)
resumen2 <- resumen2[order(resumen2$MAPE), ]
resumen2

##                              Modelo MAPE
## 5                             ARIMA 0.71
## 3          Promedio Móvil Ponderado 2.62
## 2              Promedio Móvil (k=3) 2.81
## 1                             Naive 3.01
## 4 Suavizador Exponencial (alfa=0.2) 4.36

pronostico_final2 <- data.frame(
  Mes = 37:42,
  Naive = as.numeric(naive2$mean),
  Prom_Movil = ma2$forecast,
  Prom_Movil_Pond = wma2$forecast,
  Suavizador_Exp = as.numeric(ses2$mean),
  ARIMA = as.numeric(pronostico2$mean)
)
pronostico_final2

##   Mes    Naive Prom_Movil Prom_Movil_Pond Suavizador_Exp    ARIMA
## 1  37 34846.17   35259.72        35045.23       34308.55 35498.90
## 2  38 34846.17   34968.60        34937.99       34308.55 34202.17
## 3  39 34846.17   35024.83        34958.44       34308.55 36703.01
## 4  40 34846.17   35084.38        34966.09       34308.55 36271.90
## 5  41 34846.17   35025.94        34958.85       34308.55 37121.98
## 6  42 34846.17   35045.05        34961.20       34308.55 37102.65

# Conclusión Método ARE.
mejor2 <- resumen2$Modelo[1]
cat("El modelo con el MAPE más bajo (", round(resumen2$MAPE[1],2),
    "%) es:", mejor2,
    ". Por lo tanto es el modelo más adecuado para pronosticar las ventas de los próximos 6 meses de Lechitas.\n")

## El modelo con el MAPE más bajo ( 0.71 %) es: ARIMA . Por lo tanto es el modelo más adecuado para pronosticar las ventas de los próximos 6 meses de Lechitas.