A diferencia de los métodos naive y de medias móviles, el suavizamiento exponencial le da más peso a observaciones más recientes.

library("tidyverse")
library("fpp2")
# create training and validation of the Google stock data
goog.train <- window(goog, end = 900)
goog.test <- window(goog, start = 901)

# create training and validation of the AirPassengers data
qcement.train <- window(qcement, end = c(2012, 4))
qcement.test <- window(qcement, start = c(2013, 1))

Suaviamiento exponencial simple - SES (simple exponencial smoothing)

Este procedimiento es adecuado cuando los datos no tienen tendencia o patrón estacional. Los pesos de cada observación son determinados por un parametro de suavizamiento \(\alpha\). Para un conjunto de datos con \(T\) observaciones , calculamos el valor predicho \(\hat{y_{t+1}}\), el cual estará basado en \(y_{1}\) a través de \(y_{t}\) de la siguiente forma:

\[ \hat{y_{t+1}} = \alpha y_{t} + \alpha(1- \alpha)y_{t-1} + ... +\alpha(1-\alpha)^{t-1}y_{1} \] Donde \(0 < \alpha < 1\). Este \(\alpha\) puede ser visto como una tasa de aprendizaje. Valores cercanos a cero son considerados como aprendizaje lento ya que se da más peso a información histórica, mientra que valores cercanos a 1 son considerados como aprendizaje rápido porque el algoritmo da más peso a las observaciones recientes.

ses.goog <- ses(goog.train, alpha = .2, h = 100)
autoplot(ses.goog) + theme_classic()

Como se comentó antes, los modelos SES no son utiles cuando las series tienen tendencia o comportamiento estacional, por ello el pronóstico a 100 dias de la serie es una linea plana.

Podemos probar removiendo la tendencia de la series diferenciandola

goog.dif <- diff(goog.train)
autoplot(goog.dif)

ses.goog.dif <- ses(goog.dif, alpha = .2, h = 100)
autoplot(ses.goog.dif)

Podemos evaluar el desempeño del modelo en nuestro set de validación pero primero debemos diferenciarla, al igual que set con el que se entrenó el modelo.

goog.dif.test <- diff(goog.test)
accuracy(ses.goog.dif, goog.dif.test)
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.01368221 9.317223 6.398819  99.97907 253.7069 0.7572009
## Test set      0.97219517 8.141450 6.117483 109.93320 177.9684 0.7239091
##                     ACF1 Theil's U
## Training set -0.05440377        NA
## Test set      0.12278141 0.9900678

Podemos hacer tuning de la tasa de aprendizaje \(\alpha\) buscando aquel parametro que disminuya el RMSE

# identify optimal alpha parameter
alpha <- seq(.01, .99, by = .01)
RMSE <- NA
for(i in seq_along(alpha)) {
  fit <- ses(goog.dif, alpha = alpha[i], h = 100)
  RMSE[i] <- accuracy(fit, goog.dif.test)[2,2]
}

# convert to a data frame and idenitify min alpha value
alpha.fit <- tibble(alpha, RMSE)
alpha.min <- filter(alpha.fit, RMSE == min(RMSE))

# plot RMSE vs. alpha
ggplot(alpha.fit, aes(alpha, RMSE)) +
  geom_line() +
  geom_point(data = alpha.min, aes(alpha, RMSE), size = 2, color = "blue")

Corremos el modelo con nuestro hiperparametro optimizado

# refit model with alpha = .05
ses.goog.opt <- ses(goog.dif, alpha = .05, h = 100)

# performance eval
accuracy(ses.goog.opt, goog.dif.test)
##                       ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.01188991 8.939340 6.030873 109.97354 155.7700 0.7136602
## Test set      0.30483955 8.088941 6.028383  97.77722 112.2178 0.7133655
##                    ACF1 Theil's U
## Training set 0.01387261        NA
## Test set     0.12278141 0.9998811
# plotting results
p1 <- autoplot(ses.goog.opt) +
  theme(legend.position = "bottom")
p2 <- autoplot(goog.dif.test) +
  autolayer(ses.goog.opt, alpha = .5) +
  ggtitle("Predicted vs. actuals for the test data set")

gridExtra::grid.arrange(p1, p2, nrow = 1)

Método de Holt

El modelo SES no funciona tampoco con predicciones de un plazo largo que además tiene tendencia. Una alternativa es el método de holt, que captura la tendencia de los datos mientras aplica suavizamiento exponencial

El método de holt hace predicciones sobre datos con tendencia usando dos parametros de suavizamiento \(\alpha\) y \(\beta\) que corresponenden al componente de nivel y tendencia respectivamente.

