Séries Temporais com R: Análise ARIMA do Consumo do Varejo em MS

Abstract

This is an exercise for class use. We analyse data on Retail consumption, from January/2000 to nowadays.

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 com R: Análise do Consumo do Varejo em MS. Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/arima_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é jul/2019 obtida com o pacote BETS e importada do Banco Central do Brasil. Portanto, são 235 observações mensais.

Conforme Hyndman e Athanasopoulos (2018), se combinarmos a diferenciação com autoregressão e um modelo de média móvel, obteremos um modelo ARIMA não sazonal. ARIMA é um acrônimo para Média Móvel Integrada AutoRegressiva (neste contexto, “integração” é o inverso da diferenciação). O modelo completo pode ser escrito como \[ \begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t}, \tag{1.1} \end{equation} \]

onde \(y'_{t}\) é a série diferenciada (pode ter sido diferenciada mais de uma vez). Os “preditores” no lado direito incluem valores defasados de \(y'_{t}\) assim como erros defasados. Chamamos isso de modelo ARIMA(p,d,q), em que p é a ordem da parte autoregressiva (AR(p)); d é a ordem de diferenciação envolvida (I(d)); q é a ordem da parte da média móvel (MA(q)).

2 Dados

Farei de duas formas para o leitor. Uma carrega direto do site do Banco Central do Brasil com o pacote BETS (FERREIRA, SPERANZA e COSTA, 2018) e a outra eu gerei a estrutura idêntica pela função dput() para os leitores que não conseguirem por qualquer motivo o acesso ao site do Banco Central (as vezes vejo isso ocorrer dependendo dos bloqueios da sua rede de internet). A forma pelo dput assume o nome varejoms2 enquanto a extraída pelo BETS tem nome varejoms. Esclareço ao leitor que após baixar a série pelo BETS, fiz o dput e a partir de então, desabilitei o bloco (Chunk) que acessa o BETS apenas para agilizar os cálculos.

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é jul/2019 
# 235 observações mensais
varejoms <- BETSget(1479) 
print(varejoms)
class(varejoms)
dput(varejoms)  # opção para ter os dados como na structure abaixo

suppressMessages(library(readxl))
suppressMessages(library(foreign))
suppressMessages(library(dynlm))
suppressMessages(library(car))
suppressMessages(library(lmtest))
suppressMessages(library(sandwich))
suppressMessages(library(fpp2))
suppressMessages(library(tseries))
suppressMessages(library(zoo))
suppressMessages(library(forecast))
suppressMessages(library(ggplot2))
# 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é jul/2019 Loading the dataset library(readxl) dados <-
# read_excel('dados.xlsx',sheet = 'dados') attach(dados)
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é jul/2019 (em 12.09.2019)

# varejoms <- BETSget(1479)

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, 85.2, 90.9), .Tsp = c(2000, 2019.5,
    12), class = "ts")

A rotina de dados obtidos pelo BETS já retorna a série em formato ts, ou seja, série temporal. Farei então a criação de uma série em diferenças para observar o comportamento da série em nível e em diferenças.

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.

dvarejo <- diff(varejoms)
# estatisticas basicas
summary(varejoms)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  31.60   44.05   64.50   67.28   90.00  129.60 

# Min. 1st Qu.  Median Mean 3rd Qu.  Max.  31.60 44.05 64.50 67.28 90.00 129.60

# plot basico lembrar que em class(), ele já indicou que era ts = serie
# temporal
plot(varejoms)

# pelo pacote dygraph dá mais opções
library(dygraphs)
library(htmlwidgets)
dygraph(varejoms, main = "Índice de volume de vendas no varejo total de Mato Grosso do Sul <br> (Mensal)  (2011=100) BCB 1479") %>%
    dyAxis("x", drawGrid = TRUE) %>%
    dyEvent("2005-1-01", "2005", labelLoc = "bottom") %>%
    dyEvent("2015-1-01", "2015", labelLoc = "bottom") %>%
    dyEvent("2018-1-01", "2018", labelLoc = "bottom") %>%
    dyEvent("2019-1-01", "2019", labelLoc = "bottom") %>%
    dyOptions(drawPoints = TRUE, pointSize = 2)

