Previsão de Séries Temporais com o R

Este arquivo foi utilizado como arquivo de apoio à aulas ministradas sobre previsão de demanda utilizando o R. Para um tramento mais aprofundado sobre o tema, recorrer ao excelente livro Forecating: Principles and Practice.

Este post possui uma apresentação relacionada neste link.

Sobre o R

  • O R é um ambiente de computação estatística, com mais de 13 mil pacotes publicados.
  • Cada um destes pacotes tem um fim específico. Neste curso, utilizaremos principalmente as bibliotecas ‘forecast’ e o fpp2. Ambos os pacotes são utlizados no livro Forecating: Principles and Practice.

Passos para Realizar Previsões

  1. Faça um gráfico dos dados;
  2. Selecione uma função para a previsão;
  3. Estime os parâmetros da função;
  4. Avalie a Qualidade do modelo de previsão;
  5. Selecione e implemente o melhor modelo encontrado.

Instalando Bibliotecas no R

  • Instale E carregue as bibliotecas que iremos utilizar nesta aula rodando estes comandos:
bibliotecas = c("forecast", "fpp2", "readxl")

install.packages(bibliotecas)

# Caso encontre algum erro, rode o install packages separadamente para cada biblioteca:
install.packages("forecast")
install.packages("fpp2")
install.packages("readxl")

Antes de começar, carregue as bibliotecas:

library(forecast)
library(fpp2)
library(readxl)

Trabalhando com a Primeira Série

  • Observe as série gold:
# A série temporal "gold" foi carregada pela biblioteca forecast.
autoplot(gold)

Padrões em séries temporais

  • Tendência: Os dados possuem uma tendência geral de aumento ou queda (exemplo: Crescimento da população).
  • Sazonalidade: Os dados se comportam de modo similar, obedecendo padrões com durações fixas (exemplo: venda de ovos de páscoa).
  • Ciclicidade: A série possui um comportamento variando em ciclos, porém sem duração fixa (exemplo: ciclos econômicos.).

Observando outras series temporais

# Produção de Lâ na Austrália
autoplot(woolyrnq) 

Observando outras series temporais

# Produção de Gás na Austrália
autoplot(gas) 

Observando o a Sazonalidade das Séries

# A função frequency determina a frequência da série. 
# Para dados sazonais, irá definir o período dominante da sazonalidade, 
# e para dados em ciclos, a duração média dos ciclos.
frequency(gas) 
## [1] 12

Observando Gráficos Sazonais

Observando a série temporal:

# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda Anti-Diabeticos")

Observando um gráfico sazonal:

ggseasonplot(a10) 

Observando um gráfico sazonal “polar”:

ggseasonplot(a10, polar = T)

Observando a Venda de Cerveja:

beer <- window(ausbeer, start = 1992)
autoplot(beer) + ylab("Venda de Cerveja")

Observando a Venda de Cerveja e sua Sazonalidade:

ggseasonplot(beer) 

Mais uma forma de visualizar a demanda sazonal:

ggsubseriesplot(beer)

Criando Objetos ‘ts’ para Séries Temporais a partir de dados reais

Para manipular séries temporais no R, é interessante criar objetos do tipo ‘ts’. Podemos criar objetos a partir de arquivos do excel, csv, bases de dados, ou mesmo APIs públicas.

Podemos utilizar a biblioteca ‘readxl’ para ler arquivos do excel. Obtenha aqui o arquivo de dados VendasCarros.xlsx.

library(readxl)
dados_excel = readxl::read_xlsx(path = "2018-10-10-previsao-serie-temporais/VendasCarros.xlsx", sheet = "Dados")
head(dados_excel)
## # A tibble: 6 x 4
##     Ano   Mês VendasCarros VendasCarrosAux
##   <dbl> <dbl>        <dbl>           <dbl>
## 1  2018     1         2493            2486
## 2  2018     2         2875            3003
## 3  2018     3         3504            3454
## 4  2018     4         3786            3665
## 5  2018     5         4094            4353
## 6  2018     6         4994            4842

