Ayudantía 2 Adm. de la Producción II Pronósticos con el método Holt-Winters
Eduardo Rojas y Yandary Rojas
Primer semestre 2023
Pronóstico mediante el método Holt-Winters
Para esta sección evaluaremos la precisión de los pronósticos, para lo cual primero realizaremos pronósticos basados en el método de Holt-Winters con unos parametros entregados de \(\alpha, \beta, \gamma\) para luego optimizarlos. Para esto utilizaremos el dataset Datos.xlsx que muestra la venta de Post-it durante el año \(2020-2023\) de la empresa MAS Inc.
info <- read_excel('Datos.xlsx')
head(info)Creación de la serie de tiempo
A continuación se muestra el código para inicializar la serie de tiempo que analizaremos.
dtfs <- ts(na.omit(info$Ventas), frequency = 12, start=c(2020,5)) #se debe señalizar la frecuencia
#ts(dataframe,frecuency=xxx,start=c(yyyy,m)) en x se debe colocar la frecuencia de los datos
#xxxx: mensual, trimestral, semestral, semanal, etc.
#en start se coloca la fecha inicial de los datos
info$Mes <- as.Date(info$Mes) #convierte la variable en tipo fecha
info$year <- year(info$Mes) #se agrega una nueva columna donde señalizar el año de cada fecha por separadoCabe mencionar que es muy importante que el dataframe que se entregue no contenga valores nulos o bien NA, para ello se usa el comando na.omit(df).
Por otra parte se debe señalizar como vienen los datos, en este caso trabajaremos de manera mensual, por ello se señaliza una frecuencia de 12.
Por último se debe señalizar el inicio de la serie temporal con la fecha correspondiente, donde, en este caso, se lanzó el producto del Post-it al mercado el mes de Mayo en el año 2020.
Graficar una serie de tiempo
A continuación se muestra cómo se puede graficar la serie de tiempo para analizar sus componentes tales como, Tendencia, Ciclo, Estacionalidad.
ggplot(info, aes(x = Mes, y = Ventas)) +
ggtitle ("Comportamiento de ventas de Post it \n 2020-2023") +
theme (plot.title = element_text(family="Comic Sans MS",
size=rel(2), #Tamaño relativo de la letra del título
vjust=2, #Justificación vertical, para separarlo del gráfico
face="bold", #Letra negrilla. Otras posibilidades "plain", "italic", "bold" y "bold.italic"
color="black", #Color del texto
lineheight=1.5)) + #Separación entre líneas
labs(x = "Fecha de venta",y = "Cantidad de Post it vendidos") + # Etiquetas o títulos de los ejes
geom_line() + #se agrega una linea para los datos
geom_point() + #se le agrega puntos a cada elemento del df
scale_x_date(date_labels = "%Y-%m-%d") + #formato en el que queremos la fecha
geom_vline(xintercept = as.Date("2023-01-01"), #realizar corte en una fecha especifica
linetype = 2, color = 2, linewidth = 1) +
stat_peaks(geom = "point", span = 15, color = "steelblue3", size = 2) + #destacar valores más altos
stat_peaks(geom = "label", span = 15, color = "steelblue3", angle = 0,
hjust = -0.1, x.label.fmt = "%Y-%m-%d") +
stat_peaks(geom = "rug", span = 15, color = "blue", sides = "b") + #destacar valores más bajos
stat_valleys(geom = "point", span = 11, color = "red", size = 2) +
stat_valleys(geom = "label", span = 11, color = "red", angle = 0,
hjust = -0.1, x.label.fmt = "%Y-%m-%d") +
stat_valleys(geom = "rug", span = 11, color = "red", sides = "b") 
Como se aprecia en el grafico se podría decir que presenta estacionalidad, por ello es viable utilizar el método Holt-Winters para pronosticar.
A continuación se realiza otro grafico para visualizar cómo se comportan las ventas de Post-it para cada año.
ggplot(info, aes(x = Mes, y = Ventas)) +
ggtitle ("Ventas de MAS Inc. por año") +
theme (plot.title = element_text(family="Comic Sans MS",
size=rel(2), #Tamaño relativo de la letra del título
vjust=2, #Justificación vertical, para separarlo del gráfico
face="bold", #Letra negrilla. Otras posibilidades "plain", "italic", "bold" y "bold.italic"
color="black", #Color del texto
lineheight=1.5)) + #Separación entre líneas
labs(x = "Fecha de venta",y = "Cantidad de Post it vendidos") + # Etiquetas o títulos de los ejes
geom_line() +
geom_point() +
facet_wrap(~year, scales = "free") + #separación gráfica por año (análisis detallado)
scale_x_date(date_labels = "%Y-%m-%d") + #formato en el que queremos la fecha
geom_vline(xintercept = as.Date("2023-01-01"), #destacar valores más bajos
linetype = 2, color = 2, linewidth = 1)
Como se aprecia en los graficos existe un comportamiento particularmente igual entre el año 2021 y 2022, lo que quiere decir nuevamente que podría existir estacionalidad en la serie de tiempo.
Pronóstico mediante el método Holt-Winters
A continuación se realiza un pronóstico mediante el método Holt-Winters considerando como parámetros \(\alpha=0.4,\beta=0.7,\gamma=0.45\).
HW1 <- HoltWinters(dtfs, alpha=0.4, beta=0.7,
gamma=0.45, seasonal = 'additive')#pronostico solicitado
plot(fitted(HW1)) #Gráfico que visualiza los elementos de una serie temporal
Como se aprecia en el gráfico se descompuso la serie de tiempo en sus diferentes componentes, donde no se visualiza una tendencia clara en el gráfico trend pero sí una estacionalidad debido a la similitud que existe entres dos periodos en el gráfico season.
Predicción de ventas para el año 2023
En esta sección se realizará la predicción con un nivel de confianza \(1-\alpha=0.95\) para los siguientes 9 meses del año 2023. También aprenderemos a cómo graficar la serie de tiempo junto a la predicción realizada mediante el método Holt-Winters.
HW1.pred <- predict(HW1, 9, prediction.interval = TRUE, level=0.95)#aqui se pronostica n periodos
Hw1.pred.df <- data.frame(HW1.pred) #convierte la información en un data frame
##Arreglo previo para visualizar pronóstico en un gráfico
FechaPrediccion <- as.Date(c("2023-04-01","2023-05-01","2023-06-01","2023-07-01",
"2023-08-01","2023-09-01","2023-10-01","2023-11-01",
"2023-12-01")) #fechas a pronosticar (Resto año 2023)
ValoresPrediccion <- Hw1.pred.df$fit#valores pronosticados con método Holt-Winters
PrediccionHW1 <- data.frame("Mes"=FechaPrediccion,
"Ventas"=ValoresPrediccion)#se crea data frame
PrediccionHW1$Tipo <- "Valores pronosticados" #se agrega una nueva columna que informa de cuales valores son los pronosticados
info$Tipo <- "Valores reales" #se agrega una nueva columna para indicar los valores originales
Temporalidad <- merge(x = info, y = PrediccionHW1, all = TRUE)#se juntan los valores reales y pronosticados en un solo dataframe
tail(Temporalidad,n=15)Se obtiene en un dataframe los valores originales y los valores que fueron pronosticados mediante el método Holt-Winters. Esto con la finalidad de poder graficar de manera detallada la serie de tiempo junto a su pronóstico para el resto del año 2023.
Graficando la predicción de ventas para el año 2023
ggplot(Temporalidad, aes(x = Mes, y = Ventas, color=Tipo)) +
ggtitle ("Pronóstico mediante el método Holt-Winters") +
theme (plot.title = element_text(family="Comic Sans MS",
size=rel(2), #Tamaño relativo de la letra del título
vjust=2, #Justificación vertical, para separarlo del gráfico
face="bold", #Letra negrilla. Otras posibilidades "plain", "italic", "bold" y "bold.italic"
color="black", #Color del texto
lineheight=1.5)) + #Separación entre líneas
labs(x = "Fecha de venta",y = "Cantidad de Post it vendidos") + # Etiquetas o títulos de los ejes
geom_line() + #se agrega una linea para los datos
geom_point() + #se le agrega puntos a cada elemento del df
scale_x_date(date_labels = "%Y-%m-%d") #formato en el que queremos la fecha
Como se aprecia en el grafico existe una similitud en el comportamiento estacional de la serie de tiempo en los valores pronosticados para el año 2023, pero esto no es suficiente para corroborar la exactitud de nuestro pronóstico. Para verificar que tan preciso fue el pronóstico se utilizarán las medidas de error tales como: MAPE, MAD, MSE, RMSE.
Calculando el error en el pronóstico
Para calcular las medidas de error, se creó una función denominada CalculandoErrores, la cual tiene como finalidad el calcular los errores de MAPE, MAD, MSE, RMSE.
CalculandoErrores <- function(pronostico,real) {
Comparando <- data.frame("Fecha"=real$Mes,
"Real"=real$Ventas,
"Pronostico"=pronostico$fitted)
Comparando= select(Comparando, Fecha, Real, Pronostico.xhat)
Comparando$Error <- abs(Comparando$Real-Comparando$Pronostico)
Comparando$"Error^2" <- Comparando$Error^2
Comparando$ErrorPorc <- Comparando$Error/Comparando$Real
MAD <-mean(Comparando$Error)
MAPE<-mean(Comparando$ErrorPorc)
MSE<-mean (Comparando$`Error^2`)
RMSE <- sqrt(MSE)
TipoErrores <- c("MAPE","MAD","MSE","RMSE")
Errores <- c(paste(round(MAPE*100,2),"%"),round(MAD,2),round(MSE,2),round(RMSE,2))
dfMedidasError <- data.frame(
"Error" = TipoErrores,
"Valor" = Errores
)
Lista <- list(dfMedidasError,Comparando)
return(Lista)
}Al usar esta función se nos será entregado una lista con dos dataframe en su interior, en la primera posición (Lista[1]) se encontrará la tabla con los errores calculados del pronóstico utilizando los parámetros iniciales. En el caso de la segunda posición (Lista[2]) se entrega una tabla resumen con los cálculos realizados.
A continuación se hará uso de la función anterior.
infoActualizada <- subset(info, Mes >= "2021-05-01")#se actualiza el dataframe con la fecha del 1er dato del fit realizado
ReporteErroresHW1<-CalculandoErrores(HW1,infoActualizada)
dfErroresHW1 <- data.frame(ReporteErroresHW1[1])
dfErroresHW1Se aprecia una tabla con los errores calculados con los parámetros iniciales entregados. A continuación se entregará la información detallada con los valores correspondientes de \(\alpha, \beta, \gamma\).
ErroresParametros <- data.frame("alpha"=c(0.4),"beta"=c(0.7),"gamma"=c(0.45),
"MAPE"=c(dfErroresHW1[1,2]),
"MAD"=c(dfErroresHW1[2,2]),
"MSE"=c(dfErroresHW1[3,2]),
"RMSE"=c(dfErroresHW1[4,2]))
ErroresParametrosOptimización de parametros
Se realizó una función que fuese capaz de probar las posibles combinaciones entre el parámetro \(\alpha, \beta, \gamma\) obteniendo el error de cada combinación en el pronóstico con el método Holt-Winters.
OptimizandoParametros <- function(ts,real,alcance,avance) {
dfOptimo <- data.frame(alpha = numeric(),
beta = numeric(),
gamma = numeric(),
MAPE = numeric(),
MAD = numeric(),
MSE=numeric(),
RSE=numeric())
lista <- list()
i <- 1
while(alcance < 1) {
lista[[i]] <- alcance
alcance <- alcance+avance
i<-i+1
}
lista[i] <- 1
for (a in lista){
for (b in lista){
for (c in lista){
HW <- HoltWinters(ts, alpha=a, beta=b,
gamma=c, seasonal = 'additive')#pronostico solicitado
AnalisisHW<-CalculandoErrores(HW,real)
dfErroresHW <- AnalisisHW[1]#se guarda como lista
dfErroresHW=data.frame(dfErroresHW)
ErroresHW <- dfErroresHW[1:4,2]
traspaso <- data.frame(alpha = a,beta = b,
gamma = c,
MAPE = ErroresHW[1],
MAD = as.numeric(ErroresHW[2]),
MSE= as.numeric(ErroresHW[3]),
RMSE=as.numeric(ErroresHW[4]))
dfOptimo <-rbind(dfOptimo,traspaso)
}
}
}
return(dfOptimo)
}Al final de la función se obtendrá una tabla con todas las combinaciones de parámetros según el alcance y el avance deseado, con ello se podrá encontrar la combinación que sea capaz de minimizar uno de los errores.
dfErrores<-OptimizandoParametros(dtfs,infoActualizada,0.1,0.05)
MinimoMSE <- min(dfErrores$MSE)
OptimosParametros <- filter(dfErrores,MSE==MinimoMSE)
OptimosParametrosSegún lo obtenido la combinación óptima que minimiza el MSE a 66,91 es: \(\alpha=0.1,\beta=0.35,\gamma=1\).