É possivel 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). Mais a frente, aplicarei o teste de raiz unitária na série para avaliar a estacionariedade de modo mais explícito.

3 Análise da série

Uma ressalva deve ser feita, que no presente exercício, não farei a divisão entre amostra teste e amostra treino, de modo que usarei a série toda para os ajustes. O leitor deve em geral fazer estas divisões para certificar de que o modelo é um bom preditor.

Geralmente, não é possível dizer olhando o gráfico de tempo, quais valores de lags p e q serão apropriados para os dados. Às vezes, é possível usar o gráfico ACF e o gráfico PACF para determinar valores apropriados para p e q. A primeira análise indica olhar os gráficos (correlogramas) da Função de Autocorrelação (FAC, em português, ou ACF em inglês) e Autocorrelação parcial (FACp, em português, ou PACF em inglês).

um gráfico FAC mostra as autocorrelações que medem a relação entre
\(y_t\) e \(y_{t-k}\) para diferentes valores de k. Agora, se \(y_t\) e \(y_{t-1}\) estão correlacionados, então \(y_{t}\) e \(y_{t-2}\) também deve ser correlacionado. Mas isso pode ser porque ambos estão correlacionados com \(y_{t-1}\). Para evitar essas “interferências” é que se faz a FACp. As autocorrelações parciais na FACp medem a relação entre \(y_t\) e \(y_{t-k}\) depois de remover os efeitos dos lags 1, 2, 3, … , (k-1). Portanto, a primeira autocorrelação parcial é idêntica à primeira autocorrelação, porque não há nada entre elas para remover. Cada autocorrelação parcial pode ser estimada como o último coeficiente em um modelo autoregressivo. Isto pode ser obtido com as funções acf e pacf, ou usando as funções ggAcf e ggPacf do pacote forecast. Faremos o mesmo procedimento para a série em nível e para a série em primeira diferença.

3.1 Função de Autocorrelação (FAC) e Autocorrelação parcial (FACp) com defasagem 36

3.1.1 Série em nível

Usarei a rotina do Hyndman e Athanasopoulos (2018).

require(forecast)
varejo <- varejoms
varejo %>%
    ggtsdisplay(main = "")

Neste caso, o ggtsdisplay do pacote forecast já retorna os gráficos da série, e respectivas ACF e PACF.

3.1.2 Série em primeira diferença

Como a série apresenta variações sazonais importantes, assim como tendência importante, claramente não-estacionárias, vou olhar também em primeira diferença.

varejo %>%
    diff() %>%
    ggtsdisplay(main = "Série varejo de MS em primeira diferença")

3.1.3 Primeira diferença sazonal e ACF e PACF

Farei agora a diferença sazonal.

varejo %>%
    diff(lag = 12) %>%
    ggtsdisplay(main = "Série varejo de MS em primeira diferença sazonal")

3.1.4 Ajuste sazonal e ACF e PACF

Aplicaremos o ajuste sazonal tipo STL (Seasonal Decomposition of Time Series by Loess) aos dados.

library(fpp2)
varejo %>%
    stl(s.window = "periodic") %>%
    seasadj() -> varejoadj
autoplot(varejoadj)

E agora os plots de ACF e PACF na série ajustada sazonalmente.

varejoadj %>%
    diff() %>%
    ggtsdisplay(main = "Série varejo de MS em primeira diferença e ajuste sazonal")

Parece indicar para algum AR de ordem 4 na PACF, alguma sazonalidade marcante na ACF (indicando um termo SMA), e indicacao de um ARIMA (4,1,1)(0,1,1)[12].

Uma estimação inicial desse modelo daria algo como

(fit <- Arima(varejoadj, order = c(4, 1, 1), seasonal = c(0, 1, 1)))

Series: varejoadj 
ARIMA(4,1,1)(0,1,1)[12] 

Coefficients:
          ar1      ar2      ar3      ar4     ma1     sma1
      -0.6053  -0.2438  -0.0731  -0.1830  0.1582  -0.4622
s.e.   0.3095   0.1571   0.0871   0.0689  0.3134   0.0621

