Séries Temporais com R: Análise do Consumo do Varejo em MS com fable e feasts

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 com fable e feasts. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/fable_feasts_varejoms.

QRCode

library(qrcode)
plot(qr_code("https://rpubs.com/amrofi/fable_feasts_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é jan/2020 obtida com o pacote BETS e importada do Banco Central do Brasil. Portanto, são 241 observações mensais.

Realizaremos os modelos ETS (Error Trend Seasonal) conforme Hyndman e Athanasopoulos (2021,2018), ou seja, com o pacote forecast e com o formato tsibble pelo fable e feasts.

Sugiro olhar meus videos no Youtube em: Aplicação de ETS com ações: https://youtu.be/s5S43QBllR0; Video teórico de ETS: https://youtu.be/ej8hH6g9j3Y.

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). 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é jan/2020  
# 241 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é jan/2020

library(BETS)
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, 91, 94.4, 93.4, 95.9, 102.6,
    116, 94.9), .Tsp = c(2000, 2020, 12), class = "ts")

A rotina de dados obtidos pelo BETS já retorna a série em formato ts, ou seja, série temporal. Inicialmente olharei as estatísticas descritivas da série. Em seguida farei um plot básico da série, ú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.30   67.10   68.08   90.80  129.60 

class(varejoms)

[1] "ts"

varejoms

       Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
2000  35.2  35.6  39.2  40.5  41.6  40.4  40.8  38.7  37.3  37.6  35.6  47.4
2001  34.2  32.2  38.0  37.5  38.8  35.0  38.4  40.2  38.2  39.3  36.0  46.3
2002  36.4  34.0  39.0  37.8  38.9  35.3  37.2  38.1  35.6  38.3  35.6  45.7
2003  32.2  31.6  35.2  36.8  37.5  34.8  38.4  38.1  37.0  39.0  37.3  49.0
2004  35.8  35.3  40.2  41.3  43.9  41.6  45.9  42.3  42.2  44.0  41.3  56.9
2005  38.5  38.5  45.3  43.6  46.2  44.3  47.5  46.5  46.4  46.0  44.1  61.0
2006  42.0  40.2  44.6  44.5  47.8  45.3  46.5  48.5  47.7  50.2  49.3  64.5
2007  47.0  46.8  51.0  50.5  55.0  51.3  52.8  55.3  54.8  55.6  55.3  72.2
 [ reached getOption("max.print") -- omitted 13 rows ]

Vê-se que não existem valores negativos e eles variam entre o mínimo de 31.60 e o máximo de 129.60. Graficamente, tem-se:

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

Ou pelo pacote dygraphs:

# pelo pacote dygraph dá mais opções
library(dygraphs)
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 no plot 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

Farei a divisão entre amostra teste e amostra treino para controle de acurácia, pelo tradicional 80-20 (241 observações x 80% = 193 para treino), portanto a amostra treino será de jan-2000 a jan-2016. O leitor deve em geral fazer estas divisões para certificar de que o modelo é um bom preditor. Uma análise mais detalhada da série pode ser feita fazendo o chunk abaixo. Não estaremos realizando o modelo ARIMA neste exercício, mas fica a dica para o leitor.

varejo <- varejoms
varejo %>%
    forecast::ggtsdisplay(main = "Série de Varejo MS")

library(fpp2)
train <- window(varejoms, end = c(2016, 1))  # fim das 193 observações
h <- length(varejoms) - length(train)
h  # horizonte de teste dentro da amostra observada

[1] 48

A estimação automática pelo pacote forecast será quando não especificarmos o modelo, equivalente a fazer um ETS(ZZZ), quando então ele testará todas as possibilidades entre aditivo, multiplicativo ou nenhum efeito:

ETS <- forecast::forecast(ets(train), h = h)
summary(ETS)

Forecast method: ETS(M,Ad,M)

Model Information:
ETS(M,Ad,M) 

Call:
 ets(y = train) 

  Smoothing parameters:
    alpha = 0.529 
    beta  = 0.0222 
    gamma = 1e-04 
    phi   = 0.98 

  Initial states:
    l = 39.8072 
    b = 0.0336 
    s = 1.2715 0.9724 1.0167 0.9719 1.004 1.0114
           0.957 1.0266 0.968 0.9733 0.8827 0.9445

  sigma:  0.0343

     AIC     AICc      BIC 
1297.307 1301.238 1356.035 

Error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE
Training set 0.1785147 2.112127 1.593527 0.2726068 2.624009 0.3337518
                   ACF1
Training set 0.07061013

Forecasts:
         Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