Las predicciones con esta metodología es una combinación del nivel estimado en \(t (L_{t})\) y la estimación de la tendencia en \(t (T_{t})\)

\[ \hat{y_{T+1}} = L_{t} +kT_{t} \] Los valores de \(L_{t}\) y \(T_{t}\) son actualizados de la siguiente forma:

\[ 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} \]

La primera ecuación indica que el nivel del tiempo \(t\) es un promedio ponderado del valor actual en el tiempo \(t\) y el nivel en el periodo previo, ajustado por la tendencia.

La segunda ecuación indica que la tendencia en el tiempo \(t\) es un promedio ponderado de la tendencia en el periodo previo y la información más reciente del cambio de nivel. Similar al SES, \(\alpha\) y \(\beta\) están restringidos al rango 0-1.

Para capturar el la tendencia multiplicativa se hace un ajuste menor a las ecuacions anteriores:

\[ \hat{y_{T+1}} = L_{t} kT_{t} \]

\[ L_{t} = \alpha y_{t} +\alpha(1- \alpha)(L_{t-1} T_{t-1}) \]

\[ T_{t} = \beta(\frac{L_{t}}{L_{t-1}}) +(1- \beta)T_{t-1} \]

Podemos volver a aplicar el método de Holt para la serie de Google.

holt.goog <- holt(goog.train, h = 100)
autoplot(holt.goog)

El modelo se optimiza solo minimizando los criterios AIC y BIC, encontrando los parametros de suavizamiento que generan las mejores predicciones

holt.goog$model
## Holt's method 
## 
## Call:
##  holt(y = goog.train, h = 100) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 401.1276 
##     b = 0.4091 
## 
##   sigma:  8.8149
## 
##      AIC     AICc      BIC 
## 10045.74 10045.81 10069.75
# errores de predicción
accuracy(holt.goog, goog.test)
##                        ME      RMSE       MAE         MPE     MAPE     MASE
## Training set -0.003332796  8.795267  5.821057 -0.01211821 1.000720 1.002452
## Test set      0.545744415 16.328680 12.876836  0.03013427 1.646261 2.217538
##                    ACF1 Theil's U
## Training set 0.03100836        NA
## Test set     0.87733298  2.024518

Sin embargo la optimzación está basada en el set de entrenamiento y no en el de testeo, asi que podemos realizar una optimización sobre el error de pronóstico

# identify optimal alpha parameter
beta <- seq(.0001, .5, by = .001)
RMSE <- NA
for(i in seq_along(beta)) {
  fit <- holt(goog.train, beta = beta[i], h = 100)
  RMSE[i] <- accuracy(fit, goog.test)[2,2]
}

# convert to a data frame and idenitify min alpha value
beta.fit <- data_frame(beta, RMSE)
## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
beta.min <- filter(beta.fit, RMSE == min(RMSE))

# plot RMSE vs. alpha
ggplot(beta.fit, aes(beta, RMSE)) +
  geom_line() +
  geom_point(data = beta.min, aes(beta, RMSE), size = 2, color = "blue")

# new model with optimal beta
holt.goog.opt <- holt(goog.train, h = 100, beta = 0.0601)

# accuracy of first model
accuracy(holt.goog, goog.test)
##                        ME      RMSE       MAE         MPE     MAPE     MASE
## Training set -0.003332796  8.795267  5.821057 -0.01211821 1.000720 1.002452
## Test set      0.545744415 16.328680 12.876836  0.03013427 1.646261 2.217538
##                    ACF1 Theil's U
## Training set 0.03100836        NA
## Test set     0.87733298  2.024518
# accuracy of new optimal model
accuracy(holt.goog.opt, goog.test)
##                       ME      RMSE       MAE          MPE     MAPE     MASE
## Training set -0.01493114  8.960214  6.058869 -0.005524151 1.039572 1.043406
## Test set     21.41138275 28.549029 23.841097  2.665066997 2.988712 4.105709
##                     ACF1 Theil's U
## Training set 0.009696325        NA
## Test set     0.895371665  3.435763

Optimizar beta trae algunos efectos secundarios sobre los intervalos de confianza calculados ya que al incremental hace que la tasa de aprendizaje sea más rápida y por lo tanto causando más incertidumbre.

p1 <- autoplot(holt.goog) +
  ggtitle("Original Holt's Model") +
  coord_cartesian(ylim = c(400, 1000))

p2 <- autoplot(holt.goog.opt) +
  ggtitle("Optimal Holt's Model") +
  coord_cartesian(ylim = c(400, 1000))