sigma^2 = 6.179:  log likelihood = -515.75
AIC=1045.51   AICc=1046.03   BIC=1069.33

É possível observar os AICc = 1046.03, suspeitar da não significância de ar3 e olhar os resíduos. A função autoplot retorna o círculo unitário e as raízes unitárias desta estimação, enquanto a checkresiduals fornece os correlogramas residuais.

O círculo unitário indica estabilidade das raízes unitárias do modelo, enquanto o Ljung-Box retorna a não-rejeição de H0 de que os parâmetros de autocorrelação residual são nulos (retorna um grande valor p, sugerindo também que os resíduos são ruído branco). O gráfico ACF dos resíduos do modelo ARIMA (4,1,1)(0,1,1)[12] mostra que todas as autocorrelações estão dentro dos limites à exceção de um lag 23, indicando que os resíduos estão se comportando como ruído branco.

autoplot(fit)

checkresiduals(fit)

    Ljung-Box test

data:  Residuals from ARIMA(4,1,1)(0,1,1)[12]
Q* = 13.198, df = 18, p-value = 0.7797

Model df: 6.   Total lags used: 24

library(FitAR)  # pratico para plotar graficos de Q
# LBQPlot(res, lag.max = 30) res é a série que temos interesse de realizar os
# testes lag.max é o número de lags que se deseja imprimir o gráfico
LBQPlot(residuals(fit), 36)

# Alternativa pelo LSTS
library(LSTS)
p <- Box.Ljung.Test(residuals(fit), lag = 36, main = "LSTS Ljung-Box test")
p + ggthemes::theme_economist_white() + labs(x = "Lag = k", y = "P-valor do Teste Ljung-Box para até lag = k ",
    title = "Série de consumo:", subtitle = "varejo de MS", caption = "Fonte: Elaboração própria a partir de dados do BCB.")

Checando a normalidade dos resíduos.

# checar normalidade dos resíduos
library(tseries)
jarque.bera.test(residuals(fit))

    Jarque Bera Test

data:  residuals(fit)
X-squared = 8.5422, df = 2, p-value = 0.01397

4 ARIMA

Uma estimação automatizada pela função auto.arima {forecast} padrão indica um modelo diferente.

fit2 <- auto.arima(varejo)
summary(fit2)

Series: varejo 
ARIMA(0,1,1)(0,1,1)[12] 

Coefficients:
          ma1     sma1
      -0.4558  -0.4691
s.e.   0.0592   0.0617

sigma^2 = 6.312:  log likelihood = -520.11
AIC=1046.23   AICc=1046.34   BIC=1056.44

Training set error measures:
                     ME     RMSE     MAE        MPE    MAPE      MASE
Training set 0.01107351 2.430819 1.79931 0.02282225 2.64575 0.3926471
                    ACF1
Training set -0.01460341

# p-values dos coeficientes
pvalues <- (1 - pnorm(abs(fit2$coef)/sqrt(diag(fit2$var.coef)))) * 2
pvalues

         ma1         sma1 
1.332268e-14 2.775558e-14 

library(FitAR)  # pratico para plotar graficos de Q
LBQPlot(residuals(fit2), 36)

# Alternativa pelo LSTS
library(LSTS)
p2 <- Box.Ljung.Test(residuals(fit2), lag = 36, main = "LSTS Ljung-Box test")
p2 + ggthemes::theme_economist_white() + labs(x = "Lag = k", y = "P-valor do Teste Ljung-Box para até lag = k ",
    title = "Série de consumo:", subtitle = "varejo de MS", caption = "Fonte: Elaboração própria a partir de dados do BCB.")

# checar normalidade dos resíduos
library(tseries)
jarque.bera.test(residuals(fit2))

    Jarque Bera Test

data:  residuals(fit2)
X-squared = 8.847, df = 2, p-value = 0.01199

Neste caso, o modelo teve a indicação de ARIMA(0,1,1)(0,1,1)[12], diferente do anterior (ARIMA(4,1,1)(0,1,1)[12] - definido pela intuição dos correlogramas), como reflexo do default da função admitir lags máximos e procedimentos que restringem a otimização da função verossimilhança. O AICc foi de 1046.34, pouco acima do anterior, indicando ajuste inferior que o primeiro.

