Capítulo 6 Suavização exponencial
Métodos de suavização exponencial produzem previsões a partir de médias ponderadas de observações passadas, onde o peso associado a cada observação cai a medida em que se recua mais no tempo. Ou seja, quanto mais recente a observação, maior será seu peso no modelo preditivo. Apesar de simples, é usado em larga escala nas mais diversas aplicações.
Para este capítulo, será usada como exemplo a série temporal da cotação do dólar, obtida com o pacote quantmod:
library(tidyverse)
library(forecast)
library(lubridate)
energia <- readRDS("data/ts_energia.Rda")6.1 Suavização simples
A suavização simples exponencial é considerada útil para séries sem tendência ou sazonalidade. No R, ela é implementada pela função forecast::ses()
Levando em conta que o modelo ingênuo (naive model) atribui peso 1 para a última observação e o modelo da média simples atribui peso igual para todas as observações passadas, a suavização simples poderia ser descrita como um meio termo entre ambos. Sua formulação matemática não é complexa:
\(\hat{y}_{T+1|T} = \alpha y_t + \alpha (1 - \alpha)y_{T-1} + \alpha (1 - \alpha)^2y_{T-2} + ...\)
O parâmetro \(\alpha\) é chamado de parâmetro de suavização e está definido no intervalo de 0 a 1. Por exemplo:
alpha <- 0.2
for (i in 1:5) print((1 - alpha)^i)## [1] 0.8
## [1] 0.64
## [1] 0.512
## [1] 0.4096
## [1] 0.32768alpha <- 0.8
for (i in 1:5) print((1 - alpha)^i)## [1] 0.2
## [1] 0.04
## [1] 0.008
## [1] 0.0016
## [1] 0.00032Percebe-se pelos resultados da simulação acima que quanto maior o parâmetro \(\alpha\), maior é o peso dado à observação imediatamente mais recente e menor o dado às demais.
O valor de \(\alpha\) pode ser “definido” subjetivamente, utilizando conhecimentos empíricos. Contudo, a maneira mais precisa de escolher esse valor é por meio de um algoritmo de otimização, que estimará \(\alpha\) a partir dos dados obtidos
Suponha que não façamos a mínima ideia do melhor valor de \(\alpha\) para a série temporal da cotação do dólar. Vamos testar três valores:
alpha1 <- ses(energia, alpha = 0.1, h = 6)
alpha2 <- ses(energia, alpha = 0.5, h = 6)
alpha3 <- ses(energia, alpha = 0.9, h = 6)
# calculando o erro de cada ajuste
list(alpha1, alpha2, alpha3) %>% map(accuracy)## [[1]]
## ME RMSE MAE MPE MAPE MASE
## Training set 196.7656 569.4256 454.5523 2.018782 4.542092 0.8133553
## ACF1
## Training set 0.780361
##
## [[2]]
## ME RMSE MAE MPE MAPE MASE
## Training set 41.95123 376.6313 275.7089 0.4087944 2.770339 0.4933411
## ACF1
## Training set 0.4307428
##
## [[3]]
## ME RMSE MAE MPE MAPE MASE
## Training set 23.78307 340.4885 252.1677 0.2264939 2.54179 0.4512175
## ACF1
## Training set 0.07867332plot(alpha1, plot.conf=FALSE, ylab = "", main="", fcol="white")
lines(fitted(alpha1), col="blue")
lines(fitted(alpha2), col="red")
lines(fitted(alpha3), col="green")
lines(alpha1$mean, col="blue", type="o")
lines(alpha2$mean, col="red", type="o")
lines(alpha3$mean, col="green", type="o")
legend("topleft",lty=1, col=c(1,"blue","red","green"),
c("serie original", expression(alpha == 0.1),
expression(alpha == 0.5),
expression(alpha == 9)),
pch=1)
# qual o valor otimo encontrado para alpha nesse caso?
ses(energia) %>% summary##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = energia)
##
## Smoothing parameters:
## alpha = 0.9883
##
## Initial states:
## l = 4214.9385
##
## sigma: 339.3383
##
## AIC AICc BIC
## 8262.360 8262.412 8274.780
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 21.70817 339.3383 254.3178 0.2057238 2.571917 0.4550647
## ACF1
## Training set -0.003267389
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Sep 2017 14169.41 13734.53 14604.29 13504.32 14834.50
## Oct 2017 14169.41 13558.00 14780.83 13234.33 15104.49
## Nov 2017 14169.41 13422.06 14916.77 13026.43 15312.40
## Dec 2017 14169.41 13307.29 15031.53 12850.91 15487.91
## Jan 2018 14169.41 13206.10 15132.72 12696.16 15642.67
## Feb 2018 14169.41 13114.58 15224.25 12556.18 15782.64
## Mar 2018 14169.41 13030.38 15308.44 12427.42 15911.41
## Apr 2018 14169.41 12952.00 15386.83 12307.54 16031.29
## May 2018 14169.41 12878.36 15460.46 12194.92 16143.90
## Jun 2018 14169.41 12808.71 15530.12 12088.39 16250.43O valor de \(\alpha\) encontrado ffoi de 0,9999, praticamente um modelo ingênuo.
6.2 Linear de Holt
Holt criou uma extensão ao método de suavização simples que permite prever dados com tendência que possui dois parâmetros \(\alpha\) e \(\beta\). Matematicamente, temos:
\(\hat{y}_{t+h}=l_t + hT_t\)
\(l_t = \alpha y_t + (1 - \alpha)(l_{t-1} + T_{t-1})\)
\(T_t = \beta (l_t - l_{t-1}) + (1 - \beta)T_{t-1}\)
Onde \(T_t\) corresponde a uma estimativa do componente de tendência e \(l_t\) uma estimativa do componente de nível da série no período \(t\). Assim como \(\alpha\), o parâmetro \(\beta\) também está definido no intervalo [0,1].
Recomenda-se que \(l_0\) e \(T_0\) sejam inicializados como \(y_1\) e \(y_2 - y_1\), respectivamente.
No R, a função para aplicar o modelo linear de Holt é forecast::holt(). Os parâmetros podem ser impostos manualmente ou calculados automaticamente por otimização:
mod1 <- holt(energia, alpha = 0.6, beta = 0.4)
mod2 <- holt(energia)
mod2$model## Holt's method
##
## Call:
## holt(y = energia)
##
## Smoothing parameters:
## alpha = 0.9833
## beta = 1e-04
##
## Initial states:
## l = 4222.5274
## b = 21.5091
##
## sigma: 338.6585
##
## AIC AICc BIC
## 8264.499 8264.630 8285.199plot(mod1)
lines(fitted(mod1), col = "blue")
lines(fitted(mod2), col = "red")
# calculando a qualidade de ajuste
list(mod1, mod2) %>% map(accuracy)## [[1]]
## ME RMSE MAE MPE MAPE MASE
## Training set 0.485427 417.464 305.8195 -0.0131544 3.009725 0.5472196
## ACF1
## Training set 0.2376248
##
## [[2]]
## ME RMSE MAE MPE MAPE MASE
## Training set -0.1987728 338.6585 252.7144 -0.03022627 2.557064 0.4521957
## ACF1
## Training set 0.0013771586.3 Holt-winter aditivo e multiplicativo
Uma evolução do modelo linear de Holt foi criado por Holt e Winter para possibilitar a modelagem de séries temporais por suavização exponencial que também possuam um componente sazonal. O método de Holt-Winters possui três equações para calcular os componentes \(l_t\) de nível, \(T_t\) de tendência e \(s_t\) de sazonalidade, com os parâmetros \(\alpha\), \(\beta\) e \(\gamma\).
Esse método possui duas variações, que dependem da natureza do componente sazonal. O método aditivo é preferido quand as variações sazonais são razoavelmente constantes por toda a série, enquanto o multiplicativo pode ser usado quando as variações sazonais são proporcionais à mudança do nível da série.
A formulação matemática completa, um pouco mais complexa que o modelo linear de Holt, pode ser encontrada aqui.
No R, este método é implementado pela função forecast::hw():
# vamos testar tanto o metodo aditivo quanto o multiplicativo para a serie de exemplo
ajuste_ad <- hw(energia, seasonal = "additive")
ajuste_mult <- hw(energia, seasonal = "multiplicative")
plot(energia)
lines(fitted(ajuste_ad), col = "blue")
lines(fitted(ajuste_mult), col = "red")
# calculando a qualidade de ajuste
list(ajuste_ad, ajuste_mult) %>% map(accuracy)## [[1]]
## ME RMSE MAE MPE MAPE MASE
## Training set -2.659188 262.4393 192.8352 -0.04155619 2.067298 0.3450506
## ACF1
## Training set 0.03984058
##
## [[2]]
## ME RMSE MAE MPE MAPE MASE
## Training set -23.04857 261.8311 191.8186 -0.2690185 2.012541 0.3432315
## ACF1
## Training set 0.030588346.4 Seleção automática do melhor modelo de suavização exponencial
Além dos apresentados neste capítulo, existem muito mais variações de métodos de suavização exponencial. São 15 no total, que são:
Felizmente, o pacote forecast traz uma função que automatiza internamente a seleção do melhor método de previsão, através da função ets():
# ajustando um modelo
modelo.ets <- ets(energia)
# verificando o output
summary(modelo.ets)## ETS(A,N,A)
##
## Call:
## ets(y = energia)
##
## Smoothing parameters:
## alpha = 0.8441
## gamma = 0.1104
##
## Initial states:
## l = 5072.4965
## s=-165.4613 148.9696 169.1238 228.8757 266.4089 0.6793
## 53.0921 74.6406 113.6357 -255.5424 -281.5381 -352.8838
##
## sigma: 262.6268
##
## AIC AICc BIC
## 8048.548 8049.619 8110.646
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 22.05688 262.6268 195.692 0.2044575 2.100431 0.3501625
## ACF1
## Training set 0.01929531Para a série temporal de exemplo, a função retornou um modelo ETS(A, N, A). A primeira letra se refere ao componente de erro e pode ser A (aditivo) ou M (multiplicativo), a segunda ao componente de tendência e pode ser N (não possui), A (aditivo), Ad (aditivo amortecido), M (multiplicativo) ou Md (multiplicativo amortecido) e a terceira ao componente de sazonalidade, que pode ser N, A ou M.
Ou seja, o algoritmo da função detectou que a série de exemplo possui componente de erro aditivo, não possui tendência e a sazonalidade é aditiva.
O ajuste do modelo, graficamente, é:
plot(energia)
lines(fitted(modelo.ets), col = "red")
