Séries Temporais com R: Análise do Consumo do Varejo em MS com o método THETA e 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 o método THETA e fable e feasts. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/fable_theta_varejoms.

1 Introdução

Neste arquivo utilizo a série do Índice de volume de vendas no varejo Total de Mato Grosso do Sul (média 2014=100), 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 meu Video teórico de Theta no Youtube em: Séries temporais: método Theta com RStudio https://youtu.be/MOImzhmLyQU.

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,
main = "Índice de volume de vendas no varejo total de Mato Grosso do Sul <br> 
        (Mensal)  (2014=100) BCB 1479",
sub="Fonte: Banco Central do Brasil, via pacote BETS.") 

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)  (2014=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)

Fonte: Banco Central do Brasil, via pacote BETS.

É 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.

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. A função fable::THETA realizará a estimação do método theta de Assimakopoulos e Nikolopoulos (2000), que segundo Hyndman e Billah (2003), equivale a um modelo de suavização exponencial simples com drift. A função report fornecerá a saída do modelo.

fit <- train %>%
  model(theta = THETA(value))
fit

report(fit)

Series: value 
Model: THETA 

Alpha: 0.7013 
Drift: 0.194 
sigma^2: 4.5622 

Entre os passos da estimação, a série é testada para sazonalidade pelo teste descrito em Assimakopoulos e Nikolopoulos (2000). Se acusar sazonalidade, a série será ajustada sazonalmente por uma decomposição clássica multiplicativa antes de aplicar o método theta. Os forecasts resultantes serão corrigidos (re-sazonalizados) para apresentação final (O’HARA-WILD et al, 2021).

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
theta	2016 fev	N(86, 4.6)	86.48225
theta	2016 mar	N(96, 6.8)	95.79919
theta	2016 abr	N(95, 9.1)	94.69419
theta	2016 mai	N(101, 11)	100.58242
theta	2016 jun	N(93, 14)	93.31562
theta	2016 jul	N(99, 16)	99.02754
theta	2016 ago	N(99, 18)	98.77018
theta	2016 set	N(96, 20)	96.40543
theta	2016 out	N(101, 23)	100.51970
theta	2016 nov	N(97, 25)	96.56830
theta	2016 dez	N(126, 27)	126.32847
theta	2017 jan	N(94, 29)	94.25584
theta	2017 fev	N(89, 31)	88.54820
theta	2017 mar	N(98, 34)	98.08317
theta	2017 abr	N(97, 36)	96.94734
theta	2017 mai	N(103, 38)	102.97094
theta	2017 jun	N(96, 40)	95.52720
theta	2017 jul	N(101, 43)	101.36987
theta	2017 ago	N(101, 45)	101.10183
theta	2017 set	N(99, 47)	98.67679
theta	2017 out	N(103, 49)	102.88334
theta	2017 nov	N(99, 52)	98.83459
theta	2017 dez	N(129, 54)	129.28740
theta	2018 jan	N(96, 56)	96.45924
theta	2018 fev	N(91, 58)	90.61416
theta	2018 mar	N(100, 61)	100.36715
theta	2018 abr	N(99, 63)	99.20050
theta	2018 mai	N(105, 65)	105.35947
theta	2018 jun	N(98, 67)	97.73879
theta	2018 jul	N(104, 70)	103.71220
theta	2018 ago	N(103, 72)	103.43348
theta	2018 set	N(101, 74)	100.94814
theta	2018 out	N(105, 76)	105.24699
theta	2018 nov	N(101, 79)	101.10088
theta	2018 dez	N(132, 81)	132.24632
theta	2019 jan	N(99, 83)	98.66265
theta	2019 fev	N(93, 85)	92.68011
theta	2019 mar	N(103, 88)	102.65112
theta	2019 abr	N(101, 90)	101.45365
theta	2019 mai	N(108, 92)	107.74799
theta	2019 jun	N(100, 94)	99.95037
theta	2019 jul	N(106, 97)	106.05453
theta	2019 ago	N(106, 99)	105.76512
theta	2019 set	N(103, 101)	103.21950
theta	2019 out	N(108, 103)	107.61064
theta	2019 nov	N(103, 106)	103.36718
theta	2019 dez	N(135, 108)	135.20525
theta	2020 jan	N(101, 110)	100.86605
theta	2020 fev	N(95, 112)	94.74606
theta	2020 mar	N(105, 115)	104.93510
theta	2020 abr	N(104, 117)	103.70681
theta	2020 mai	N(110, 119)	110.13652
theta	2020 jun	N(102, 121)	102.16195
theta	2020 jul	N(108, 123)	108.39686
theta	2020 ago	N(108, 126)	108.09677
theta	2020 set	N(105, 128)	105.49085
theta	2020 out	N(110, 130)	109.97429
theta	2020 nov	N(106, 132)	105.63347
theta	2020 dez	N(138, 135)	138.16418
theta	2021 jan	N(103, 137)	103.06945

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

varejo.ts %>%
  model(theta = THETA(value)) %>% 
  forecast(h = "2 years") %>% 
  autoplot(varejo.ts)

E os nossos valores de forecast são:

varejo.ts %>%
  model(theta = THETA(value)) %>% 
  forecast(h = "2 years")

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.

O’HARA-WILD, Mitchell; HYNDMAN, Rob J.; WANG, Earo. (2021). fable: Forecasting Models for Tidy Time Series. R package version 0.3.1. Disponível em: https://CRAN.R-project.org/package=fable. Accessed on 20 May 2021.

O’HARA-WILD, Mitchell; HYNDMAN, Rob J.; WANG, Earo. (2021). feasts: Feature Extraction and Statistics for Time Series. R package version 0.2.1. Disponível em: https://CRAN.R-project.org/package=feasts. Accessed on 20 May 2021.