Feb 2016       86.21033  82.41812  90.00255  80.41064  92.01002
Mar 2016       95.09920  90.32250  99.87590  87.79386 102.40453
Apr 2016       94.61642  89.29945  99.93338  86.48482 102.74802
May 2016      100.37904  94.15827 106.59981  90.86519 109.89289
Jun 2016       93.61309  87.28337  99.94282  83.93262 103.29357
Jul 2016       98.96709  91.72670 106.20748  87.89387 110.04031
Aug 2016       98.27817  90.55068 106.00566  86.45999 110.09635
Sep 2016       95.17126  87.17298 103.16954  82.93894 107.40357
Oct 2016       99.59857  90.69357 108.50356  85.97956 113.21757
Nov 2016       95.29284  86.26461 104.32107  81.48535 109.10032
Dec 2016      124.64406 112.17374 137.11438 105.57235 143.71577
Jan 2017       92.61328  82.85825 102.36831  77.69425 107.53231
Feb 2017       86.58377  77.00783  96.15972  71.93863 101.22892
Mar 2017       95.50274  84.43897 106.56651  78.58216 112.42331
Apr 2017       95.00972  83.50544 106.51401  77.41544 112.60401
May 2017      100.78780  88.05705 113.51854  81.31780 120.25779
Jun 2017       93.98653  81.62447 106.34859  75.08039 112.89267
Jul 2017       99.35384  85.76824 112.93945  78.57645 120.13124
Aug 2017       98.65441  84.65124 112.65758  77.23841 120.07042
Sep 2017       95.52819  81.47248 109.58390  74.03183 117.02454
 [ reached 'max' / getOption("max.print") -- omitted 28 rows ]

Portanto, o modelo escolhido acima foi um ETS(M,Ad,M) cujos parâmetros de suavização foram:

• alpha = 0.529

• beta = 0.0222

• gamma = 1e-04

• phi = 0.98

Ou seja, os parâmetros para o nível (alpha), tendência (beta), sazonalidade (gamma = 0.0001) e amortecimento (phi). Entretanto, ressaltamos que se trata de uma série com volatilidade e não estacionária. Recomendamos realizar uma transformação Box-Cox previamente à realização do ETS. Deixo o modelo acima em standby portanto.

3.1 Transformação da série

A transformação Box-Cox pode ser realizada diretamente pelo comando BoxCox.lambda:

(BoxCox.lambda(varejoms))

[1] 0.05583817

Portanto, rodaremos novo ETS com essa transformação, que praticamente reproduz a logaritmica:

ETS <- ets(train, lambda = 0.05583817)
summary(ETS)

ETS(A,A,A) 

Call:
 ets(y = train, lambda = 0.05583817) 

  Box-Cox transformation: lambda= 0.0558 

  Smoothing parameters:
    alpha = 0.5126 
    beta  = 0.0173 
    gamma = 1e-04 

  Initial states:
    l = 4.0782 
    b = -0.0028 
    s = 0.3035 -0.0249 0.021 -0.0268 0.0068 0.0123
           -0.0598 0.0376 -0.034 -0.025 -0.1421 -0.0688

  sigma:  0.0424

      AIC      AICc       BIC 
-187.2360 -183.7389 -131.7703 

Training set error measures:
                     ME     RMSE      MAE       MPE     MAPE     MASE      ACF1
Training set 0.02273207 2.086761 1.560396 0.1122085 2.572966 0.326813 0.1377376

Agora temos o resultado: ETS(A,A,A),

• alpha = 0.5126
• beta = 0.0173
• gamma = 1e-04

O plot do forecast será:

autoplot(forecast.ets(ETS, h = h)) + xlab("Ano") + ylab("Índice 2011=100") + ggtitle("Índice de volume de vendas no varejo total de Mato Grosso do Sul <br> (Mensal)  (2011=100) BCB 1479")

A investigação da acurácia será do forecast estimado com a série treino:

ETS.fc <- forecast::forecast(ETS, h = h)
forecast::accuracy(ETS.fc, varejoms)["Test set", "RMSE"]  # apenas RMSE

[1] 13.85666

forecast::accuracy(ETS.fc, varejoms)  # todas as medidas

                       ME      RMSE       MAE         MPE      MAPE     MASE
Training set   0.02273207  2.086761  1.560396   0.1122085  2.572966 0.326813
Test set     -12.85049406 13.856660 12.850494 -13.8986653 13.898665 2.691437
                  ACF1 Theil's U
Training set 0.1377376        NA
Test set     0.3017299  1.546331

Ou seja, o modelo se comporta melhor com o treino que com o teste, isso porque a série exibe uma quebra importante após 2015, o período de recessão, diferente do forte crescimento antes de 2015.

A seguir, faremos o mesmo exercicio com fable e feasts.

4 Estimação com fable e feasts

Primeiro faremos a conversão de objeto ts para tsibble, ou time series table. Converterei a série inteira varejoms e depois farei o split entre train e test.