gridExtra::grid.arrange(p1, p2, nrow = 1)

Método holt winters

Para hacer predicciones sobre datos con tendencia y comportamiento estacional podemos usar el método holt winters. Este método se pude implementar de forma aditiva o mutliplicativa , donde la elección dependerá del dataset. El modelo aditivo es mejor cuando la tendencia estacional es de la misma magnitud a lo largo de todos los datos, mientras que el modelo aditivo es preferido cuando la magnitud de la estacionalidad cambia cuando el tiempo incrementa.

autoplot(decompose(qcement))

modelo aditivo

Para el modelo aditivo la ecuación es la siguiente:

\[ \hat{y_{t+1}} = L_{t} +kT_{t} + S_{t+k-m} \] El nivel \(L_{t}\), tendencia \(T_{t}\) y la estacionalidad \(S_{t}\) son actualizados a partir de las siguientes ecuaciones con sus correspondientes parametros:

\[ L_{t} = \alpha (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 (y_{t} - L_{t}) + (1-\gamma)S_{t-m} \]

ETS (Error, trend, seasonality)

qcement.hw <- ets(qcement.train, model = "AAA")
autoplot(forecast(qcement.hw))

En este modelo podemos especificar el error, tendencia y estacionalidad de la serie de la siguiente forma:

En R, por ejemplo model = "AAN, indica que el error y tencia es aditiva y no existe estacionalidad. Si queremos buscar el mejor modelo model = "ZZZ

summary(qcement.hw)
## ETS(A,A,A) 
## 
## Call:
##  ets(y = qcement.train, model = "AAA") 
## 
##   Smoothing parameters:
##     alpha = 0.6418 
##     beta  = 1e-04 
##     gamma = 0.1988 
## 
##   Initial states:
##     l = 0.4511 
##     b = 0.0075 
##     s = 0.0049 0.0307 9e-04 -0.0365
## 
##   sigma:  0.0854
## 
##      AIC     AICc      BIC 
## 126.0419 126.8676 156.9060 
## 
## Training set error measures:
##                       ME       RMSE       MAE          MPE     MAPE      MASE
## Training set 0.001463693 0.08393279 0.0597683 -0.003454533 3.922727 0.5912949
##                    ACF1
## Training set 0.02150539
checkresiduals(qcement.hw)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,A,A)
## Q* = 20.288, df = 3, p-value = 0.0001479
## 
## Model df: 8.   Total lags used: 11

si revisamos los errores del modelo, tienden a crecer a medida que pasa el tiempo, por lo que la medida más apropiada podría ser de error multiplicativo.

qcement.f1 <- forecast(qcement.hw, h = 5)

# check accuracy
accuracy(qcement.f1, qcement.test)
##                       ME       RMSE        MAE          MPE     MAPE      MASE
## Training set 0.001463693 0.08393279 0.05976830 -0.003454533 3.922727 0.5912949
## Test set     0.031362775 0.07144211 0.06791904  1.115342984 2.899446 0.6719311
##                     ACF1 Theil's U
## Training set  0.02150539        NA
## Test set     -0.31290496 0.2112428

modelo multiplicativo

\[ \hat{y_{t+1}} = (L_{t} +kT_{t})S_{t+k-m} \] El nivel \(L_{t}\), tendencia \(T_{t}\) y la estacionalidad \(S_{t}\) son actualizados a partir de las siguientes ecuaciones con sus correspondientes parametros:

\[ L_{t} = \alpha (\frac{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 (\frac{y_{t}}{L_{t}}) + (1-\gamma)S_{t-m} \]

qcement.hw2 <- ets(qcement.train, model = "MAM")
checkresiduals(qcement.hw2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,A,M)
## Q* = 23.433, df = 3, p-value = 3.281e-05
## 
## Model df: 8.   Total lags used: 11

Comparación de distintos modelos

# additive error, trend and seasonality
qcement.hw1 <- ets(qcement.train, model = "AAA")
qcement.f1 <- forecast(qcement.hw1, h = 5)
accuracy(qcement.f1, qcement.test)
##                       ME       RMSE        MAE          MPE     MAPE      MASE
## Training set 0.001463693 0.08393279 0.05976830 -0.003454533 3.922727 0.5912949
## Test set     0.031362775 0.07144211 0.06791904  1.115342984 2.899446 0.6719311
##                     ACF1 Theil's U
## Training set  0.02150539        NA
## Test set     -0.31290496 0.2112428
# multiplicative error, additive trend and seasonality
qcement.hw2 <- ets(qcement.train, model = "MAA")
qcement.f2 <- forecast(qcement.hw2, h = 5)
accuracy(qcement.f2, qcement.test)
##                        ME       RMSE        MAE        MPE     MAPE      MASE
## Training set -0.002700842 0.08387446 0.05912640 -0.3689856 3.811958 0.5849446
## Test set      0.018046198 0.06438073 0.05763575  0.5614035 2.493544 0.5701974
##                      ACF1 Theil's U
## Training set  0.001490172        NA
## Test set     -0.344335948 0.1722734
# additive error and trend and multiplicative seasonality
qcement.hw3 <- ets(qcement.train, model = "AAM", restrict = FALSE)
qcement.f3 <- forecast(qcement.hw3, h = 5)
accuracy(qcement.f3, qcement.test)
##                       ME       RMSE        MAE       MPE     MAPE      MASE
## Training set 0.007542179 0.07818508 0.05640602 0.3985457 3.693964 0.5580315
## Test set     0.046793468 0.09456029 0.09241462 1.6559245 3.882433 0.9142688
##                     ACF1 Theil's U
## Training set -0.01243122        NA
## Test set     -0.05070111 0.2737314
# multiplicative error, additive trend, and multiplicative seasonality
qcement.hw4 <- ets(qcement.train, model = "MAM")
qcement.f4 <- forecast(qcement.hw4, h = 5)
accuracy(qcement.f4, qcement.test)
##                        ME       RMSE        MAE         MPE     MAPE      MASE
## Training set 0.0009645886 0.07807481 0.05595930 -0.02540554 3.657772 0.5536120
## Test set     0.0164979314 0.08297241 0.07162921  0.37381187 3.064455 0.7086363
##                     ACF1 Theil's U
## Training set -0.03038457        NA
## Test set     -0.06185361 0.1930288
qcement.hw5 <- ets(qcement.train, model = "ZZZ")
summary(qcement.hw5)
## ETS(M,A,M) 
## 
## Call:
##  ets(y = qcement.train, model = "ZZZ") 
## 
##   Smoothing parameters:
##     alpha = 0.7664 
##     beta  = 0.0131 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 0.49 
##     b = 0.0064 
##     s = 1.0298 1.0488 1.0151 0.9063
## 
##   sigma:  0.0477
## 
##        AIC       AICc        BIC 
##  0.3161397  1.1418277 31.1802503 
## 
## Training set error measures:
##                        ME       RMSE       MAE         MPE     MAPE     MASE
## Training set 0.0009645886 0.07807481 0.0559593 -0.02540554 3.657772 0.553612
##                     ACF1
## Training set -0.03038457

Podemos optimizar el parametro gamma, dado que la función ets como vimos antes con otras funciones minimizar el RMSE, aIC, BIC de los errores del set de entrenamiento

gamma <- seq(0.01, 0.85, 0.01)
RMSE <- NA

for(i in seq_along(gamma)) {
  hw.expo <- ets(qcement.train, "AAA", gamma = gamma[i])
  future <- forecast(hw.expo, h = 5)
  RMSE[i] = accuracy(future, qcement.test)[2,2]
}

error <- data_frame(gamma, RMSE)
minimum <- filter(error, RMSE == min(RMSE))
ggplot(error, aes(gamma, RMSE)) +
  geom_line() +
  geom_point(data = minimum, color = "blue", size = 2) +
  ggtitle("gamma's impact on forecast errors",
          subtitle = "gamma = 0.21 minimizes RMSE")

# previous model with additive error, trend and seasonality
accuracy(qcement.f1, qcement.test)
##                       ME       RMSE        MAE          MPE     MAPE      MASE
## Training set 0.001463693 0.08393279 0.05976830 -0.003454533 3.922727 0.5912949
## Test set     0.031362775 0.07144211 0.06791904  1.115342984 2.899446 0.6719311
##                     ACF1 Theil's U
## Training set  0.02150539        NA
## Test set     -0.31290496 0.2112428
# new model with optimal gamma parameter
qcement.hw6 <- ets(qcement.train, model = "AAA", gamma = 0.21)
qcement.f6 <- forecast(qcement.hw6, h = 5)
accuracy(qcement.f6, qcement.test)
##                        ME       RMSE        MAE        MPE     MAPE      MASE
## Training set -0.001312025 0.08377557 0.05905971 -0.2684606 3.834134 0.5842847
## Test set      0.033492771 0.07148708 0.06775269  1.2096488 2.881680 0.6702854
##                     ACF1 Theil's U
## Training set  0.04832198        NA
## Test set     -0.35877010 0.2202448
autoplot(qcement.f6)