5 Modelo ARIMA automático mais geral

A recomendação do Hyndman e Athanasopoulos (2018), nestas situações, é fazer um modelo auto.arima mais geral, que investigue mais possibilidades de modelos,fazendo os argumentos da auto.arima da forma stepwise=FALSE e approximation=FALSE. Se estiver usando a série com ajuste sazonal (dessazonalizada), será interessante usar tambem o argumento seasonal=FALSE.

require(fpp2)
fit3 <- auto.arima(varejoms, stepwise = FALSE, approximation = FALSE)
summary(fit3)

Series: varejoms 
ARIMA(4,1,0)(0,1,1)[12] 

Coefficients:
          ar1      ar2      ar3      ar4     sma1
      -0.4529  -0.1741  -0.0504  -0.1874  -0.4610
s.e.   0.0662   0.0726   0.0725   0.0660   0.0618

sigma^2 = 6.158:  log likelihood = -515.88
AIC=1043.75   AICc=1044.14   BIC=1064.17

Training set error measures:
                      ME     RMSE      MAE        MPE     MAPE      MASE
Training set 0.009480306 2.384598 1.779725 0.02458176 2.619097 0.3883732
                    ACF1
Training set 0.004907171

# p-values dos coeficientes
pvalues <- (1 - pnorm(abs(fit3$coef)/sqrt(diag(fit3$var.coef)))) * 2
pvalues

         ar1          ar2          ar3          ar4         sma1 
7.652989e-12 1.642129e-02 4.869277e-01 4.524327e-03 8.726353e-14 

Neste cenário, a rotina auto.arima encontra o ARIMA(4,1,0)(0,1,1)[12], com AICc de 1044.14, portanto, melhor (menor) que as duas primeiras estimações. como intuido anteriormente, ele detecta o AR(4), a Integração (I(1)), a integração sazonal e o SMA(1). O termo ar3 não foi significativo, mas os demais foram.

5.0.1 Checando resíduos

Os modelos podem ser considerados white noise, mas não normais.

checkresiduals(fit3)

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0)(0,1,1)[12]
Q* = 13.381, df = 19, p-value = 0.8185

Model df: 5.   Total lags used: 24

require(FitAR)  # pratico para plotar graficos de Q
LBQPlot(residuals(fit3), 36)

# Alternativa pelo pacote LSTS
library(LSTS)
p3 <- Box.Ljung.Test(residuals(fit3), lag = 36, main = "LSTS Ljung-Box test")
p3 + ggthemes::theme_economist_white() + labs(x = "Lag = k", y = "P-valor do Teste Ljung-Box para até lag = k ",
    title = "Série de consumo:", subtitle = "varejo de MS", caption = "Fonte: Elaboração própria a partir de dados do BCB.")

# checar normalidade dos resíduos
require(tseries)
jarque.bera.test(residuals(fit3))

    Jarque Bera Test

data:  residuals(fit3)
X-squared = 8.0982, df = 2, p-value = 0.01744

Pode-se fazer o gráfico dos forecasts deste modelo:

autoplot(forecast(fit3, h = 24), title = "Forecasts de volume de vendas no varejo total de Mato Grosso do Sul",
    xlab = "Ano", ylab = "Índice (2011=100)")

5.0.2 Teste de raiz unitária

O teste de raiz unitária deveria ser feito antes de rodar o ARIMA, a fim de investigar e/ou confirmar a não-estacionariedade da série, que foi intuida a partir da ACF e PACF. A rotina auto.arima já faz esse teste em seu processo. Existem diferentes testes de raiz unitária possíveis (ADF, PP, KPSS) assim como diferentes testes para a sazonalidade.

5.0.2.1 Teste ADF

O teste de Dickey-Fuller aumentado (ADF) é nossa primeira escolha. O padrão é que H0: série não estacionária ou série tem raiz unitária.

ADF.test <- adf.test(varejoms)
ADF.test

    Augmented Dickey-Fuller Test

data:  varejoms
Dickey-Fuller = -1.7075, Lag order = 6, p-value = 0.6988
alternative hypothesis: stationary