Transformando um Data Frame em uma Série Temporal

  • Para trabalhar com séries temporais, iremos transformar esta tabela:
ts_vendas = ts(data = dados_excel$VendasCarros, start = c(2018,1),frequency = 12)
ts_vendas
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2018 2493 2875 3504 3786 4094 4994 4910 5575 5839 6202 6833 7382
## 2019 2470 3059 3537 4003 4110 4516 5364 5854 5752 6089 7025 7801
## 2020 2473 2884 3328 3689 4133 4932 5097 5397 6205 6374 6697 8024
  • Agora que temos uma série temporal, vamos plotar:
forecast::autoplot(ts_vendas)

Previsões com Modelos de Suavização Exponencial

Suavização Exponencial Simples

  • Método da média simples: Leva em consideração todos os períodos;
  • Método “Naive”: Leva em consideração apenas o último período;
  • Método da Suavização Exponencial: Leva em consideração a demanda observada nos períodos anteriores, porém faz com que o impacto dos períodos anteriores caia progressivamente.

Notação da Previsão: \(\hat{y}_{t+h|t}\)

Previsão para a demanda \(y_{t+h}\), considerando dados disponiveis até \(y_{t}\)

Equação da Previsão: \[\hat{y}_{t+h|t} = \alpha y_{t} + \alpha(1-\alpha)y_{t-1} + \alpha(1-\alpha)^2y_{t-2} + ... \ para \ 0 \leq \alpha \leq 1\]

Modelo de Suavização Exponencial Simples

Previsão: \[\hat{y}_{t+h|t} = l_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)l_{t-1}\] O valor \(\alpha\) é obtido minimizando os erros ao rodar o modelo em um set de teste.

Estimando o Modelo com o R

  • A função ‘ses’ (Simple Exponetial Smoothing) estima um modelo de suavização exponencial simples:
dados_petroleo = window(oil, start= 1996)
previsao_petroleo = ses(dados_petroleo, h = 5) # Previsao para os próximos 5 anos
summary(previsao_petroleo)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = dados_petroleo, h = 5) 
## 
##   Smoothing parameters:
##     alpha = 0.8339 
## 
##   Initial states:
##     l = 446.5868 
## 
##   sigma:  29.8282
## 
##      AIC     AICc      BIC 
## 178.1430 179.8573 180.8141 
## 
## Error measures:
##                    ME     RMSE     MAE      MPE     MAPE      MASE
## Training set 6.401975 28.12234 22.2587 1.097574 4.610635 0.9256774
##                     ACF1
## Training set -0.03377748
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2014       542.6806 504.4541 580.9070 484.2183 601.1429
## 2015       542.6806 492.9073 592.4539 466.5589 618.8023
## 2016       542.6806 483.5747 601.7864 452.2860 633.0752
## 2017       542.6806 475.5269 609.8343 439.9778 645.3834
## 2018       542.6806 468.3452 617.0159 428.9945 656.3667
  • Visualizando a Previsão:
autoplot(previsao_petroleo)

Modelo de Suavização Exponencial com Tendência Linear - Holt

Previsão: \[\hat{y}_{t+h|t} = l_t + h b_t\] Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+b_{t-1})\] Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta)b_{t-1}\]

Estimando o Modelo com o R

  • A função ‘holt’ (Simple Exponetial Smoothing) estima o modelo de Holt (que considera a tendência):
previsao_petroleo_holt = holt(dados_petroleo, h = 5) # Previsao para os próximos 5 anos
autoplot(previsao_petroleo_holt)

  • Também podemos observar os parâmetros estimados para o modelo:
summary(previsao_petroleo_holt)
## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = dados_petroleo, h = 5) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 428.4833 
##     b = 5.5771 
## 
##   sigma:  27.8684
## 
##      AIC     AICc      BIC 
## 177.2927 182.2927 181.7446 
## 
## Error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.2851768 24.57757 20.53231 -0.3285466 4.337459 0.8538816
##                   ACF1
## Training set 0.3429178
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2014       534.4382 498.7234 570.1529 479.8172 589.0591
## 2015       540.0148 504.3000 575.7295 485.3938 594.6357
## 2016       545.5913 509.8766 581.3061 490.9704 600.2123
## 2017       551.1679 515.4532 586.8827 496.5469 605.7889
## 2018       556.7445 521.0298 592.4592 502.1235 611.3655

Modelo de Holt com “Damp”

  • Podemos utilizar o modelo Holt com “Damp”, pressupondo que o crescimento não será linear no longo prazo.

Previsão: \[\hat{y}_{t+h|t} = l_t + (\phi+\phi^2+...+\phi^h) b_t\]

Nível: \[l_t = \alpha y_t + (1-\alpha)(l_{t-1}+\phi b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) \phi b_{t-1}\] O parâmetro \(\phi\) é entre 0 e 1. Se o parâmetro for igual a 1, o crescimento será linear.

Utilizando o Modelo de Holt com “Damp”:

previsao_petroleo_holt = holt(dados_petroleo, h = 15, PI = F)
previsao_petroleo_holt_damp = holt(dados_petroleo, h = 15, damped = T, PI = F) # Previsao para os próximos 5 anos
autoplot(dados_petroleo) + 
  autolayer(previsao_petroleo_holt, series="Holt Linear") + 
  autolayer(previsao_petroleo_holt_damp, series="Holt com Damp")

Modelo de Suavização Exponencial com Tendência e Sazonalidade - Holt-Winters

  • Método Aditivo: Adequado quando a amplitude dos ciclos de sazonalidade não está correlacionada ao tempo.
  • Método Multiplicativo: Adequado quando a amplitude dos ciclos de sazonalidade está correlacionadao ao tempo.

Modelo Holt-Winters Aditivo

Previsão: \[\hat{y}_{t+h|t} = l_t + hb_t + s_{t-m+h_m^+}\]

