Meu primeiro pacote R: aplicação para série temporal do varejo

Depois de um longo hiato devido à falta de tempo, o blog está de volta à ativa.

Um dos (muitos) motivos de minha ausência tem sido a elaboração do meu TCC, que é sobre previsão de demanda. Eu desenvolvi um sistema que seleciona automaticamente o melhor modelo de previsão dentre os disponíveis no pacote forecast para uma dada série temporal de acordo com a métrica de erro escolhida pelo usuário. O nome do pacote é mafs e já está disponível em meu Github para ser baixado e instalado gratuitamente. Notem, porém, que este é meu primeiro pacote R e eu provavelmente acabei cometendo muitos erros de principiante. Por isso, o pacote ainda é limitado e pode não funcionar em algumas situações que eu não vislumbrei. Uma possível limitação do pacote, por exemplo, é que ele só foi testado para séries mensais e não de outros períodos, como semanais, diárias ou trimestrais.

Demonstração do pacote

Apresentação e análise exploratória dos dados

Para demonstrar na prática como funciona o pacote, irei analisar neste post uma série temporal de periodicidade mensal referente ao volume de vendas do varejo, tema que tenho pesquisado recentemente e obtida no site do IBGE. O dataset será disponibilizado no repositório do blog.

# carregar bibliotecas importantes
library(mafs)
library(magrittr)
library(forecast)
library(ggplot2)
# importar dados
varejo <- read.csv2("https://raw.githubusercontent.com/sillasgonzaga/sillasgonzaga.github.io/master/data/varejo.csv",
                    stringsAsFactors = FALSE)
# exibir dados
head(varejo); tail(varejo)
##   Período Moveis.e.eletrodomesticos
## 1  jan/00                      29.5
## 2  fev/00                      28.4
## 3  mar/00                      30.1
## 4  abr/00                      28.8
## 5  mai/00                      34.8
## 6  jun/00                      30.5
##     Período Moveis.e.eletrodomesticos
## 185  mai/15                     103.8
## 186  jun/15                      91.6
## 187  jul/15                      94.3
## 188  ago/15                      91.1
## 189  set/15                      90.2
## 190                                NA
# retirar última linha, que veio em branco
varejo <- varejo[1:(nrow(varejo)-1), ]

Como pode-se ver, a série temporal vai desde Janeiro de 2000 a Setembro de 2015. Essa informação é importante para criar um objeto da class ts que será usado como input das funções do pacote mafs.

# transformar para série temporal
varejo <- ts(varejo[, 2], start = c(2000, 1), frequency = 12)
# Visualizar série
plot(varejo)

# Visualizar decomposição sazonal da série
varejo %>% decompose %>% plot

O gráfico da série decomposta mostra que há fortes componentes de tendência e sazonalidade na série. O componente aleatório possui média de 0,13, o que, por ser próxima a zero, nos leva a acreditar que a decomposição foi bem sucedida. O elemento sazonal da série também pode ser analisado nos gráficos a seguir.

# função ggmonthplot do pacote forecast
ggmonthplot(varejo)

# estratificação por mês
ggseasonplot(varejo, year.labels = TRUE) + geom_point() + theme_bw()

A partir dos dois gráficos é possível fazer uma observação interessante: A tendência é praticamente alternada. A série sempre cai de Janeiro a Fevereiro, sobe em Março, cai em Abril, sobe em Maio, cai em Junho, sobe ou se mantém estável em Julho, sobe em Agosto, cai ou se mantém estável em Setembro, e sobe de Outubro a Dezembro. A diferença mais evidente ocorrente entre os meses de Novembro e Dezembro.

Poderiam ser feitas mais algumas análises exploratórias, mas eu acabaria fugindo do escopo do post.

Aplicação do modelo.

O pacote mafs é um wrapper de diversos modelos presentes no pacote forecast, que são:

available_models()
##  [1] "auto.arima" "ets"        "nnetar"     "tbats"      "bats"      
##  [6] "stlm_ets"   "stlm_arima" "StructTS"   "meanf"      "naive"     
## [11] "snaive"     "rwf"        "rwf_drift"  "splinef"    "thetaf"    
## [16] "croston"    "tslm"       "hybrid"

Cada um desses modelos pode ser aplicado à série temporal analisada por meio da função mafs::apply_selected_model(). Por exemplo, para o modelo de redes neurais, temos:

apply_selected_model(varejo, "nnetar", horizon = 6) %>% forecast(h  = 6) %>% plot

Imagine-se agora na situação onde vocë é um analista de previsão e precisa realizar, periodicamente, projeções de centenas ou milhares de séries temporais. Seria impraticável testar todos esses 18 modelos disponíveis, não seria? Pensando nisso, a principal função do mafs, chamada select_forecast() automatiza esse processo. Ela depende de quatro parâmetros:
* x, que é a série temporal de input;
* test_size, que é o tamanho da série de teste a ser usado para mensurar a acurácia das previsões obtidas;
* horizon, o tamanho do horizonte de previsão; * error, a métrica de erro para definir o melhor modelo.

O código da função pode ser conferido aqui. Resumidamente, ela separa a série de input em duas: a série de treino, usada para construir o ajuste dos modelos, e a série de teste, usada para mensurar a previsão obtida com os ajustes nas séries de treino em comparação com a série original. A partir das previsões obtidas, a de melhor acurácia (de acordo com a métrica escolhida pelo usuário) é selecionada para prever os valores futuros da série.