library(fable)
library(feasts)
library(fabletools)
library(tsibble)
library(tsibbledata)
library(lubridate)
library(dplyr)
varejo.ts <- as_tsibble(varejoms)

Faremos o split com uso da função filter_index e especificando filter_index(~ "2016-01"):

train <- varejo.ts %>%
    filter_index(~"2016-01")
train

Podemos plotar e estimar. Já sabemos que precisa de um lambda = 0.05583817. A função report fornecerá a saída do modelo.

fit <- train %>%
    model(ets = ETS(box_cox(value, 0.05583817)))
fit

# A mable: 1 x 1 ets <model> 1 <ETS(A,A,A)>
report(fit)

Series: value 
Model: ETS(A,A,A) 
Transformation: box_cox(value, 0.05583817) 
  Smoothing parameters:
    alpha = 0.5126422 
    beta  = 0.01728316 
    gamma = 0.000137223 

  Initial states:
     l[0]         b[0]      s[0]       s[-1]      s[-2]      s[-3]       s[-4]
 4.078226 -0.002772345 0.3034788 -0.02490655 0.02104653 -0.0267539 0.006841962
      s[-5]       s[-6]     s[-7]       s[-8]       s[-9]     s[-10]
 0.01227541 -0.05981128 0.0376117 -0.03398304 -0.02495835 -0.1420682
      s[-11]
 -0.06877305

  sigma^2:  0.0018

      AIC      AICc       BIC 
-187.2360 -183.7389 -131.7703 

Agora o forecast à frente 60 meses. Coloquei 48 meses do treino mais 12 meses fora da série observada.

fc <- fit %>%
    forecast(h = "5 years")
fc

E a acurácia da parte de teste foi:

fit %>%
    forecast(h = "4 years", bootstrap = FALSE) %>%
    accuracy(varejo.ts, measures = point_accuracy_measures)

Segue o plot da estimação:

fc %>%
    autoplot(varejo.ts)

E as nossas estimativas para até Jan/2021 são:

knitr::kable(fc)

.model	index	value	.mean
ets	2016 fev	t(N(5.1, 0.0018))	87.84039
ets	2016 mar	t(N(5.2, 0.0023))	96.50878
ets	2016 abr	t(N(5.2, 0.0028))	96.13526
ets	2016 mai	t(N(5.3, 0.0034))	101.92931
ets	2016 jun	t(N(5.2, 0.004))	94.82510
ets	2016 jul	t(N(5.3, 0.0047))	100.58323
ets	2016 ago	t(N(5.3, 0.0053))	100.47879
ets	2016 set	t(N(5.2, 0.0061))	98.21280
ets	2016 out	t(N(5.3, 0.0068))	102.23612
ets	2016 nov	t(N(5.2, 0.0076))	98.98283
ets	2016 dez	t(N(5.6, 0.0085))	127.78863
ets	2017 jan	t(N(5.2, 0.0094))	96.29544
ets	2017 fev	t(N(5.1, 0.01))	91.26773
ets	2017 mar	t(N(5.3, 0.011))	100.26624
ets	2017 abr	t(N(5.2, 0.012))	99.89263
ets	2017 mai	t(N(5.3, 0.013))	105.91364
ets	2017 jun	t(N(5.2, 0.014))	98.56220
ets	2017 jul	t(N(5.3, 0.016))	104.54821
ets	2017 ago	t(N(5.3, 0.017))	104.45549
ets	2017 set	t(N(5.3, 0.018))	102.12089
ets	2017 out	t(N(5.3, 0.019))	106.31080
ets	2017 nov	t(N(5.3, 0.021))	102.95220
ets	2017 dez	t(N(5.6, 0.022))	132.85355
ets	2018 jan	t(N(5.2, 0.024))	100.19681
ets	2018 fev	t(N(5.2, 0.025))	94.99414
ets	2018 mar	t(N(5.3, 0.027))	104.35402
ets	2018 abr	t(N(5.3, 0.029))	103.98450
ets	2018 mai	t(N(5.4, 0.03))	110.25586
ets	2018 jun	t(N(5.3, 0.032))	102.64038
ets	2018 jul	t(N(5.3, 0.034))	108.87809
ets	2018 ago	t(N(5.3, 0.036))	108.80236
ets	2018 set	t(N(5.3, 0.038))	106.39732
ets	2018 out	t(N(5.4, 0.04))	110.77283
ets	2018 nov	t(N(5.3, 0.042))	107.30364
ets	2018 dez	t(N(5.7, 0.044))	138.40412
ets	2019 jan	t(N(5.3, 0.046))	104.48207
ets	2019 fev	t(N(5.2, 0.049))	99.09199
ets	2019 mar	t(N(5.3, 0.051))	108.85108
ets	2019 abr	t(N(5.3, 0.054))	108.49006
ets	2019 mai	t(N(5.4, 0.056))	115.03972
ets	2019 jun	t(N(5.3, 0.059))	107.13877
ets	2019 jul	t(N(5.4, 0.061))	113.65663
ets	2019 ago	t(N(5.4, 0.064))	113.60358
ets	2019 set	t(N(5.4, 0.067))	111.12507
ets	2019 out	t(N(5.4, 0.07))	115.70862
ets	2019 nov	t(N(5.4, 0.073))	112.12165
ets	2019 dez	t(N(5.7, 0.076))	144.54630
ets	2020 jan	t(N(5.3, 0.079))	109.23467
ets	2020 fev	t(N(5.3, 0.082))	103.64138
ets	2020 mar	t(N(5.4, 0.086))	113.84480
ets	2020 abr	t(N(5.4, 0.089))	113.49686
ets	2020 mai	t(N(5.5, 0.092))	120.35782
ets	2020 jun	t(N(5.4, 0.096))	112.14487
ets	2020 jul	t(N(5.4, 0.1))	118.97638
ets	2020 ago	t(N(5.4, 0.1))	118.95213
ets	2020 set	t(N(5.4, 0.11))	116.39581
ets	2020 out	t(N(5.5, 0.11))	121.21360
ets	2020 nov	t(N(5.4, 0.12))	117.49951
ets	2020 dez	t(N(5.8, 0.12))	151.39695
ets	2021 jan	t(N(5.4, 0.12))	114.54669