Nível: \[l_t = \alpha (y_t-s_{t-m}) + (1-\alpha)(l_{t-1}+b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]

Sazonalidade: \[b_t = \gamma(y_t - l_{t-1 - b_{t-1}}) + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a zero.

Modelo Holt-Winters Multiplicatio

Previsão: \[\hat{y}_{t+h|t} = (l_t + hb_t) s_{t-m+h_m^+}\]

Nível: \[l_t = \alpha \frac{y_t}{s_{t-m}} + (1-\alpha)(l_{t-1}+b_{t-1})\]

Tendência: \[b_t = \beta *(l_t-l_{t-1}) + (1-\beta) b_{t-1}\]

Sazonalidade: \[b_t = \gamma \frac{y_t}{l_{t-1} + b_{t-1}} + (1-\gamma)s_{t-m} \] \(s_{t-m+h_m^+}\): componente de sazonalidade do último ano de dados disponíveis. \(m\): Período de sazonalidade. A média do componente de sazonalidade tende a 1.

Exemplo - Holt-Winters Aditivo

Relembrando a série de anti-glicêmicos.

# Produção de Gás na Austrália
autoplot(a10) + ylab("Demanda AntiDiabeticos")

Realizando a Previsão:

# Produção de Gás na Austrália
previsao_anti_diab_aditivo = hw(a10, seasonal = "additive", PI= F)
previsao_anti_diab_multiplicativo = hw(a10, seasonal = "multiplicative", PI = F)
autoplot(a10) + ylab("Demanda de Remédios") +
  autolayer(previsao_anti_diab_aditivo, series="HW Add.") + 
  autolayer(previsao_anti_diab_multiplicativo, series="HW Mult.")

Desafio:

  • Colete dados de demanda de diferentes produts em sua empresa;
  • Verifique quais são as técnicas de previsão aplicadas a estes produtos;
  • Identifique se há ou não tendência e sazonalidade;
  • Aplique os métodos aprendidos;
  • Proponha um modelo de previsão novo;
  • Compare a acurácia do novo modelo em relação ao que existe atualmente na empresa.

ARIMA

O modelo ARIMA não será coberto nesta aula, porém a biblioteca forecast possui uma função para estimar modelos Arima de modo automático.

modelo_arima = forecast::auto.arima(a10)

autoplot(forecast(modelo_arima))

Mais Aspectos Técnicos

O modelo “Naive”

  • O modelo “naive” simplesmente pressupõe que o futuro repetirá o passado, logo a previsão é correspondente ao último valor observado.
previsao = naive(oil)
autoplot(oil, series="Dados") + 
  xlab("Ano") +
  autolayer(fitted(previsao), series = "Previsao") + 
  ggtitle("Produção de Petróleo na Arábia Saudita")

Observando os Erros

  • Os erros também podem ser plotados. Espera-se que os erros sejam “white noise”. Se isso é verdade, é provável que o modelo tenha captado toda a informação disponível nos dados.
autoplot(residuals(previsao))

Pressupostos sobre os Erros

  • Os erros deveriam ser não-correlacionados;
  • Os erros devem ter média zero; Propriedades úteis:
  • Variância Constante;
  • São normalmente distribuídos.

Verificando Pressupostos sobre os Erros

  • A função ‘checkresiduals’ verifica os pressupostos indicados anteriormente. Espera-se que o resultado do p-valor do Ljung-box seja acima de 0.05, indicando que os erros não são correlacionados.
checkresiduals(previsao)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 11.814, df = 9.8, p-value = 0.2824
## 
## Model df: 0.   Total lags used: 9.8

O que fazer se os erros não forem normais?

  • Ainda assim, as previsões podem ser boas e podem ser utilizadas, porém os intervalos de predição podem ser muito justos ou amplos em função deste problema.

Training and Test Sets

  • Traning Set: Parte dos dados que você utiliza para construir o modelo.
  • Test Set: Parte dos dados que você utiliza para testar o modelo.
  • O test set não pode ser usado para calcular a previsão.
  • Um modelo que se ajusta bem aos dados de treinamento não necessáriamente terá uma boa previsão;
  • É comum construirmos um modelo altamente complexo que possui poucos erros nos dados de treinamento, porém gera péssimas previsões. Isto é chamado de overfitting, e deve ser evitado.

Medidas de Acurácia da Previsão

  • Mean Absolute Error: \(MAE = avg(|e_t|)\)
  • Mean Square Error: \(MAE = avg(|e_t^2|)\)
  • Mean Absolute Percentage Error: \(MAPE = 100 * avg(|e_t/y_t|)\)
  • Mean Absolute Scaled Error: \(MASE = MAE / Q\)

Exemplo com um Modelo Naive

  • Separamos o modelo em duas partes, e geramos um modelo naive para ilustrar este ponto:
treinamento = window(oil, end=2003)
teste = window(oil, start= 2004)
previsao = naive(treinamento,h = 10) # h = Número de períodos a prever
autoplot(previsao) + 
  ylab("Vendade Petróleo") + 
  autolayer(teste, series = "Dados de Teste")

O Comando Accuracy

  • Usar o Comando Accuracy para observar o modelo:
accuracy(previsao, teste)
##                    ME     RMSE      MAE      MPE      MAPE      MASE
## Training set  9.87358 52.56156 39.42504 2.506565 12.570647 1.0000000
## Test set     21.60250 35.09832 29.97666 3.963914  5.777875 0.7603458
##                   ACF1 Theil's U
## Training set 0.1801528        NA
## Test set     0.4029519  1.184862

Obtendo Mais Séries de Dados reais:

  • O agregador de dados Quandl possui milhares de séries temporais disponíveis diretamente no R, pela biblioteca Quandl.
  • Exemplo: obtendo a série do índice Bovespa
library(Quandl)
ts_bovespa = Quandl("BCB/7", type = "ts")
autoplot(ts_bovespa)