5.0.2.2 Teste ADF invertendo a hipótese nula

Neste caso, o padrão é que H0: a série é estacionária ou a série não tem raiz unitária. Neste caso, a hipótese alternativa é que a série é “explosiva”.

ADF.test12 <- adf.test(varejoms, k = 12, alternative = c("explosive"))
ADF.test12

    Augmented Dickey-Fuller Test

data:  varejoms
Dickey-Fuller = -1.0897, Lag order = 12, p-value = 0.07762
alternative hypothesis: explosive

6 ARIMA - Codigo do Hyndman e Athanasopoulos, secao 8.1

A função ndiffs usa um teste de raiz unitária para determinar o número de diferenças necessárias para que séries temporais xsejam feitas estacionárias.

Se test="kpss", o teste KPSS é usado com a hipótese nula que x possui uma raiz estacionária contra uma alternativa de raiz unitária. Em seguida, o teste retorna o menor número de diferenças necessárias para passar no teste ao nível alpha.

Se test="adf", o teste ADF é utilizado e se test="pp", o teste de Phillips-Perron é usado. Em ambos os casos, a hipótese nula é que x tem uma raiz unitária contra uma alternativa de série estacionária. Em seguida, o teste retorna o menor número de diferenças necessárias para falhar o teste no nível alpha.

A função nsdiffs usa testes de raiz unitária sazonal para determinar o número de diferenças sazonais necessárias para que as séries temporais x sejam feitas estacionárias (possivelmente com algum diferencial de atraso também). Se test="ch", o teste de Canova-Hansen (1995) é usado (com hipótese nula de sazonalidade determinística) e se test="ocsb" o teste de Osborn-Chui-Smith-Birchenhall (1988) é usado (com hipótese nula de existir uma raiz unitária sazonal).

x <- varejoms
ns <- nsdiffs(x)  #numero de diferencas sazonais
ns

[1] 1

if (ns > 0) {
    xstar <- diff(x, lag = frequency(x), differences = ns)
} else {
    xstar <- x
}

Portanto, temos uma diferença sazonal necessária (ns=1) e a série na diferença sazonal está armazenada em xstar.

nd <- ndiffs(xstar)  # numero de diferencas
nd

[1] 1

if (nd > 0) {
    xstar <- diff(xstar, differences = nd)
}
# xstar será a série devidamente em diferenças sazonais e diferenças normais
autoplot(xstar)

auto.arima(xstar)  # ver que o resultado é igual, mas agora nao tem diferencas

Series: xstar 
ARIMA(0,0,1)(0,0,1)[12] with zero mean 

Coefficients:
          ma1     sma1
      -0.4558  -0.4691
s.e.   0.0592   0.0617

sigma^2 = 6.312:  log likelihood = -520.11
AIC=1046.23   AICc=1046.34   BIC=1056.44

modelo <- auto.arima(xstar, stepwise = FALSE, approximation = FALSE)
summary(modelo)

Series: xstar 
ARIMA(4,0,0)(0,0,1)[12] with zero mean 

Coefficients:
          ar1      ar2      ar3      ar4     sma1
      -0.4529  -0.1741  -0.0504  -0.1874  -0.4610
s.e.   0.0662   0.0726   0.0725   0.0660   0.0618

sigma^2 = 6.158:  log likelihood = -515.88
AIC=1043.75   AICc=1044.14   BIC=1064.17

Training set error measures:
                     ME     RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.01044227 2.453404 1.882938 -Inf  Inf 0.5055836 0.004680592

autoplot(modelo)

Comentário: uma vez que a série utilizada em auto.arima é a resultante xstar, o mmodelo não identifica a necessidade de diferenciação nem diferenciação sazonal. Com esse procedimento, ele encontra o ARIMA(4,0,0)(0,01)[12] para xstar que já teve as diferenças aplicadas sobre varejoms. Da mesma forma que em fit e fit3, o modelo com argumentos stepwise=FALSE, approximation=FALSE retorna os melhores ajustes.

Referências

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.

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, R.J., & ATHANASOPOULOS, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. Disponível em: https://otexts.com/fpp2/. Accessed on 12 Set 2019.