Após fazer tudo isso, a função retorna como output três objetos, como pode ser conferido em sua documentação (help("select_forecast")).

output <- select_forecast(varejo, test_size = 6, horizon = 6, error = "MAPE")
# output com resultado de modelos
output$df_models
##         model        ME     RMSE       MAE       MPE     MAPE     MASE
## 1  auto.arima -11.87629 12.25375 11.876285 -12.63291 12.63291 1.683951
## 2        bats -14.05779 14.38923 14.057791 -15.05612 15.05612 1.993269
## 3     croston -22.86110 23.31736 22.861101 -24.61078 24.61078 3.241499
## 4         ets -10.50871 11.15770 10.508709 -11.27795 11.27795 1.490041
## 5      hybrid -14.06093 14.34473 14.060933 -15.01467 15.01467 1.993714
## 6       meanf  25.49699 25.90687 25.496995  26.97627 26.97627 3.615245
## 7       naive  -9.45000 10.50579  9.583333 -10.30420 10.43265 1.358831
## 8      nnetar -19.77968 22.52728 19.779681 -20.68638 20.68638 2.804582
## 9         rwf  -9.45000 10.50579  9.583333 -10.30420 10.43265 1.358831
## 10  rwf_drift -10.87115 11.95716 10.871154 -11.83788 11.83788 1.541432
## 11     snaive -18.16667 18.49784 18.166667 -19.29829 19.29829 2.575871
## 12    splinef -26.14144 26.54988 26.141436 -28.11256 28.11256 3.706622
## 13 stlm_arima -21.31549 21.63862 21.315487 -22.84688 22.84688 3.022345
## 14   stlm_ets -17.70250 18.02916 17.702498 -18.98920 18.98920 2.510056
## 15   StructTS -14.82386 15.16935 14.823861 -15.77388 15.77388 2.101891
## 16      tbats -13.34142 13.85944 13.341419 -14.26644 14.26644 1.891694
## 17     thetaf -14.75199 15.18143 14.751994 -15.76853 15.76853 2.091701
## 18       tslm -25.84932 26.10141 25.849322 -27.70415 27.70415 3.665202
##           ACF1 best_model runtime_model
## 1  -0.24342789      naive         1.471
## 2   0.36081896      naive         3.064
## 3  -0.21035874      naive         0.975
## 4   0.40036466      naive         1.547
## 5   0.11501661      naive         9.041
## 6  -0.21035874      naive         0.003
## 7  -0.21035874      naive         0.004
## 8  -0.09009255      naive         0.984
## 9  -0.21035874      naive         0.002
## 10 -0.06220684      naive         0.004
## 11 -0.25244021      naive         0.003
## 12 -0.18889000      naive         0.437
## 13  0.50840456      naive         0.103
## 14  0.50385963      naive         0.029
## 15 -0.26397115      naive         0.430
## 16  0.30026636      naive         6.824
## 17  0.19311316      naive         0.011
## 18  0.50145846      naive         0.003
# output com valores previstos e reais
output$df_comparison
##         time forecasted observed
## 1 2015-04-03      103.4     92.7
## 2 2015-05-03      103.4    103.8
## 3 2015-06-02      103.4     91.6
## 4 2015-07-03      103.4     94.3
## 5 2015-08-02      103.4     91.1
## 6 2015-09-02      103.4     90.2
# output com valores previstos, incluindo o intervalo de confiança de 80 e de 95%
output$best_forecast
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Oct 2015           90.2 72.52035 107.8796 63.16131 117.2387
## Nov 2015           90.2 65.19720 115.2028 51.96152 128.4385
## Dec 2015           90.2 59.57795 120.8220 43.36762 137.0324
## Jan 2016           90.2 54.84070 125.5593 36.12262 144.2774
## Feb 2016           90.2 50.66710 129.7329 29.73965 150.6603
## Mar 2016           90.2 46.89388 133.5061 23.96901 156.4310

O output de output$df_models mostra que o modelo de menor MAPE foi curiosamente o naive, que corresponde simplesmente a usar o último valor observado como previsão dos próximos valores. Tal previsão pode ser conferida visualmente com outra função do mafs, chamada gg_fit()

gg_fit(varejo, 6, "naive") + theme_bw() + theme(legend.position = "bottom")

Para avaliar a eficiência do meu método, pode-se calcular o MAPE real, isto é, o erro relativo médio entre os valores previstos e os reais, presentes no objeto output$df_comparison

x <- output$df_comparison
# Calcular MAPE real
mape_real <- 100 * abs(x$forecasted - x$observed)/x$observed
# mostrar mape mês a mês
mape_real
## [1] 11.5426106  0.3853565 12.8820961  9.6500530 13.5016465 14.6341463
# mostrar mape médio
mean(mape_real)
## [1] 10.43265

Obtivemos um MAPE médio de 10,43%.

Ideias para o futuro

Devido à automatização possibilitada pelo pacote, é possível pensar em diversas outras análises e testes de hipóteses. Por exemplo: o tamanho da série influencia o desempenho do sistema? Isso poderia ser feito variando o argumento test_size, calculando o MAPE real para cada valor do argumento e depois comparando os resultados. Talvez isso tema de um futuro post.

 
comments powered by Disqus