Modelos de suavizamiento
Pablo Molineros
1/5/2021
Medias móviles
Los métodos de suavizado son una familia de métodos de pronóstico que promedian los valores durante múltiplos períodos para reducir el ruido y descubrir patrones en los datos.
Los promedios móviles son un enfoque de suavizado que promedian los valores de la serie en una ventana de períodos de tiempo consecutivos, generando una serie de promedios. Los enfoques difieren principalmente en función al número de valores promediados. Entre las principales técnicas tenemos:
- Media móvil centrada: técnica para calcular y mostrar una media móvil simple.
- Media móvil de medias móviles: uso del concepto de medias móviles simples para realizar un suavizado de varios pasos.
- Media móvil ponderada: suavizado mediante el uso de pesos elegidos específicamente por sus porpiedades matemáticas.
1. Media móvil centrada
Para este método, elegimos un número de puntos cercanos y promediamos para estimar la tendencia. Al calcular una media móvil simple, es beneficioso utilizar un número impar de puntos para que el cálculo sea simétrico. Por ejemplo, para calcular un promedio móvil con un horizante de 5, la fórmula sería:
\[ \hat{y}_t=\dfrac{y_{t-2}+y_{t-1}+y_t+y_{t+1}+y_{t+2}}{5} \]
donde \(t\) es el tiempo en el instante en el que se está suavizando y 5 el número de puntos que se utilizan para calcular el promedio.
Ejemplo 1.
library(tidyverse) #Manipulación y visualización de datos
library(lubridate) #Manejo de fechas
library(vars) #Métodos de suavizamiento
library(forecast) #Modelos predictivos
library(ggExtra) #Manipulación de gráficas
library(knitr) #Visualización de tablas#Carga de la serie de datos
load("serie_1.RData")
#Arreglo de fechas para posterior gráfica
r1_t <- r1_t %>% mutate(id_date=as.Date(paste0(id_date,"-01")))
#Visualización de los 10 primeros registros
head(r1_t,10)## # A tibble: 10 x 2
## id_date total
## <date> <dbl>
## 1 2018-01-01 2
## 2 2018-02-01 8
## 3 2018-03-01 5
## 4 2018-04-01 11
## 5 2018-05-01 8
## 6 2018-06-01 18
## 7 2018-07-01 6
## 8 2018-08-01 7
## 9 2018-09-01 8
## 10 2018-10-01 12#Gráfico de la serie
p1_t <- ggplot(r1_t, aes(id_date, total)) +
geom_line() + #Gráfico de lineas
geom_point(size=.6, colour="#2072B2") + #coloca los puntos
theme_bw() + #tema del gráfico en blanco
theme(axis.title.x = element_blank(), #remuevo el titulo del eje x
axis.text.x = element_text(size = 5, colour = "black", angle = 60, vjust = .5), #texto en eje x
axis.ticks.y = element_blank(), #remueve los marcadores del eje y
axis.title.y = element_text(size = 6), #tamaño titulo en y
axis.text.y = element_text(size = 5, colour = "black"), #tamaño del texto en eje y
plot.title = element_text(size=8, hjust = .5), #tamaño del titulo centrado
panel.grid.major = element_line(colour = "#2072B2", #color, tamaño y diseño de
size = .05, #la malla de la gráfica
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", limits = c(0,30), breaks = seq(0,30,5))+ #escalo el eje y
ggExtra::removeGridX() + #remueve grid del eje x
ggtitle("Serie repuesto AC/DCDV1") + #titulo de la gráfica
scale_x_date(date_labels = "%b%Y", #etiquetas en las fechas, limites y saltos
limits = as.Date(c('2018-01-01','2020-01-01')),
date_breaks = "2 month") +
geom_hline(yintercept = mean(r1_t$total), color = "blue", linetype = "dashed", size=.3)
p1_t 
#Transformación en serie de tiempo para calculos
r1_ts <- ts(r1_t[,2], frequency = 12, start = c(2018,1))
#Descomposición de la serie
autoplot(decompose(r1_ts)) + theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black"),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed"))+
ggExtra::removeGridX() +
ggtitle("Descomposición de la serie repuesto AC/DCDV1") 
La serie de tiempo que se estudia en este ejemplo parece tener una ligera tendencia y estacionalidad a través del tiempo. Se emplea una media móvil simple como indicador de la tendencia para esta serie.
#Definimos el horizonte para el promedio
k <- 7
#Calculamos la media movil simple
sma <- rollmean(r1_ts, k, fill = NA)
sma## Jan Feb Mar Apr May Jun Jul
## 2018 NA NA NA 8.285714 9.000000 9.000000 10.000000
## 2019 15.000000 15.857143 15.571429 17.428571 16.714286 17.714286 18.428571
## 2020 13.428571 12.142857 11.857143 12.142857 12.285714 11.857143 12.571429
## Aug Sep Oct Nov Dec
## 2018 9.857143 10.285714 10.285714 12.285714 13.857143
## 2019 18.428571 17.000000 17.571429 16.000000 16.000000
## 2020 14.285714 13.714286 NA NA NA#Transformamos el suavizamiento a un data frame
sma <- data.frame(sma)
#Añadimos las fechas correspondientes
sma <- sma %>% mutate(id_date=as.Date(paste0(rep(c(2018, 2019, 2020), each=12),"-",rep(1:12, 3),"-",1)))
#Creamos una base con las dos series
p2_t <- merge(r1_t, sma, by="id_date")
p2_tMelted <- reshape2::melt(p2_t, id.var='id_date')#Graficamente
ggplot(p2_tMelted, aes(x=id_date, y=value, col=variable)) + geom_line() +
theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", angle = 60, vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
legend.title=element_blank(),
legend.text = element_text(size=5),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", limits = c(0,30), breaks = seq(0,30,5))+
ggExtra::removeGridX() +
ggtitle("Comparativa series AC/DCDV1 y SMA") +
scale_x_date(date_labels = "%b%Y",
date_breaks = "2 month") +
scale_color_discrete(labels=c("Serie AC/DCDV1", "SMA"))
1.1 Media móvil asimétrica (Predicción)
Los promedios móviles centrados se calculan promediando los datos tanto del pasado como en el futuro de un punto determinado de tiempo. En este sentido, no se pueden utilizar para realizar pronósticos porque, en el momento de la predicción, el futuro es incierto. Por lo tanto, se usan medias móviles asimétricas, donde la ventana de \(k\) períodos se coloca sobre los valores \(k\) disponibles más recientes de la serie. Por ejemplo, si tenemos datos hasta el período de tiempo \(t\), podemos predecir el valor de \(t+1\) promediando sobre \(k\) períodos antes de \(t+1\), la fórmula sería:
\[ \hat{y}_{t+1}=\dfrac{y_{t-4}+y_{t-3}+y_{t-2}+y_{t-1}+y_{t}}{5} \]
Ejemplo 1.1.
Continuando con el ejemplo propuesto, se aplica la media móvil asimétrica con una ventana de 12 días promedio para pronosticar la demanda en los siguientes períodos.
#Media móvil asimétrica horizonte de 12 días
sma_a <- rollmean(r1_ts, 12, align = "right", fill = NA)
#Transformamos a un data frame y calculamos la predicción para el 2020-12-01
sma_a <- data.frame(sma_a)
#Agregamos las fechas
sma_a <- sma_a %>% mutate(id_date=as.Date(paste0(rep(c(2018, 2019, 2020), each=12),"-",rep(1:12, 3),"-",1)))
#Creamos una base con las dos series
p3_t <- merge(r1_t, sma_a, by="id_date")
p3_tMelted <- reshape2::melt(p3_t, id.var='id_date')ggplot(p3_tMelted, aes(x=id_date, y=value, col=variable)) + geom_line() +
theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", angle = 60, vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
legend.title=element_blank(),
legend.text = element_text(size=5),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", limits = c(0,30), breaks = seq(0,30,5))+
ggExtra::removeGridX() +
ggtitle("Comparativa series AC/DCDV1 y SMA Asimétrico") +
scale_x_date(date_labels = "%b%Y",
date_breaks = "2 month") +
scale_color_discrete(labels=c("Serie AC/DCDV1", "SMA A."))
#Predicción para ene2021
ene21 <- sum(p3_t[25:36,2])/12; ene21## [1] 13.08333#Predicción para los próximos 6 meses
val <- p3_t$total.x
fit <- NA
k <- 1:6
for(i in seq_along(k)){
fit[i] <- sum(val[length(val):(length(val)-11)])/12
val <- c(val,fit[i])
}
val[37:42]## [1] 13.08333 13.17361 13.18808 13.20375 13.55406 13.51690
2. Media móvil de medias móviles
Esta técnica se emplea a menudo con un número par de puntos de datos para que el producto final sea simétrico alrededor de cada punto. Por ejemplo, un promedio de media móvil de un número par está desiquilibrado y, para nuestro propósito, el desequilibrio estará a favor de observaciones más recientes.
\[ \hat{y}_t=\dfrac{y_{t-1}+y_t+y_{t+1}+y_{t+2}}{4} \]
Para hacer que la media móvil sea simétrica (más precisa), tomamos un 2-MA del 4-MA para crear un 2x4-MA. Para el paso 2-MA, promediamos las medias móviles actuales y anteriores.
\[ \hat{y}_t=\frac{1}{8}y_{t-2}+\frac{1}{4}y_{t-1}+\frac{1}{4}y_{t}+\frac{1}{4}y_{t+1}+\frac{1}{8}y_{t+2} \]
3. Media móvil poderada
Un promedio móvil de un promedio móvil se puede considerar como un MA simétrico que tiene diferentes pesos en cada observación cercana. Por ejemplo, el 2x4-MA es equivalente a un 5-MA ponderado por \([\frac{1}{8}, \frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{8}]\). En general, un m-MA ponderado se puede escribir como:
\[ \hat{T}_t=\sum_{j=-k}^{k}a_{j}y_{t+j} \] donde \(k=\frac{(m-1)}{2}\) y los pesos vienen dados por \([a_{-k}, \dots, a_{k}]\). Es importante que todos los pesos sumen uno y que sean simétricos para que \(a_j=a_{-j}\). Este simple m-MA es un caso especial donde todos los pesos son iguales a \(\frac{1}{m}\). Una de las principales ventajas es que ofrecen una estimación más suave del ciclo de tendencia. En lugar de que las observaciones ingresen y salgan del cálculo con el peso total, sus pesos se aumentan lentamente y luego disminuyen lentamente, lo que da como resultado una curva más suave.
Suavizado exponencial simple
Este método es adecuado para datos sin tendencia o patron estacional. Se pondera las observaciones recientes más que las observaciones anteriores. El peso de cada ponderación se determina mediante el uso de un parámetro de suacvizado, que se denota como \(\alpha\). Para un conjunto con \(T\) observaciones, se calcula el valor predicho.
\[ \hat{y}_t=\alpha y_t+\alpha(1-\alpha)y_{t-1}+\cdots+\alpha(1-\alpha)^{t-1}y_t \]
donde \(0 < \alpha \leq 1\). También es común utilizar la forma componente de este modelo dada por: \[ \hat{y}_{t+1}=l_t \]
\[ l_t=\alpha y_t +(1-\alpha)l_{t-1} \]
En ambas ecuaciones podemos observar que el maypr peso se coloca a la observación más reciente. En la práctica, \(\alpha\) es igual a \(0.1 - 0.2\) tiende a funcionar bastante bien. Cuando \(\alpha\) está más cerca de cero, consideramos este como aprendizaje lento porque el algoritmo da más peso a los datos históricos. Cuando \(\alpha\) está más cerca de 1 consideramos este aprendizaje rápido porque el algoritmo da más peso a la observación más reciente; por lo tanto, los cambios recientes en los datos tendrá un mayor impacto en los valores pronosticados.
Continuando con el ejemplo, se aplica el método de suavizamiento para la serie de los repuestos.
#Muestra de entrenamiento
train.ts <- ts(r1_t[1:30,2], frequency = 12, start = c(2018,1))
#Muestra de validacion
test.ts <- ts(r1_t[30:36,2], frequency = 12, start = c(2020,8))
test.ts1 <- ts(r1_t[31:36,2], frequency = 12, start = c(2020,8))
#Suavizado exponencial simple
ses <- ses(train.ts, h = 4, alpha = .2)
ses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 13.59617 6.384798 20.80755 2.567326 24.62502
## Aug 2020 13.59617 6.241985 20.95036 2.348912 24.84344
## Sep 2020 13.59617 6.101892 21.09046 2.134659 25.05769
## Oct 2020 13.59617 5.964371 21.22798 1.924339 25.26801#Grafico de la predicción del modeloo ses
p4_t <- autoplot(ses) + theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
legend.title=element_blank(),
legend.text = element_text(size=5),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", limits = c(0,30), breaks = seq(0,30,5))+
ggExtra::removeGridX() +
ggtitle("Proyecciones de la serie AC/DCDV1 con SES")
p4_t
El método mejora en sus predicciones cuando la serie no presenta ninguna tendencia. Por lo tanto, procedemos a diferenciar la serie original de los datos
dr1_ts <- diff(train.ts)
dr1_ts.test <- diff(test.ts)#Gráfico de la serie diferenciada
p4_t <- autoplot(dr1_ts) +
theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", breaks = seq(-10,15,5))+
ggExtra::removeGridX() +
ggtitle("Serie diferenciada repuesto AC/DCDV1")
p4_t 
#Calculamos el SES de la serie diferenciada
ses_dif <- ses(dr1_ts, alpha=.2, h = 4)
ses_dif## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 0.1051907 -9.594044 9.804426 -14.72851 14.93889
## Aug 2020 0.1051907 -9.786127 9.996508 -15.02227 15.23266
## Sep 2020 0.1051907 -9.974550 10.184931 -15.31044 15.52082
## Oct 2020 0.1051907 -10.159515 10.369896 -15.59332 15.80370#Grafico de la predicción del modeloo ses
p5_t <- autoplot(ses_dif) + theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
legend.title=element_blank(),
legend.text = element_text(size=5),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n",breaks = seq(-10,15,5))+
ggExtra::removeGridX() +
ggtitle("Proyecciones de la serie diferenciada AC/DCDV1 con SES")
p5_t
Para comparar que tan bien predice el modelo, podemos comparar nuestros pronósticos con nuestro conjunto de datos de validación.
kable(accuracy(ses_dif,dr1_ts.test), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.2439326 | 7.302714 | 5.64375 | NaN | Inf | 0.7159981 | -0.6005006 | NA |
| Test set | -0.6051907 | 3.551937 | 3.50000 | 99.56171 | 99.56171 | 0.4440299 | -0.5000000 | 0.4135442 |
Método Holt
El método Holt realiza predicciones para datos con tendencia pero sin estacionalidad utilizando dos parámetros de suavizamiento, \(\alpha\) y \(\beta\), que corresponden a los componentes de nivel y tendencia, respectivamente. Para el método de Holt, la predicción será una línea de una pendiente distinta de cero que se extiende desde el paso del tiempo despúes del útlimo punto de datos recopilados en adelante.
El pronóstico de k-step-forward se obtiene combinando el nivel estimado en el tiempo \(t(L_t)\) y la estimación de tendencia en el tiempo \(t(T_t)\).
\[ \hat{y}_{T+1}=L_{T}+kT_{t} \]
El nivel \(L_t\) y tendencia \(T_t\) se actualizan a través de un par de ecuaciones de actualización, que es donde entran los parámetros de suavizamiento.
\[ L_t=\alpha y_t+ \alpha(1-\alpha)(L_{t-1}+T_{t-1}) \] \[ T_t=\beta(L_t-L_{t-1})+(1-\beta)T_{t-1} \] En estas ecuaciones, la primera significa que el nivel en el momento \(t\) es un promedio ponderado del valor real en el momento \(t\) y el nivel en el período anterior, ajustado por la tendencia. La segunda ecuación significa que la tendencia en el momento \(t\) es un promedio ponderado de la tendencia en el período anterior y la información más reciente sobre el cambio en el nivel. Similar al SES, \(\alpha\) y \(\beta\) se limitan a 0 - 1 con valores más altos proporcionan un aprendizaje más rápido y valores más bajos que proporcionan un aprendizaje más lento.
Para capturar una tendencia multiplicativa, se hace un pequeño ajuste en las ecuaciones anteriores: \[ \hat{y}_{t+1}=L_t \times kT_t \] \[ L_t = \alpha y_t + \alpha(1-\alpha) (L_{t-1}\times T_{t-1}) \] \[ T_t = \beta \left( \dfrac{L_t}{L_{t-1}} \right) +(1-\beta)T_{t-1} \]
Aplicando el método de Holt a la serie en estudio:
hlt <- holt(train.ts, h=6)
hlt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 14.13846 6.848548 21.42837 2.989501 25.28742
## Aug 2020 14.40264 6.816103 21.98918 2.800031 26.00526
## Sep 2020 14.66683 6.794632 22.53902 2.627345 26.70631
## Oct 2020 14.93101 6.782982 23.07903 2.469678 27.39234
## Nov 2020 15.19519 6.780187 23.61019 2.325554 28.06483
## Dec 2020 15.45937 6.785429 24.13332 2.193720 28.72502p6_t <- autoplot(hlt) +
theme_bw() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(size = 5, colour = "black", vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("Cantidad enviada\n", breaks = seq(-10,30,5))+
ggExtra::removeGridX() +
ggtitle("Serie repuesto AC/DCDV1 método de Holt")
p6_t 
hlt$model## Holt's method
##
## Call:
## holt(y = train.ts, h = 6)
##
## Smoothing parameters:
## alpha = 0.2881
## beta = 1e-04
##
## Initial states:
## l = 5.9636
## b = 0.2642
##
## sigma: 5.6883
##
## AIC AICc BIC
## 212.0481 214.5481 219.0541kable(accuracy(hlt,test.ts1), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.0053028 | 5.295565 | 4.139699 | -17.05496 | 40.56978 | 0.4543572 | -0.0233594 | NA |
| Test set | -0.5310076 | 4.124122 | 3.869183 | -12.83840 | 29.62489 | 0.4246664 | 0.1790143 | 0.6522865 |
Podemos intentar mejorar la predicción ajustando el modelo para diferentes valores de \(\beta\)
beta <- seq(.001,.5, .001)
RMSE <- NA
for(i in seq_along(beta)){
fit <- holt(train.ts, beta = beta[i], h = 6)
RMSE[i] <- accuracy(fit, test.ts1)[2,2]
}
beta.fit <- data.frame(beta, RMSE)
beta.min <- filter(beta.fit, RMSE==min(RMSE))
beta.min ## beta RMSE
## 1 0.22 3.693711ggplot(beta.fit, aes(beta, RMSE)) +
geom_line() +
geom_point(data=beta.min, aes(beta, RMSE), size=1, color="blue") +
theme_bw() +
theme(axis.title.x = element_text(size=6, colour = "black"),
axis.text.x = element_text(size = 5, colour = "black", vjust = .5),
axis.ticks.y = element_blank(),
axis.title.y = element_text(size = 6),
axis.text.y = element_text(size = 5, colour = "black"),
plot.title = element_text(size=8, hjust = .5),
panel.grid.major = element_line(colour = "#2072B2",
size = .05,
linetype = "dashed")) +
scale_y_continuous("RMSE", breaks = seq(0,16,2))+
scale_x_continuous("Beta") +
ggExtra::removeGridX() +
ggtitle("Valor óptimo del parámetro beta")
Probamos estimar nuevamente el modelo con el valor de \(\beta\) óptimo.
hlt_1 <- holt(train.ts, h=6, beta = .22)
hlt_1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 13.12214 4.636840 21.60743 0.1449953 26.09928
## Aug 2020 13.79062 3.431815 24.14942 -2.0518055 29.63304
## Sep 2020 14.45910 1.486902 27.43131 -5.3801665 34.29837
## Oct 2020 15.12759 -1.055409 31.31058 -9.6221702 39.87735
## Nov 2020 15.79607 -4.081267 35.67341 -14.6036942 46.19584
## Dec 2020 16.46456 -7.512626 40.44174 -20.2053790 53.13449hlt_1$model## Holt's method
##
## Call:
## holt(y = train.ts, h = 6, beta = 0.22)
##
## Smoothing parameters:
## alpha = 0.4802
## beta = 0.22
##
## Initial states:
## l = 2.8158
## b = 1.8178
##
## sigma: 6.6211
##
## AIC AICc BIC
## 219.1587 220.7587 224.7635kable(accuracy(hlt_1,test.ts1), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.1741343 | 6.163920 | 4.896282 | -18.08756 | 47.13131 | 0.5373968 | -0.1243392 | NA |
| Test set | -0.7275878 | 3.693711 | 3.423337 | -13.31291 | 26.84947 | 0.3757321 | 0.1218431 | 0.5968626 |
Método de Holt Winters
Para realizar predicciones utilizando datos con tendencia y estacionalidad, recurrimos al método estacional de Holt Winters. Este método se puede implementar con una estructura “aditiva” o “multiplicativa”, donde la elección del modelo depende del conjunto de datos. El modelo aditivo se utiliza mejor cuando la tendencia estacional es de la misma magnitud en el conjunto de datos, mientras que se prefiere el modelo multiplicativo cuando la magnitud de la estacionalidad cambia a medida que aumenta el tiempo.
Modelo Aditivo
Para el modelo aditivo, la forma de la ecuación regular es: \[ \hat{y}_{t+1}=L_t+kT_t+S_{t+k-m} \] El nivel \(L_t\), la tendencia \(T_t\) y la estacionalidad \(S_t\) se actualizan a través de un par de ecuaciones de actualización donde intervienen los parámetros de suavizamiento.
\[ L_t=\alpha(y_t-S_{t-m})+(1-\alpha)(L_{t-1}+ T_{t-1}) \]
\[ T_t=\beta(L_{t-1}+ T_{t-1})+(1-\beta)T_{t-1} \] \[ S_t=\gamma(y_t-L_t)+(1-\gamma)S_{t-m} \] De manera similar al SES y al método de Holt,los tres parámetros están restringidos al intervalo \((0,1]\).
Aplicando el método de Holt-Winters a la serie en estudio:
hw <- ets(train.ts, model = "AAA")
summary(hw)## ETS(A,A,A)
##
## Call:
## ets(y = train.ts, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.2162
## beta = 1e-04
## gamma = 0.6296
##
## Initial states:
## l = 5.5737
## b = 0.3134
## s = -2.4555 -4.1799 3.2755 -2.5057 -8.2977 -6.5536
## 10.6981 0.8942 2.4752 -0.5081 -2.4323 9.5898
##
## sigma: 10.3776
##
## AIC AICc BIC
## 253.5505 304.5505 277.3709
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.2222448 7.089231 5.203259 -26.05345 64.58367 0.5710894
## ACF1
## Training set -0.01327052hw_f <- forecast(hw, h=6)
hw_f## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 20.097107 6.797713 33.39650 -0.2425603 40.43677
## Aug 2020 10.892772 -2.714115 24.49966 -9.9171659 31.70271
## Sep 2020 20.253089 6.345228 34.16095 -1.0171498 41.52333
## Oct 2020 17.716596 3.513862 31.91933 -4.0046104 39.43780
## Nov 2020 11.761679 -2.730194 26.25355 -10.4017288 33.92509
## Dec 2020 7.359078 -7.416542 22.13470 -15.2382823 29.95644kable(accuracy(hw_f, test.ts1), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.2222448 | 7.089231 | 5.203259 | -26.053454 | 64.58367 | 0.5710894 | -0.0132705 | NA |
| Test set | 0.8033573 | 8.221212 | 7.148340 | -8.454511 | 48.30016 | 0.7845739 | 0.3160046 | 1.506997 |
checkresiduals(hw_f)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 15.458, df = 3, p-value = 0.001464
##
## Model df: 16. Total lags used: 19
Modelo Multiplicativo
Si tenemos estacionalidad multiplicativa, entonces nuestra forma de la ecuación cambia:
\[ \hat{y}_{t+1}= (L_{t} +kT_{t})S_{t+k-m} \] El nivel \(L_t\), la tendencia \(T_t\) y la estacionalidad \(S_t\) se actualizan a través de un par de ecuaciones de actualización donde intervienen los parámetros de suavizamiento.
\[ L_t=\alpha \dfrac{y_t}{S_{t-m}}+(1-\alpha)(L_{t-1}+T_{t-1}) \]
\[ T_t=\beta(L_{t}-L_{t-1})+(1-\beta)T_{t-1} \] \[ S_t = \gamma \dfrac{y_t}{L_t}+(1-\gamma)S_{t-m} \]
hw_m <- ets(train.ts, model = "MMM")
summary(hw_m)## ETS(M,M,M)
##
## Call:
## ets(y = train.ts, model = "MMM")
##
## Smoothing parameters:
## alpha = 0.1458
## beta = 0.1424
## gamma = 0.0029
##
## Initial states:
## l = 8.0605
## b = 1.0239
## s = 0.7042 0.8237 1.1328 1.1125 0.6453 0.9269
## 1.6749 1.0641 1.1795 0.8384 0.985 0.9126
##
## sigma: 0.4895
##
## AIC AICc BIC
## 223.5533 274.5533 247.3736
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.839297 8.45185 4.416584 -29.94618 48.01181 0.4847471 0.2102473hwm_f <- forecast(hw_m, h=6)
hwm_f## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 7.640882 2.8523952 12.488628 0.2845318 15.02008
## Aug 2020 4.822279 1.7984342 8.104849 0.2837716 10.01423
## Sep 2020 7.525233 2.7181165 13.229900 0.3739944 16.78250
## Oct 2020 6.944160 2.2358446 13.628637 0.2274005 18.74688
## Nov 2020 4.574542 1.3278182 9.918623 0.2261158 14.73759
## Dec 2020 3.540794 0.8629064 8.921983 0.1880326 14.30815kable(accuracy(hwm_f, test.ts1), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -1.839296 | 8.45185 | 4.416584 | -29.94618 | 48.01181 | 0.4847471 | 0.2102473 | NA |
| Test set | 8.918598 | 10.49940 | 8.918598 | 56.43483 | 56.43483 | 0.9788706 | 0.2941233 | 1.954509 |
checkresiduals(hwm_f)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,M,M)
## Q* = 41.631, df = 3, p-value = 4.805e-09
##
## Model df: 16. Total lags used: 19hw_m1 <- ets(train.ts, model = "MAM")
summary(hw_m1)## ETS(M,A,M)
##
## Call:
## ets(y = train.ts, model = "MAM")
##
## Smoothing parameters:
## alpha = 0.1388
## beta = 0.1388
## gamma = 1e-04
##
## Initial states:
## l = 8.6143
## b = 0.0336
## s = 0.6905 0.7984 1.1889 1.0968 0.6865 0.9147
## 1.8665 1.057 1.1208 0.7631 0.9099 0.907
##
## sigma: 0.5168
##
## AIC AICc BIC
## 223.6023 274.6023 247.4227
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.050813 8.343387 4.624479 -24.93227 49.10704 0.5075647 0.1036556hwm1_f <- forecast(hw_m1, h=6)
hwm1_f## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2020 6.9661575 2.352342 11.579973 -0.0900650 14.022380
## Aug 2020 4.2457772 1.232768 7.258787 -0.3622233 8.853778
## Sep 2020 5.2139557 0.666400 9.761511 -1.7409305 12.168842
## Oct 2020 3.9505050 -1.325990 9.227000 -4.1191984 12.020208
## Nov 2020 1.5105357 -2.649404 5.670475 -4.8515425 7.872614
## Dec 2020 0.3181932 -4.043977 4.680363 -6.3531698 6.989556kable(accuracy(hwm1_f, test.ts1), align = "c")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -1.050813 | 8.343387 | 4.624479 | -24.93227 | 49.10704 | 0.5075647 | 0.1036556 | NA |
| Test set | 11.352207 | 12.883862 | 11.352207 | 73.60455 | 73.60455 | 1.2459739 | 0.3208982 | 2.448427 |
checkresiduals(hwm1_f)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,A,M)
## Q* = 36.453, df = 3, p-value = 6.007e-08
##
## Model df: 16. Total lags used: 19Nota: Optimizar el modelo AAA para el parámetro beta