Com base nesse modelo, se entendermos ser o correto, ETS(A,A,A), então posso usar agora para a série completa e estimar forecasts até fim de 2021:

varejo.ts %>%
    model(ets = ETS(box_cox(value, 0.05583817))) %>%
    forecast(h = "2 years") %>%
    autoplot(varejo.ts)

E os nossos valores de forecast são:

varejo.ts %>%
    model(ets = ETS(box_cox(value, 0.05583817))) %>%
    forecast(h = "2 years")

Se desejarem ser mais precisos e forçar o modelo para executar o ETS(A,A,A), sem deixar livre para ele especificar o melhor ETS, faça o seguinte código. Observar que os termos entrarão como fórmula e o sinal de + é para indicar o acréscimo do termo e não que é um método “aditivo”. Nosso forecast então será:

set.seed(1)
varejo.ts %>%
    model(ets = ETS(box_cox(value, 0.05583817) ~ error(method = c("A")) + trend(method = c("A")) +
        season(method = c("A")))) %>%
    forecast(h = "2 years")

Aparece uma pequena diferença comparada ao resultado do pacote forecast:

set.seed(1)
forecast::forecast(ets(varejoms, lambda = 0.05583817, model = "AAA"), h = 24)

         Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
Feb 2020       85.65140  82.04965  89.40207  80.20137  91.44979
Mar 2020       93.03698  88.69205  97.58233  86.47026 100.07257
Apr 2020       88.19589  83.65127  92.97291  81.33592  95.59958
May 2020       93.76184  88.51928  99.29659  85.85732 102.34997
Jun 2020       90.45723  84.98494  96.26103  82.21614  99.47395
Jul 2020       94.65794  88.53002 101.18480  85.43968 104.80966
Aug 2020       96.86060  90.17883 104.00790  86.82031 107.99020
Sep 2020       95.09669  88.12289 102.58928  84.62952 106.77787
Oct 2020       97.75182  90.18002 105.92110  86.39941 110.50245
Nov 2020      101.06224  92.82444 109.98690  88.72447 115.00753
Dec 2020      119.91457 109.73493 130.98125 104.68209 137.22314
Jan 2021       95.20096  86.62441 104.57488  82.38494 109.88336
Feb 2021       86.44426  77.94412  95.81437  73.77024 101.15490
Mar 2021       93.87700  84.31768 104.45320  79.63701 110.49769
Apr 2021       88.97855  79.54273  99.46430  74.93841 105.47701
May 2021       94.57439  84.20816 106.13733  79.16458 112.78637
Jun 2021       91.22694  80.85272 102.84911  75.82230 109.55399
Jul 2021       95.44521  84.24680 108.03881  78.83293 115.32495
Aug 2021       97.64901  85.82856 110.99488  80.13163 118.73920
Sep 2021       95.85597  83.87193 109.44409  78.11528 117.35405
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]

Lembrar que usamos dados apenas até Jan/2020 na série original. Se usar a série atualizada até agora, você terá mais observações que as 241 iniciais.

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, Rob. (2019). fpp3: Data for “Forecasting: Principles and Practice” (3rd Edition). R package. Disponível em: https://github.com/robjhyndman/fpp3-package, https://OTexts.org/fpp3/.

HYNDMAN, R.J.; ATHANASOPOULOS, G. (2020) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. Disponível em: https://otexts.com/fpp3/. Accessed on 02 Apr 2020.