Séries Temporais: suavização exponencial e acurácia
Adriano Marcos Rodrigues Figueiredo, *e-mail: adriano.figueiredo@ufms.br*
28 agosto 2019
Abstract
This is an undergrad student level instruction for class use.Licença
This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.
License: CC BY-SA 4.0
Citação
Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Séries Temporais: suavização exponencial e acurácia. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/smoothing_varejoms.
1 Introdução
Neste arquivo utilizo a série do Índice de volume de vendas no varejo Total de Mato Grosso do Sul, série mensal a partir de jan/2000 até jun/2019 obtida com o pacote BETS e importada do Banco Central do Brasil. Portanto, são 234 observações mensais.
É um exercício rápido de suavização exponencial e acurácia, similar ao realizado para decomposição clássica e x11 (em http://rpubs.com/amrofi/decompose_x11_varejoms).
2 Dados
O leitor pode carregar direto do site do Banco Central do Brasil com o pacote BETS (FERREIRA, SPERANZA e COSTA, 2018) como abixo ou usar a saída do dput(varejoms) conforme colocado abaixo.
library(BETS)
# Pegando as séries a partir do site do Banco Central do Brasil
# Índice de volume de vendas no varejo Total de Mato Grosso do Sul
# mensal a partir de jan/2000 até jun/2019
# 234 observações mensais
varejoms <- BETSget(1479)
print(varejoms)
class(varejoms)
dput(varejoms) # opção para ter os dados como na structure abaixo# uso a opção abaixo para quando a internet por algum motivo falha
varejoms <- structure(c(35.2, 35.6, 39.2, 40.5, 41.6, 40.4, 40.8, 38.7, 37.3,
37.6, 35.6, 47.4, 34.2, 32.2, 38, 37.5, 38.8, 35, 38.4, 40.2, 38.2, 39.3,
36, 46.3, 36.4, 34, 39, 37.8, 38.9, 35.3, 37.2, 38.1, 35.6, 38.3, 35.6,
45.7, 32.2, 31.6, 35.2, 36.8, 37.5, 34.8, 38.4, 38.1, 37, 39, 37.3, 49,
35.8, 35.3, 40.2, 41.3, 43.9, 41.6, 45.9, 42.3, 42.2, 44, 41.3, 56.9, 38.5,
38.5, 45.3, 43.6, 46.2, 44.3, 47.5, 46.5, 46.4, 46, 44.1, 61, 42, 40.2,
44.6, 44.5, 47.8, 45.3, 46.5, 48.5, 47.7, 50.2, 49.3, 64.5, 47, 46.8, 51,
50.5, 55, 51.3, 52.8, 55.3, 54.8, 55.6, 55.3, 72.2, 54.5, 52.1, 56.2, 57.2,
60.8, 56.1, 61.8, 61.6, 59.8, 63.3, 57.7, 77.4, 61.4, 51.9, 57.3, 57.9,
61.9, 57.3, 61.1, 61.1, 60.6, 65.6, 63.5, 83.1, 64.1, 60.2, 67.8, 67.1,
72.8, 68.5, 71.1, 69.3, 69.9, 71.1, 67.9, 92.7, 67.5, 64.8, 69.1, 69.4,
79.6, 70.2, 73.8, 72.5, 71.3, 75.6, 74.7, 100.8, 79.5, 75.7, 82.4, 78, 84.8,
83.2, 84.8, 88.5, 86.3, 91.7, 92.8, 111.4, 92.8, 83.7, 92.5, 88.3, 93.9,
88.8, 96, 95.9, 93.2, 98.3, 100.5, 128.8, 97.2, 90.2, 94.3, 94.4, 101.1,
92, 96.4, 98.2, 97.6, 105.8, 103.1, 129.6, 99.6, 87.8, 97, 94.8, 98.6, 93.4,
98.4, 96.4, 92.5, 100.6, 97.2, 124.5, 91.5, 85.1, 91.6, 88.5, 92.2, 87.4,
90.5, 88.1, 85.2, 89.4, 93.4, 116.9, 90.8, 84, 89.7, 86.3, 90, 87.3, 90.8,
93.5, 93.7, 91.4, 93.5, 114.1, 87.8, 81.1, 94.5, 83.2, 89.9, 88.8, 89.3,
93.7, 93.5, 96.3, 101.3, 118.3, 93.8, 85.2, 90, 86.6, 90, 84.5), .Tsp = c(2000,
2019.41666666667, 12), class = "ts")
class(varejoms)[1] "ts"3 Primeiras análises
Inicialmente olharei as estatísticas descritivas da série. Em seguida farei um plot básico da série e o plot pelo pacote dygraphs, útil para ver os pontos de picos e momentos específicos.
# estatisticas basicas
summary(varejoms) Min. 1st Qu. Median Mean 3rd Qu. Max.
31.60 44.02 64.30 67.17 90.00 129.60 # Min. 1st Qu. Median Mean 3rd Qu. Max. 31.60 44.02 64.30 67.17 90.00
# 129.60
# pelo pacote dygraph dá mais opções
library(dygraphs)
dygraph(varejoms)É possível visualizar nos plots acima: sazonalidade (por exemplo, picos em dezembro de cada ano); a tendência aparentemente crescente até 2014 e decresce com a “crise” brasileira; e uma aparente não-estacionariedade (média e variância mudam no tempo). Em outra postagem aplicarei o teste de raiz unitária na série para avaliar a estacionariedade de modo mais explícito.
4 Suavização (ou alisamento) Exponencial Simples (SES) ou ainda, ‘Holt 1 parâmetro’
Este método é apropriado para séries que se movem aleatoriamente em torno de uma média ou nível constante, sem tendência e sem sazonalidade: \(Z_t=μ_t+a_t\) em que \(a_t\) é um ruído branco de média zero e variância constante e \(μ_t\) é um parâmetro localmente constante (nível); \(Z_t\) é a série em análise.
\[ {\bar Z_t} = \alpha {Z_t} + \left( {1 - \alpha } \right){\bar Z_{t - 1}} \]
Em que α é a constante de suavização, 0<α<1. Portanto, similarmente,
\[ {\bar Z_{t - 1}} = \alpha {Z_{t - 1}} + \left( {1 - \alpha } \right){\bar Z_{t - 2}} \]
O que substituindo (2) em (1) e desenvolvendo os termos tem-se:
\[ {\bar Z_t} = \alpha {Z_t} + \alpha \left( {1 - \alpha } \right){Z_{t - 1}} + \alpha {\left( {1 - \alpha } \right)^2}{Z_{t - 2}} + \ldots \]
Ou seja, a SES faz uma previsão que é a média ponderada por \(α(1-α)^k\) que, quanto mais recente, maior é o peso.
Se \(0<α<1\), então \(α>α(1-α)>α(1-α)^2>\ldots>α(1-α)^k\)
Na MMS tem-se o mesmo peso entre observações recentes e antigas. Na SES, a vantagem sobre MMS é que aplica pesos maiores às observações recentes. A questão será achar α que minimiza a soma de quadrados de ajustamento:
\[ S = \mathop \sum \limits_{t = l + 1}^N {\left( {{Z_t} - {{\hat Z}_{t - 1}}} \right)^2} \]
O forecast será constante para todas as observações futuras!!!
x <- varejoms
# simple exponential smoothing (SES)
x.ses <- ses(x, h = 12, level = c(95))
x.ses # exibe o forecast por holt 1 parametro Point Forecast Lo 95 Hi 95
Jul 2019 90.3704 75.71475 105.0260
Aug 2019 90.3704 75.36817 105.3726
Sep 2019 90.3704 75.02941 105.7114
Oct 2019 90.3704 74.69797 106.0428
Nov 2019 90.3704 74.37340 106.3674
Dec 2019 90.3704 74.05528 106.6855
Jan 2020 90.3704 73.74325 106.9975
Feb 2020 90.3704 73.43696 107.3038
Mar 2020 90.3704 73.13612 107.6047
Apr 2020 90.3704 72.84045 107.9004
May 2020 90.3704 72.54967 108.1911
Jun 2020 90.3704 72.26357 108.4772summary(x.ses) # fornece o relatorio da estimacao
Forecast method: Simple exponential smoothing
Model Information:
Simple exponential smoothing
Call:
ses(y = x, h = 12, level = c(95))
Smoothing parameters:
alpha = 0.2188
Initial states:
l = 37.9823
sigma: 7.4775
AIC AICc BIC
2222.106 2222.210 2232.471
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set 1.0234 7.445484 4.746198 0.8232233 6.805751 1.031984
ACF1
Training set 0.03393904
Forecasts:
Point Forecast Lo 95 Hi 95
Jul 2019 90.3704 75.71475 105.0260
Aug 2019 90.3704 75.36817 105.3726
Sep 2019 90.3704 75.02941 105.7114
Oct 2019 90.3704 74.69797 106.0428
Nov 2019 90.3704 74.37340 106.3674
Dec 2019 90.3704 74.05528 106.6855
Jan 2020 90.3704 73.74325 106.9975
Feb 2020 90.3704 73.43696 107.3038
Mar 2020 90.3704 73.13612 107.6047
Apr 2020 90.3704 72.84045 107.9004
May 2020 90.3704 72.54967 108.1911
Jun 2020 90.3704 72.26357 108.4772plot(x.ses, col = 1, main = "Índice de volume de vendas no varejo Total
de Mato Grosso do Sul",
sub = "modelo SES")
5 Suavização Holt-Winters - 3 parâmetros
# Holt-Winters
x.hwm <- hw(x, seasonal = "multiplicative") # multiplicativo
plot(x, ylab = "Indice 2011=100", main = "Índice de volume de vendas no varejo Total de Mato Grosso do Sul",
type = "o", xlab = "Mês/Ano")
lines(fitted(x.ses), col = "red", lty = 2)
lines(x.ses$mean, type = "o", col = "red")
lines(x.hwm$mean, type = "o", col = "green")
legend("topleft", lty = 1, pch = 1, col = 1:3, c("original", "Holt Simples",
"Holt Winters' Multiplicativo"))
summary(x.hwm)
Forecast method: Holt-Winters' multiplicative method
Model Information:
Holt-Winters' multiplicative method
Call:
hw(y = x, seasonal = "multiplicative")
Smoothing parameters:
alpha = 0.4597
beta = 0.0192
gamma = 0.2855
Initial states:
l = 41.5931
b = -4e-04
s = 1.2461 0.9587 1.0374 1.0132 1.0348 1.0325
0.9809 1.0255 0.9906 0.968 0.8487 0.8637
sigma: 0.0348
AIC AICc BIC
1658.395 1661.228 1717.135
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.04501476 2.32154 1.76375 0.0133938 2.665347 0.383499
ACF1
Training set 0.03646486
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jul 2019 88.88771 84.92468 92.85075 82.82678 94.94865
Aug 2019 90.62694 86.14540 95.10848 83.77302 97.48086
Sep 2019 88.92541 84.10245 93.74837 81.54933 96.30149
Oct 2019 91.47520 86.08043 96.86997 83.22461 99.72579
Nov 2019 93.50672 87.55113 99.46230 84.39844 102.61500
Dec 2019 113.71305 105.93485 121.49125 101.81732 125.60878
Jan 2020 88.48556 82.01555 94.95558 78.59053 98.38060
Feb 2020 80.83252 74.53942 87.12562 71.20805 90.45698
Mar 2020 87.79739 80.54427 95.05050 76.70471 98.89006
Apr 2020 82.93359 75.68515 90.18204 71.84805 94.01913
May 2020 87.87547 79.77122 95.97972 75.48108 100.26985
Jun 2020 84.38425 76.19165 92.57686 71.85475 96.91376
Jul 2020 87.80739 78.44802 97.16675 73.49348 102.12130
Aug 2020 89.52432 79.55703 99.49161 74.28067 104.76797
Sep 2020 87.84235 77.63994 98.04476 72.23911 103.44559
Oct 2020 90.35991 79.42497 101.29484 73.63637 107.08345
Nov 2020 92.36545 80.73216 103.99874 74.57387 110.15703
Dec 2020 112.32369 97.61559 127.03180 89.82959 134.81780
Jan 2021 87.40329 75.51629 99.29029 69.22370 105.58288
Feb 2021 79.84280 68.57524 91.11036 62.61056 97.07504
[ reached 'max' / getOption("max.print") -- omitted 4 rows ]6 Avaliação da simulação
Para uma só equação, podem-se olhar os indicadores tradicionais como o R², estatísticas F e t, critérios de informação de Akaike, Schwarz entre outros. Para modelos de previsão, devo ter um erro de previsão mínimo, definido pelo pesquisador de modo a julgá-lo adequado.
As vezes pode-se ter equações bem estimadas junto de outras com poucos parâmetros significativos, num mesmo sistema. Dependendo do objetivo a ser alcançado, se simular ou explicar, busca-se ajustamento com alta significância (para explicar) ou busca-se melhor poder de explicação. Por melhor poder de explicação entende-se: valores previstos mais próximos dos observados.
Alguns critérios para avaliação da previsão podem ser citados: Raiz do Erro de Simulação Quadrático Médio; Raiz do Erro Quadrático Médio Percentual; Erro de Simulação Médio; Erro Percentual Médio; Erro Absoluto Médio; Erro Absoluto Percentual Médio; Raiz do Erro de Simulação Ex-post Quadrático Médio; Coeficiente de desigualdade de Theil (U).
6.1 Raiz do Erro de Simulação Quadrático Médio (RMSE)
Este critério é dado pela expressão: \[ RMSE = \;\sqrt {\frac{1}{T}\mathop \sum \limits_{t = 1}^T {{\left( {Y_t^s - Y_t^a} \right)}^2}} \]
Em que: \(Y_t^s\) é o valor simulado (previsto) de \(Y_t\); \(Y_t^a\) é o valor efetivo (observado); e T é o número de períodos na simulação.
O resultado é um escalar com o valor da raiz quadrada da média da diferença quadrática entre x e y, ou seja, a raiz do erro quadrático médio. Normalmente deve ser calculado pela fórmula, mas é reportado pelo programa quando se faz um “forecast”.
Para o cálculo da Raiz do Erro de Simulação Ex-post Quadrático Médio, utilizam-se os resultados do modelo de previsão aplicado a um conjunto restrito (histórico) para checar os erros no período de simulação. Exemplo: usando os dados de 1950:01 a 1985:04, históricos, simulam-se os resultados para 1986:01 a 1988:01 e calcula-se o erro ex-post.
6.2 Raiz do Erro Quadrático Médio Percentual (PRMSE)
\[ PRMSE = \;\sqrt {\frac{1}{T}\mathop \sum \limits_{t = 1}^T {{\left( {\frac{{Y_t^s - Y_t^a}}{{Y_t^a}}} \right)}^2}} \]
Em que: \(Y_t^s\) é o valor simulado (previsto) de \(Y_t\); \(Y_t^a\) é o valor efetivo (observado); e T é o número de períodos na simulação. Neste caso, a diferença entre o previsto e o observado (o erro) é dividida pelo valor observado, obtendo-se assim em base percentual relativo ao valor observado.
6.3 Erro de Simulação Médio (ME) e o Erro Percentual Médio (MPE)
O erro de simulação é basicamente o termo \((Y_t^s - Y_t^a)\), como definido anteriormente. O erro médio (ME) é, portanto, a média dos erros para T observações, enquanto o erro percentual médio (MPE) é obtido relativizando-se pelos valores observados, como nas expressões a seguir:
\[ ME= \frac{1}{T}\sum\limits_{t = 1}^T {\left( {Y_t^s - Y_t^a} \right)} \]
\[ MPE = 100.\frac{1}{T}\sum\limits_{t = 1}^T {\frac{{\left( {Y_t^s - Y_t^a} \right)}}{{Y_t^a}}} \]
Uma ressalva necessária é que: “[…] O problema com os erros médios é que eles podem ser quase zero, quando elevados erros positivos cancelam elevados erros negativos” (PINDYCK e RUBINFELD, 2004, p.447).
Alguns autores utilizam estas duas últimas medidas com o termo entre parênteses $ {( {Y_t^s - Y_t^a} )} $ na forma de ‘valor absoluto’ \(\left| {\left( {Y_t^s - Y_t^a} \right)} \right|\) . Tal procedimento auxilia a investigar o erro de previsão do modelo evitando os problemas de ‘erros positivos’ cancelarem ‘erros negativos’. Desta forma, surgem as outras duas opções: MAE e MAPE.
6.4 Erro Absoluto Médio (MAE) e Erro Absoluto Percentual Médio (MAPE)
O erro absoluto médio (MAE) considera a expressão análoga do erro de simulação médio com os valores absolutos como na expressão:
\[ MAE = \frac{1}{T}\sum\limits_{t = 1}^T {\left| {Y_t^s - Y_t^a} \right|} \]
Similarmente ao erro percentual médio, o erro absoluto percentual médio (MAPE) será:
\[ MAPE = 100.\frac{1}{T}\sum\limits_{t = 1}^T {\left| {\frac{{\left( {Y_t^s - Y_t^a} \right)}}{{Y_t^a}}} \right|} \]
Desejam-se valores baixos das medidas. Entretanto, é importante observar se o modelo está bem ajustado para os ‘pontos de inflexão’ nos dados (oscilações, vales e picos que se observam nos gráficos previstos e realizados).
6.5 Erro Absoluto Percentual Médio Simétrico (sMAPE)
As medidas de acurácia tipo MAPE têm a desvantagem de penalizar mais os erros negativos que os positivos. Assim, surge a medida sMAPE ou MAPE “simétrico” proposto por Armstrong (1978 , p. 348), citado por Hyndman e Athanasopoulos (2018). Esta medida foi utilizada para escolha do melhor modelo na competição de previsão do M3 de Makridakis (\(M^3\) competition). A expressão para cálculo é do tipo:
\[ \text{sMAPE} = \text{mean}\left(200|{Y_t^s - Y_t^a}|/(Y_t^s+Y_t^a)\right) \]
Esta medida não está disponibilizada diretamente nas saídas do pacote forecast nem no fpp2, em virtude da crítica de Hyndman e Athanasopoulos (2018) em que citam Hyndman e Koehler (2006). A medida sMAPE pode ser obtida pelos pacotes Metrics ou TSPred.
6.6 critérios de informação, Akaike (AIC), Schwarz (SIC ou BIC) e Hannan-Quin (HQ)
A ideia é que nos critérios de informação, acrescentando lags se reduz a soma dos resíduos (pois mais variação é explicada no modelo), ou seja, foi acrescentada mais informação. Mas ainda deve-se penalizar a adição do regressor, uma ideia de redução dos graus de liberdade do modelo.
Em geral, os critérios têm a seguinte expressão: \(C=lnσ ̂^2 (T)+c_T φ(T)\).
O primeiro termo associa a variância, e seu valor deve reduzir com o aumento de variáveis. No segundo termo tem-se: \(c_T\) é o número de parâmetros estimados; e \(φ(T)\) é a ordem do processo que penaliza o acréscimo de variáveis. Assim, quanto menor C melhor é o modelo, ou seja, tem menos resíduos, menos parte não explicada. O Critério de Informação de Akaike (ou AIC de Akaike’s Information Criterion) ou o Critério de Informação de Schwarz ou Bayesiano (ou SIC de Schwarz’s Information Criterion ou em alguns livros BIC de Bayesian Information Criterion) são expressos no Eviews da forma já logaritmizada como:
\[AIC = - \frac{{2l}}{n} + \frac{{2k}}{n}\]
\[SIC = - \frac{{2l}}{n} + \frac{{k.\log n}}{n}\]
\[ l = - {\frac{n}{2}} {\left( {1 + \log(2\pi)+\log \left( {\frac{{\hat \varepsilon ^\prime \hat \varepsilon }}{n}} \right)} \right)} \]
Em que: k é o número de regressores incluindo-se o intercepto; n é o número de observações; l é o log Verossimilhança da regressão; e \(\hat \varepsilon\) são os resíduos estimados do modelo. No formato mais simplificado exposto por Greene (2002), tem-se igual ao de Bueno (2008, p.47):
\[ \begin{array}{l} SIC = BIC\left( {p,q} \right) = \ln {{\hat \sigma }^2} + k\frac{{\ln T}}{T}\\ {{\hat \sigma }^2} = \frac{{\sum\nolimits_{t = 1}^T {\hat \varepsilon _t^2} }}{T} \end{array} \]
Em que k=p+q+1 se houver intercepto.
\[ AIC\left( {p,q} \right) = \ln {{\hat \sigma }^2} + k\frac{{2}}{T} \]
\[ HQ\left( {p,q} \right) = \ln {{\hat \sigma }^2} + k\frac{2}{T} \ln T \]
Quanto mais parâmetros são estimados, menor será o erro estimado, mas isto é penalizado no segundo termo dos critérios. Desejo menores valores dos critérios. O BIC é consistente assintoticamente, e o AIC funciona melhor em pequenas amostras, mas é pior com mais parâmetros. Já o HQ é menos forte que o BIC, ou seja, o preferido é o BIC.
A função window auxiliará na extração de uma parte de uma série temporal, e criação dos subsets de treino e teste no forecast. Na função window especificamos o início (start) (e/ou fim: end) de modo a definir o subset. O treino é onde será obtido o modelo estimado e o teste é onde verificaremos o acerto, ou acurácia, do modelo. Identificado o melhor modelo, faz-se então o forecast fora da amostra. A amostra de teste é em geral 20% da amostra existente.
library(fpp2)
# farei as janelas de teste e treino, estimarei os modelos para a amostra
# treino modelo 1: cons.ses modelo 2: cons.holttrend.EXP modelo 3:
# cons.holttrend.DAMPm modelo 4: conshw2
# amostra treino até julho 2015 (187 meses)
cons <- window(varejoms, start = c(2000, 1), end = c(2015, 7))
# amostra teste (47 meses até junho 2019)
teste <- window(varejoms, start = c(2015, 8))Vamos checar se separamos corretamente:
autoplot(varejoms) + autolayer(cons, series = "Training") + autolayer(teste,
series = "Test")
# estimar modelos
cons.ses <- ses(cons, h = 60, level = c(95)) # forecast até jul/2020
cons.holttrend.EXP <- holt(cons, h = 60, level = c(95), exponential = TRUE)
cons.holttrend.DAMPm <- holt(cons, h = 60, level = c(95), exponential = TRUE,
damped = TRUE)
conshw2 <- hw(cons, h = 60, seasonal = "multiplicative") # multiplicativo
autoplot(varejoms) + forecast::autolayer(cons.ses, series = "SES", PI = FALSE) +
forecast::autolayer(cons.holttrend.EXP, series = "Holt Expon. Trend", PI = FALSE) +
forecast::autolayer(cons.holttrend.DAMPm, series = "Holt Trend Damped mult",
PI = FALSE) + forecast::autolayer(conshw2, series = "Holt-Winters mult",
PI = FALSE) + xlab("Ano") + ylab("Consumo") + ggtitle("Forecasts para consumo do varejo MS") +
guides(colour = guide_legend(title = "Forecast"))
summary(cons.ses)
Forecast method: Simple exponential smoothing
Model Information:
Simple exponential smoothing
Call:
ses(y = cons, h = 60, level = c(95))
Smoothing parameters:
alpha = 0.239
Initial states:
l = 37.8842
sigma: 6.8763
AIC AICc BIC
1703.310 1703.441 1713.003
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set 1.339005 6.839444 4.260657 1.242669 6.788559 0.895096
ACF1
Training set -0.02999918
Forecasts:
Point Forecast Lo 95 Hi 95
Aug 2015 97.73313 84.25580 111.2105
Sep 2015 97.73313 83.87616 111.5901
Oct 2015 97.73313 83.50666 111.9596
Nov 2015 97.73313 83.14651 112.3197
Dec 2015 97.73313 82.79504 112.6712
Jan 2016 97.73313 82.45165 113.0146
Feb 2016 97.73313 82.11581 113.3504
Mar 2016 97.73313 81.78705 113.6792
Apr 2016 97.73313 81.46492 114.0013
May 2016 97.73313 81.14905 114.3172
Jun 2016 97.73313 80.83909 114.6272
Jul 2016 97.73313 80.53471 114.9315
Aug 2016 97.73313 80.23563 115.2306
Sep 2016 97.73313 79.94157 115.5247
Oct 2016 97.73313 79.65229 115.8140
Nov 2016 97.73313 79.36757 116.0987
Dec 2016 97.73313 79.08720 116.3791
Jan 2017 97.73313 78.81098 116.6553
Feb 2017 97.73313 78.53874 116.9275
Mar 2017 97.73313 78.27030 117.1960
Apr 2017 97.73313 78.00552 117.4607
May 2017 97.73313 77.74424 117.7220
Jun 2017 97.73313 77.48633 117.9799
Jul 2017 97.73313 77.23167 118.2346
Aug 2017 97.73313 76.98013 118.4861
Sep 2017 97.73313 76.73161 118.7346
Oct 2017 97.73313 76.48599 118.9803
Nov 2017 97.73313 76.24318 119.2231
Dec 2017 97.73313 76.00308 119.4632
Jan 2018 97.73313 75.76561 119.7006
Feb 2018 97.73313 75.53068 119.9356
Mar 2018 97.73313 75.29820 120.1680
Apr 2018 97.73313 75.06811 120.3981
[ reached 'max' / getOption("max.print") -- omitted 27 rows ]summary(cons.holttrend.EXP)
Forecast method: Holt's method with exponential trend
Model Information:
Holt's method with exponential trend
Call:
holt(y = cons, h = 60, level = c(95), exponential = TRUE)
Smoothing parameters:
alpha = 0.139
beta = 1e-04
Initial states:
l = 38.4459
b = 1.0059
sigma: 0.1038
AIC AICc BIC
1650.810 1651.141 1666.965
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.06616242 6.613945 4.188691 -1.273058 6.816974 0.879977
ACF1
Training set 0.03938391
Forecasts:
Point Forecast Lo 95 Hi 95
Aug 2015 102.6680 82.21325 123.6675
Sep 2015 103.2643 82.12224 125.3514
Oct 2015 103.8640 83.62881 125.2875
Nov 2015 104.4671 83.27280 126.5133
Dec 2015 105.0738 83.27614 127.6882
Jan 2016 105.6840 84.00109 128.7852
Feb 2016 106.2978 84.62298 130.1737
Mar 2016 106.9151 83.86320 131.0781
Apr 2016 107.5360 84.29447 131.9080
May 2016 108.1605 85.15207 132.6212
Jun 2016 108.7886 85.37154 133.6896
Jul 2016 109.4204 85.84993 136.1859
Aug 2016 110.0558 85.83369 135.9517
Sep 2016 110.6950 86.05295 135.8222
Oct 2016 111.3378 85.86740 138.2137
Nov 2016 111.9844 87.25056 139.5991
Dec 2016 112.6347 87.23104 140.2454
Jan 2017 113.2888 87.74143 140.6656
Feb 2017 113.9467 87.84125 143.3614
Mar 2017 114.6085 88.14357 144.0116
Apr 2017 115.2740 89.17252 143.4174
May 2017 115.9435 89.18443 146.0490
Jun 2017 116.6168 89.56823 145.7710
Jul 2017 117.2940 90.09394 147.9986
Aug 2017 117.9752 90.15779 148.5323
Sep 2017 118.6603 91.26682 150.9042
Oct 2017 119.3494 90.97347 150.8392
Nov 2017 120.0426 91.37486 152.2017
Dec 2017 120.7397 93.07605 153.2288
Jan 2018 121.4409 92.49209 153.3283
Feb 2018 122.1461 93.33783 155.9765
Mar 2018 122.8555 93.22401 157.5944
Apr 2018 123.5689 94.07154 158.8590
[ reached 'max' / getOption("max.print") -- omitted 27 rows ]summary(cons.holttrend.DAMPm)
Forecast method: Damped Holt's method with exponential trend
Model Information:
Damped Holt's method with exponential trend
Call:
holt(y = cons, h = 60, damped = TRUE, level = c(95), exponential = TRUE)
Smoothing parameters:
alpha = 0.0835
beta = 0.0169
phi = 0.9572
Initial states:
l = 38.9916
b = 1.0001
sigma: 0.1054
AIC AICc BIC
1652.699 1653.165 1672.085
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.4716645 6.607187 4.182636 0.09037218 6.669775 0.878705
ACF1
Training set 0.06912782
Forecasts:
Point Forecast Lo 95 Hi 95
Aug 2015 100.23830 79.59387 120.7889
Sep 2015 99.95370 79.00841 119.8228
Oct 2015 99.68205 79.08449 120.2223
Nov 2015 99.42273 79.01535 120.4963
Dec 2015 99.17515 78.39902 120.4894
Jan 2016 98.93875 78.12389 120.1448
Feb 2016 98.71300 77.72787 119.5006
Mar 2016 98.49740 77.60502 120.8873
Apr 2016 98.29148 76.60354 121.9379
May 2016 98.09479 76.75972 121.9921
Jun 2016 97.90688 75.93438 121.6110
Jul 2016 97.72737 75.41087 123.2243
Aug 2016 97.55585 74.90334 122.3694
Sep 2016 97.39196 74.49509 123.3998
Oct 2016 97.23534 74.80426 123.2504
Nov 2016 97.08567 73.43970 124.3794
Dec 2016 96.94263 73.25322 124.7477
Jan 2017 96.80591 72.13365 126.3009
Feb 2017 96.67523 72.26753 125.9150
Mar 2017 96.55031 70.99320 126.0826
Apr 2017 96.43090 70.63147 127.9588
May 2017 96.31673 70.12923 127.9068
Jun 2017 96.20758 69.92206 129.8022
Jul 2017 96.10323 69.31071 129.5336
Aug 2017 96.00345 69.44874 130.7782
Sep 2017 95.90804 67.73914 131.9182
Oct 2017 95.81680 67.82567 134.2973
Nov 2017 95.72956 67.24986 134.0381
Dec 2017 95.64613 66.54705 134.8991
Jan 2018 95.56633 65.73807 134.8337
Feb 2018 95.49002 65.50880 136.6154
Mar 2018 95.41704 65.25910 137.3355
Apr 2018 95.34723 64.50632 138.0112
[ reached 'max' / getOption("max.print") -- omitted 27 rows ]summary(conshw2)
Forecast method: Holt-Winters' multiplicative method
Model Information:
Holt-Winters' multiplicative method
Call:
hw(y = cons, h = 60, seasonal = "multiplicative")
Smoothing parameters:
alpha = 0.3785
beta = 0.02
gamma = 0.298
Initial states:
l = 39.9041
b = -0.0938
s = 1.2461 0.9641 1.0256 1.0036 1.0186 1.0309
0.9832 1.0446 0.9987 0.9644 0.8469 0.8733
sigma: 0.0356
AIC AICc BIC
1258.334 1261.956 1313.263
Error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.07561897 2.121197 1.632621 0.2012019 2.748482 0.3429876
ACF1
Training set 0.1300521
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Aug 2015 99.83153 95.27346 104.38959 92.86056 106.80249
Sep 2015 98.12611 93.30452 102.94769 90.75213 105.50009
Oct 2015 104.02632 98.54510 109.50754 95.64352 112.40912
Nov 2015 101.94575 96.20443 107.68707 93.16516 110.72634
Dec 2015 129.15806 121.40670 136.90942 117.30337 141.01274
Jan 2016 99.39022 93.05112 105.72932 89.69540 109.08504
Feb 2016 90.18348 84.08603 96.28093 80.85824 99.50873
Mar 2016 97.95304 90.94864 104.95743 87.24074 108.66534
Apr 2016 95.91397 88.67569 103.15226 84.84397 106.98397
May 2016 102.45680 94.31291 110.60069 90.00180 114.91180
Jun 2016 96.33800 88.28751 104.38849 84.02583 108.65016
Jul 2016 101.66983 92.75319 110.58647 88.03301 115.30666
Aug 2016 102.91979 92.95323 112.88636 87.67725 118.16234
Sep 2016 101.15393 90.94955 111.35831 85.54768 116.76019
Oct 2016 107.22809 95.97031 118.48588 90.01080 124.44539
Nov 2016 105.07558 93.60513 116.54603 87.53303 122.61813
Dec 2016 133.11337 118.01810 148.20863 110.02715 156.19958
Jan 2017 102.42630 90.37071 114.48189 83.98887 120.86373
Feb 2017 92.93143 81.58845 104.27442 75.58383 110.27903
Mar 2017 100.93028 88.16543 113.69514 81.40812 120.45245
[ reached 'max' / getOption("max.print") -- omitted 40 rows ]Acurácia
As previsões acima foram feitas com a amostra treino até julho/2015, mas temos dados até julho/2019. Faremos a acurácia para este período. A função accuracy nos fornece várias medidas de acurácia, e para a sMAPE precisaremos do pacote Metrics e a função smape.
teste <- window(varejoms, start = c(2015, 8))
kable(forecast::accuracy(cons.ses, teste), caption = "SES")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | 1.339006 | 6.839444 | 4.260657 | 1.242669 | 6.788559 | 0.895096 | -0.0299992 | NA |
| Test set | -4.884189 | 10.160182 | 8.684242 | -6.091654 | 9.325527 | 1.824421 | 0.1721972 | 0.9967548 |
kable(forecast::accuracy(cons.holttrend.EXP, teste), caption = "Holt-Trend Exponencial")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | -0.0661624 | 6.613945 | 4.188691 | -1.273058 | 6.816974 | 0.879977 | 0.0393839 | NA |
| Test set | -24.8072759 | 28.235483 | 25.815423 | -27.810865 | 28.630100 | 5.423408 | 0.5854286 | 2.821569 |
kable(forecast::accuracy(cons.holttrend.DAMPm, teste), caption = "Holt-Trend Amortecida Multiplicativo")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | 0.4716645 | 6.607187 | 4.182636 | 0.0903722 | 6.669775 | 0.878705 | 0.0691278 | NA |
| Test set | -3.7017105 | 9.566188 | 7.775848 | -4.7858372 | 8.272414 | 1.633582 | 0.1613323 | 0.9301276 |
kable(forecast::accuracy(conshw2, teste), caption = "Holt-Winters Multiplicativo")| ME | RMSE | MAE | MPE | MAPE | MASE | ACF1 | Theil’s U | |
|---|---|---|---|---|---|---|---|---|
| Training set | 0.075619 | 2.121197 | 1.632621 | 0.2012019 | 2.748482 | 0.3429876 | 0.1300521 | NA |
| Test set | -13.116922 | 14.078510 | 13.116922 | -14.2562156 | 14.256216 | 2.7556560 | 0.5765315 | 1.399466 |
Para o pacote Metrics, observar que ele tem funções de mesmo nome que o pacote forecast, e portanto, deves especificar exatamente qual o pacote e cuidar para alterar a ordem das séries (actual x predicted). No pacote Metrics existe uma diferença da ordem de \(10^2\) comparando com os resultados do pacote forecast!
library(Metrics)
teste <- window(varejoms, start = c(2015, 8))
smape(teste, cons.ses$mean)[1] 0.09066301smape(teste, cons.holttrend.EXP$mean)[1] 0.2447362smape(teste, cons.holttrend.DAMPm$mean)[1] 0.08124576smape(teste, conshw2$mean)[1] 0.1318151Metrics::mape(teste, cons.ses$mean)[1] 0.09325527Metrics::mape(teste, cons.holttrend.EXP$mean)[1] 0.286301Metrics::mape(teste, cons.holttrend.DAMPm$mean)[1] 0.08272414Metrics::mape(teste, conshw2$mean)[1] 0.1425622Comentário O modelo estimado utilizando a amostra de treino perde qualidade em virtude da forte queda após 2015, que a amostra não consegue perceber. Isso é ocasionado em quebras na série, de modo que o forecast ficou prejudicado pela escolha da amostra treino. Observa-se isso no item 5, quando o forecast foi realizado sem delimitar a amostra, realizando a previsão com todo o conjunto de informações.
autoplot(varejoms) + forecast::autolayer(x.hwm, series = "Holt-Winters mult",
PI = FALSE) + xlab("Ano") + ylab("Consumo - Indice 2011=100") + ggtitle("Forecasts para Índice de volume de vendas no varejo Total
de Mato Grosso do Sul") +
guides(colour = guide_legend(title = "Forecast"))
Outra opção pelo pacote stats e função HoltWinters:
hw <- stats::HoltWinters(varejoms)
predicted <- predict(hw, n.ahead = 24, prediction.interval = TRUE)
dygraph(predicted, main = "Previsão do Varejo de MS") %>% dyAxis("x", drawGrid = TRUE) %>%
dySeries(c("lwr", "fit", "upr"), label = "Varejo") %>% dyOptions(colors = RColorBrewer::brewer.pal(3,
"Set1"))# combine the time series actual and forcasted values
combined <- cbind(predicted, actual = varejoms)
# plot the different values as different series
dygraph(combined, main = "Previsão do Varejo de MS - 24 meses") %>% dyAxis("x",
drawGrid = TRUE) %>% dySeries("actual", label = "Observado") %>% dySeries(paste0("predicted.",
c("lwr", "fit", "upr")), label = "Previsto") %>% dyOptions(colors = RColorBrewer::brewer.pal(3,
"Set1"), drawPoints = TRUE, pointSize = 2) %>% dyRangeSelector() %>% dyHighlight(highlightCircleSize = 4,
highlightSeriesBackgroundAlpha = 0.5, hideOnMouseOut = TRUE) %>% dyEvent("2018-01-01",
"2018", labelLoc = "bottom") %>% dyEvent("2019-01-01", "2019", labelLoc = "bottom") %>%
dyLegend(show = "follow")Referências
BOX, G.E.P.; JENKINS, G.M. Time series analysis: forecasting and control. Revised edition, San Francisco: Holden-Day, 1976.
FERREIRA, Pedro Costa; SPERANZA, Talitha; COSTA, Jonatha (2018). BETS: Brazilian Economic Time Series. R package version 0.4.9. Disponível em: https://CRAN.R-project.org/package=BETS.
GRANGER, C.W.J.; NEWBOLD, P. Forecasting economic time series. Academic Press, 1977.edition, 1987).
HAMILTON, James D. Time Series Analysis. Princeton University Press, 1994.
HENDRY, David F. Dynamic econometrics, Oxford University Press, 1995.
HYNDMAN, Rob. (2018). fpp2: Data for “Forecasting: Principles and Practice” (2nd Edition). R package version 2.3. Disponível em: https://CRAN.R-project.org/package=fpp2.
HYNDMAN, Rob J.; ATHANASOPOULOS, George . Forecasting: principles and practice. Otexts, 2018. Second Edition. Disponível em: https://www.otexts.org/fpp2.
KENNEDY, Peter. A guide to econometrics. 4.ed. Cambridge: MIT Press, 1998. P. 278-279.
MATTOS, Rogério Silva. Decomposição com regressão (Apostila). Juiz de Fora: UFJF, 2016. Disponível em http://www.ufjf.br/rogerio_mattos/files/2009/06/Decomposição-com-Regressão.pdf. Acesso em 22/fev./2018.
MILLS, Terence C. time series techniques for economists. Cambridge University Press, 1990.
MORETTIN, Pedro A.; TOLOI, Clélia M.C. Análise de Séries Temporais. São Paulo: Edgard Blucher/ABE, 2006.
SAX C.; EDDELBUETTEL, D. (2018). Seasonal Adjustment by X-13ARIMA-SEATS in R. Journal of Statistical Software,87(11), 1-17. doi: 10.18637/jss.v087.i11 (URL: https://doi.org/10.18637/jss.v087.i11).
US BUREAU OF THE CENSUS (2013). X-13ARIMA-SEATS Reference Manual Accessible HTML Output Version. Staff Statistical Research Division, US Bureau of the Census, disponível em: http://www.census.gov/ts/x13as/docX13ASHTML.pdf.