Capítulo 2 Os fundamentos de séries temporais

Para que uma criança aprenda a identificar um gato, é preciso que os pais apontem para vários gatos e digam “gato”. Essa tarefa pode ser repetida algumas poucas vezes por crianças, mas para um modelo de previsão, essa tarefa não é tão simples, e exige uma grande quantidade de dados para que o modelo seja capaz de generalizar bem para novos casos. Para séries temporais, novos casos são as observações futuras.

É importante também que os dados sejam representativos dos novos casos e que eles sejam de boa qualidade. Informações com erros de medidas, presença de valores discrepantes (chamados de outliers) ou cheio de ruídos, podem gerar previsões de baixa qualidade mesmo que o modelo utilizado seja “adequado”.

Modelos podem realizar previsões imprecisas quando quando ele funciona muito bem para os dados de treinamento, mas não é capaz de prever com precisão novas informações.

2.1 O que pode ser previsto?

2.2 Principais modelos de previsão

2.3 O passo-a-passo de realizar uma previsão

2.4 Um Projeto de Previsão de Ponta a Ponta

O fluxo de trabalho necessário para realizar previsões de uma série temporal seguem uma sequência mais ou menos fixa que começa sempre com a visualização dos dados.

2.4.1 Carregando os pacotes necessários

Todo projeto de previsão é iniciado carregando os principais pacotes de manipulação de dados e de previsão. Com o tidyverse temos acesso a uma série de funções de manipulação de dados e de visualização. Muitas vezes os dados que temos a disposição não estão em um formato correto ou precisam de pequenos ajustes para serem utilizados pelos modelos preditivos.

O pacote timetk contém um conjunto de funções para manipulação, visualização e diagnóstico em séries temporais. Já o pacote modeltime possui uma série de funções que tornam conveniente estimar modelos de séries temporais de uma maneira simples e intuitiva.

# Pacotes de manipulação de dados
library(tidyverse)

# Coleção de pacotes para análise de machine learning e estatística
library(tidymodels)

# Pacotes específicos de séries temporais
library(modeltime)
library(timetk)

2.4.2 Dados

Abaixo temos dados diários de acidentes em rodovias federais brasileiras. Vemos que os dados estão no formato desejado: cada linha representa um valor numérico associado a uma data específica.

acidentes <- read_rds("resources/acidentes_estradas_brasil.rds")
acidentes %>% 
  head() %>% 
  knitr::kable(caption = "Base de acidentes em rodovias federais, 2007 a 2021")
Table 2.1: Base de acidentes em rodovias federais, 2007 a 2021
Date acidentes feriado semana_natal covid
2007-01-01 421 1 0 0
2007-01-02 516 0 0 0
2007-01-03 367 0 0 0
2007-01-04 315 0 0 0
2007-01-05 358 0 0 0
2007-01-06 334 0 0 0

A função skim do pacote skimr é uma forma excelente de resumir um banco de dados.

skimr::skim(acidentes)
Table 2.2: Data summary
Name acidentes
Number of rows 5448
Number of columns 5
_______________________
Column type frequency:
Date 1
numeric 4
________________________
Group variables None

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
Date 0 1 2007-01-01 2021-11-30 2014-06-16 5448

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
acidentes 0 1 351.71 147.00 71 218 349 464 1101 ▇▇▃▁▁
feriado 0 1 0.03 0.18 0 0 0 0 1 ▇▁▁▁▁
semana_natal 0 1 0.02 0.12 0 0 0 0 1 ▇▁▁▁▁
covid 0 1 0.12 0.32 0 0 0 0 1 ▇▁▁▁▁

Vemos que a coluna de datas está no formato datetime, e que a base possui informações de acidentes a partir de primeiro de janeiro de 2007 e fim em 30 de novembro de 2021. A variável dependente (acidentes) não possui valores faltantes (missing values), assim como valores negativos ou números muito fora do esperado, o que nos leva a crer que a base não precisa de nenhuma manipulação adicional. Em termos de estatísticas descritivas, temos uma média de 352 acidentes por dia; valor mínimo de 71 ocorrências, e máximo de 1101.

A figura 2.1 mostra a séries de acidentes em rodovias federais durante 2007 e 2021.

acidentes %>%
  plot_time_series(Date,
                   acidentes,
                   .facet_scales = "free_y",
                   .smooth = F,
                   .title = "")

Figure 2.1: Acidentes diários em rodovias federais brasileiras, 2007-2021

Podemos fazer algumas observações sobre a série:

  • É possível identificar dias com um número muito elevado de acidentes (outliers). Estes valores estão associados a datas como Natal e carnaval.

  • É possível identificar algumas mudanças de padrão que representam quebras estruturais na série. A primeira ocorreu em 2015, e representou uma mudança na forma como os acidentes eram reportados. Após 2015, acidentes sem vítimas não precisavam da presença de agentes rodoviários federais para serem registrados.

  • Em março de 2020 temos outra mudança abrupta no comportamento da série, que foi associada às primeiras medidas de prevenção contra o COVID-19.

Para não termos que lidar com as mudanças estruturais ocorridas na série, vamos utilizar apenas os últimos 12 meses de dados para realizar as previsões.

A figura 2.2 mostra a série limitada ao período entre novembro de 2020 e novembro de 2021. A série possui 366 observações, o que parece suficiente para realizar uma previsão simples.

df <- acidentes %>% 
  filter_by_time(
    .start_date = last(Date) %-time% "12 month",
    .end_date = "end"
  ) 
  

df %>%
  plot_time_series(Date,
                   acidentes,
                   .facet_scales = "free_y",
                   .smooth = F)

Figure 2.2: Total de acidentes em rodovias federais no Brasil

Observando apenas os últimos 12 meses temos uma visão melhor do comportamento dos dados. É possível observar que existe um padrão sazonal na série, com certos dias da semana apresentando um número elevado de acidentes.

Para melhor entender estes padrões sazonais, a figura 2.3 mostra o resultado da função plot_seasonal_diagnostics().

df %>%
  plot_seasonal_diagnostics(Date, acidentes)
## Warning: The following aesthetics were dropped during statistical transformation:
## y_plotlyDomain
## ℹ This can happen when ggplot fails to infer the correct grouping structure
##   in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation:
## y_plotlyDomain
## ℹ This can happen when ggplot fails to infer the correct grouping structure
##   in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation:
## y_plotlyDomain
## ℹ This can happen when ggplot fails to infer the correct grouping structure
##   in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation:
## y_plotlyDomain
## ℹ This can happen when ggplot fails to infer the correct grouping structure
##   in the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

Figure 2.3: Padrão sazonal da série de acidentes de trânsito

É possível observar que os sábados e domingos possuem um número de acidentes maior do que a média. Assim, como as terças-feiras estão associadas com um menor número de acidentes.

A figura 2.4 mostra a correlação entre as observações atuais e suas defasagens. Assim, o dia anterior representa o lag 1, assim como o mesmo dia da semana passada representa o lag 7. Este tipo de gráfico se chama correlogramo e será explicado melhor na seção xxx.

df %>%
  plot_acf_diagnostics(Date, acidentes, .lag = 100, .show_white_noise_bars = T)
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the plotly package.
##   Please report the issue at <]8;;https://github.com/plotly/plotly.R/issueshttps://github.com/plotly/plotly.R/issues]8;;>.

Figure 2.4: Correlograma da série de acidentes de trânsito

A figura confirma nossas suspeitas de que a série de acidentes tem um padrão bem distinto de sazonalidade. Especificamente, uma sazonalidade de 7 dias, que é associada com séries que possuem efeito de dia da semana.

2.4.3 Criando uma base de treinamento

Antes de iniciar a modelagem da série temporal, é boa prática dividir a série em dois blocos: o primeiro bloco (chamado de base de treinamento) é onde o modelo aprenderá o comportamento dos acidentes. O segundo bloco (chamado de base de teste) é onde veremos se o modelo fez um bom trabalho em reproduzir o comportamento. Uma explicação mais profunda sobre como realizar essa divisão para séries temporais pode ser encontrada na seção xxxx.

A figura 2.5 mostra a base de treinamento e teste. Utilizaremos os primeiros 8 meses para treinar o modelo e testaremos sua capacidade preditiva ao comparar a previsão com o que foi observado nos últimos 4 meses.

splits_br <- df %>%
  time_series_split(date_var = Date,
                    assess = "4 months",
                    cumulative = TRUE
  )

# visualizar base de treino e teste
splits_br %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(Date, acidentes)

Figure 2.5: Série de treinamento e teste

Um grande perigo de qualquer modelo de série temporal é a de que ele esteja sobre-ajustando os dados, ou seja, simplesmente imitando o comportamento visto no período de treinamento e sendo incapaz de generalizar para períodos futuros. Podemos minimizar este risco ao criar várias bases de treino e teste. Este processo de reamostragem será melhor explicado na seção xxx.

A figura 2.6 mostra que criamos 4 reamostragens diferentes.

splits_reamostra_br <- df %>%
  time_series_cv(Date,
                 assess = "4 months",
                 skip = "2 months",
                 cumulative = TRUE,
                 slice_limit =4)
splits_reamostra_br %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(Date, acidentes)

Figure 2.6: Cross-validation para série de acidentes

2.4.4 Fórmula

Vimos que nossa série possui um padrão sazonal bastante distinto, além de ser influenciada por dias importantes do calendário como Natal e carnaval. Portanto, um modelo que pretenda capturar bem o comportamento dos acidentes de trânsito deve incorporar essas informações.

Assim, vamos criar variáveis indicativas para o dia da semana, o semestre, o trimestre, a semana do ano, o dia do ano e o dia do trimestre, além das variáveis indicativas de feriado. Na seção xxx será explicado em detalhes o processo de criação destas variáveis com o pacote recipe.

# formula para modelos de séries temporais
# precisa uma coluna de data como covariada

receita_ts <- recipe(acidentes ~ Date + ., data = training(splits_br)) %>%
  step_timeseries_signature(Date) %>%
  step_rm(matches("(xts$)|(iso$)|(^.pm)")) %>%
  step_holiday(Date) %>%
  step_zv(all_predictors()) %>%
  step_mutate(Date_month = factor(Date_month, ordered = TRUE)) %>%
  step_dummy(all_nominal_predictors())

receita_ts %>% 
  prep() %>% 
  bake(new_data = NULL) %>% 
  head() %>% 
  knitr::kable(caption = "Base de dados com novas variáveis")
## Warning: `terms_select()` was deprecated in recipes 0.1.17.
## ℹ Please use `recipes_eval_select()` instead.
## ℹ The deprecated feature was likely used in the timetk package.
##   Please report the issue at
##   <]8;;https://github.com/business-science/timetk/issueshttps://github.com/business-science/timetk/issues]8;;>.
Table 2.3: Base de dados com novas variáveis
Date feriado semana_natal acidentes Date_index.num Date_year Date_half Date_quarter Date_day Date_wday Date_mday Date_qday Date_yday Date_mweek Date_week Date_week2 Date_week3 Date_week4 Date_mday7 Date_LaborDay Date_NewYearsDay Date_ChristmasDay Date_month_1 Date_month_2 Date_month_3 Date_month_4 Date_month_5 Date_month_6 Date_month_7 Date_month_8 Date_month_9 Date_month.lbl_01 Date_month.lbl_02 Date_month.lbl_03 Date_month.lbl_04 Date_month.lbl_05 Date_month.lbl_06 Date_month.lbl_07 Date_month.lbl_08 Date_month.lbl_09 Date_month.lbl_10 Date_month.lbl_11 Date_wday.lbl_1 Date_wday.lbl_2 Date_wday.lbl_3 Date_wday.lbl_4 Date_wday.lbl_5 Date_wday.lbl_6
2020-11-30 0 0 167 1606694400 2020 2 4 30 2 30 61 335 5 48 0 0 0 5 0 0 0 0.3853373 0.1740777 -0.1511417 -0.4113767 -0.5012804 -0.4281744 -0.2751787 -0.1308926 -0.0408164 0.3763089 0.2281037 -0.0418121 -0.3017184 -0.4518689 -0.4627381 -0.3701419 -0.2388798 -0.1236175 -0.0491049 -0.0130968 -0.3779645 0.0000000 0.4082483 -0.5640761 0.4364358 -0.1973855
2020-12-01 0 0 161 1606780800 2020 2 4 1 3 1 62 336 5 48 0 0 0 1 0 0 0 0.4954337 0.5222330 0.4534252 0.3365809 0.2148345 0.1167748 0.0526938 0.0186989 0.0045352 0.4599331 0.5018282 0.4599331 0.3687669 0.2616083 0.1641974 0.0904791 0.0430767 0.0172126 0.0054561 0.0011906 -0.1889822 -0.3273268 0.4082483 0.0805823 -0.5455447 0.4934638
2020-12-02 0 0 150 1606867200 2020 2 4 2 4 2 63 337 1 49 1 1 1 1 0 0 0 0.4954337 0.5222330 0.4534252 0.3365809 0.2148345 0.1167748 0.0526938 0.0186989 0.0045352 0.4599331 0.5018282 0.4599331 0.3687669 0.2616083 0.1641974 0.0904791 0.0430767 0.0172126 0.0054561 0.0011906 0.0000000 -0.4364358 0.0000000 0.4834938 0.0000000 -0.6579517
2020-12-03 0 0 164 1606953600 2020 2 4 3 5 3 64 338 1 49 1 1 1 1 0 0 0 0.4954337 0.5222330 0.4534252 0.3365809 0.2148345 0.1167748 0.0526938 0.0186989 0.0045352 0.4599331 0.5018282 0.4599331 0.3687669 0.2616083 0.1641974 0.0904791 0.0430767 0.0172126 0.0054561 0.0011906 0.1889822 -0.3273268 -0.4082483 0.0805823 0.5455447 0.4934638
2020-12-04 0 0 203 1607040000 2020 2 4 4 6 4 65 339 1 49 1 1 1 1 0 0 0 0.4954337 0.5222330 0.4534252 0.3365809 0.2148345 0.1167748 0.0526938 0.0186989 0.0045352 0.4599331 0.5018282 0.4599331 0.3687669 0.2616083 0.1641974 0.0904791 0.0430767 0.0172126 0.0054561 0.0011906 0.3779645 0.0000000 -0.4082483 -0.5640761 -0.4364358 -0.1973855
2020-12-05 0 0 279 1607126400 2020 2 4 5 7 5 66 340 1 49 1 1 1 1 0 0 0 0.4954337 0.5222330 0.4534252 0.3365809 0.2148345 0.1167748 0.0526938 0.0186989 0.0045352 0.4599331 0.5018282 0.4599331 0.3687669 0.2616083 0.1641974 0.0904791 0.0430767 0.0172126 0.0054561 0.0011906 0.5669467 0.5455447 0.4082483 0.2417469 0.1091089 0.0328976

2.4.5 Especificando o modelo

Com o pacote modeltime podemos especificar o modelo a ser utilizado. O primeiro modelo de previsão será o SNAIVE, ou modelo Naïve Sazonal. Este é um modelo extremamente simples que repete o último valor observado em um dia da semana para todos os valores futuros. Assim, se a última segunda-feira observada produziu 200 acidentes, ele irá repetir 200 em todas as segundas-feiras futuras. Este é um modelo simples, mas serve como um modelo de comparação (baseline) para modelos mais complexos.

O nosso segundo modelo é um modelo de Regressão com Erros Arima, que será explicado em detalhes na seção xxx. Para melhorar a previsão do modelo, realizaremos um procedimento de boosting dos resíduos do modelo (explicado na seção xxx).

O modelo de Regressão com Erros Arima Boosted possui uma série de hiperparâmetros que precisam ser definidos antes do modelo ser ajustado. A escolha destes hiperparâmetros pode produzir modelos com melhor ou pior performance.

# modelo regressao com erros arima boosted
spec_snaive <- naive_reg() %>% 
  set_engine("snaive")


spec_regarima <- arima_boost(  mtry           = tune(),
                               trees          = 1000,
                               min_n          = tune(),
                               tree_depth     = tune(),
                               learn_rate     = tune(),
                               loss_reduction = tune()
                               #sample_size    = tune()
) %>%
  set_engine("auto_arima_xgboost")

wflw_regarima <- workflow() %>%
  add_model(spec_regarima) %>%
  add_recipe(receita_ts)

Para escolhermos o melhor conjuntos de hiperparâmetros precisamos adotar um processo de tuning, que será melhor explicado na seção xxx. Abaixo utilizamos as quatro reamostragens construídas na seção xxx para rodar um número bastante elevado de modelos: para cada um das quatro reamostragens, vamos estimar um modelo para cada combinação de hiperparâmetro.

# Encontrar melhores hiperparâmetros ----
# tune regarima
set.seed(3333)
tune_results_regarima <- tune_grid(
  object     = wflw_regarima,
  resamples  = splits_reamostra_br,
  param_info = parameters(wflw_regarima),
  grid       = 5,
  control    = control_grid(verbose = FALSE,
                            allow_par = TRUE,
                            parallel_over = "everything")
)

Abaixo vemos que para a primeira reamostragem (Slice1) ajustamos o modelo cinco vezes para diferentes valores dos hiperparâmetros mtry, min_n, tree_depth e learn_rate.

tune_results_regarima %>% 
  unnest(cols = ".metrics") %>% 
  filter(.metric == "rmse") %>% 
  head()
## # A tibble: 6 × 12
##   splits            id      mtry min_n tree_de…¹ learn…² loss_…³ .metric .esti…⁴
##   <list>            <chr>  <int> <int>     <int>   <dbl>   <dbl> <chr>   <chr>  
## 1 <split [246/120]> Slice1    35    30        15 2.75e-9 5.33e-9 rmse    standa…
## 2 <split [246/120]> Slice1    14    11         7 7.31e-9 8.76e-8 rmse    standa…
## 3 <split [246/120]> Slice1     7    37         6 5.77e-2 1.49e+1 rmse    standa…
## 4 <split [246/120]> Slice1    23     7        10 4.44e-5 9.67e-3 rmse    standa…
## 5 <split [246/120]> Slice1    44    24         2 7.96e-7 2.98e-4 rmse    standa…
## 6 <split [185/120]> Slice2    35    30        15 2.75e-9 5.33e-9 rmse    standa…
## # … with 3 more variables: .estimate <dbl>, .config <chr>, .notes <list>, and
## #   abbreviated variable names ¹​tree_depth, ²​learn_rate, ³​loss_reduction,
## #   ⁴​.estimator

Podemos tomar uma média do erro de previsão das 4 reamostragens. Ao encontrarmos os hiperparâmetros que produzem o modelo com menor erro (diferença entre o observado e o previsto) podemos ajustar o modelo sobre toda a base de treinamento. Podemos chamar este procedimento de finalização do modelo.

regarima_tuned_best <- tune_results_regarima %>%
  select_best("rmse")

regarima_tuned_best %>% 
  knitr::kable(caption = "Conjunto de hiperparâmetros que produziu modelo com melhor performance")
Table 2.4: Conjunto de hiperparâmetros que produziu modelo com melhor performance
mtry min_n tree_depth learn_rate loss_reduction .config
23 7 10 4.44e-05 0.0096693 Preprocessor1_Model4

Abaixo temos a saída do melhor modelo de Regressão com Erros Arima. Uma explicação mais profunda sobre a interpretação da saída do modelo pode ser encontrada na seção xxx.

wflw_fit_regarima <- wflw_regarima %>%
  finalize_workflow(parameters = regarima_tuned_best) %>%
  fit(training(splits_br))
## frequency = 7 observations per 1 week
wflw_fit_regarima
## ══ Workflow [trained] ════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: arima_boost()
## 
## ── Preprocessor ──────────────────────────────────────────────────────────────
## 6 Recipe Steps
## 
## • step_timeseries_signature()
## • step_rm()
## • step_holiday()
## • step_zv()
## • step_mutate()
## • step_dummy()
## 
## ── Model ─────────────────────────────────────────────────────────────────────
## ARIMA(1,0,0)(0,1,1)[7] w/ XGBoost Errors
## ---
## Model 1: Auto ARIMA
## Series: outcome 
## ARIMA(1,0,0)(0,1,1)[7] 
## 
## Coefficients:
##          ar1     sma1
##       0.4699  -0.8331
## s.e.  0.0634   0.0623
## 
## sigma^2 = 441.4:  log likelihood = -1070.12
## AIC=2146.25   AICc=2146.35   BIC=2156.68
## 
## ---
## Model 2: XGBoost Errors
## 
## xgboost::xgb.train(params = list(eta = 4.43507367237661e-05, 
##     max_depth = 10L, gamma = 0.00966926536984491, colsample_bytree = 1, 
##     colsample_bynode = 0.5, min_child_weight = 7L, subsample = 1), 
##     data = x$data, nrounds = 1000, watchlist = x$watchlist, verbose = 0, 
##     objective = "reg:squarederror", nthread = 1)

Abaixo estimamos o modelo SNAIVE. Este modelo não possui hiperparâmetros e pode ser ajustado diretamente sobre a base de treinamento.

wflw_fit_snaive <- workflow() %>% 
  add_recipe(receita_ts) %>% 
  add_model(spec_snaive) %>% 
  fit(training(splits_br))
## frequency = 7 observations per 1 week

2.4.6 Previsão

Por fim, podemos realizar as previsões dos modelos estimados. É sempre interessante realizar uma primeira previsão para o período de teste para verificar como o modelo se comporta quando comparamos sua previsão com dados que foram observados mas não utilizados durante o treinamento. A seção xxx explica em detalhes como podemos utilizar modelos para realizar previsões.

A figura 2.7 mostra uma comparação entre os valores previstos e observados.

tbl_calibracao <- modeltime_table(wflw_fit_regarima, wflw_fit_snaive) %>%
  modeltime_calibrate(testing(splits_br)) 

tbl_calibracao %>% 
  modeltime_forecast(
    new_data = testing(splits_br),
    actual_data = df,
    keep_data = TRUE
  ) %>% 
  filter_by_time(Date, 
                 .start_date = last(Date) %-time% "3 month",
                 .end_date = "end") %>% 
  plot_modeltime_forecast(.legend_max_width = 10, .title = "")

Figure 2.7: Previsão de acidentes em rodovias federais brasileiras

2.4.7 Performance do modelo

Ambos os modelos parecem ter capturado bem o padrão sazonal da série, mas inspeções visuais das previsões podem ser enganosas. Assim, podemos utilizar medidas de desempenho do modelo que são construídas a partir da comparação entre os valores previstos e observados. A seção xxx traz uma explicação detalhada das medidas de desempenho mais utilizadas, assim como as vantagens e desvantagens de cada uma.

A tabela abaixo mostra diferentes medidas de performance obtida pelos modelos. Valores menores das medidas representam modelos com menor erro (com exceção do rsq). Assim, para todas as medidas utilizadas, o modelo mais complexo (Regressão com Erros Arima), a previsão foi melhor do que a do modelo baseline SNAIVE.

tbl_calibracao %>%
  modeltime_accuracy() %>% 
  table_modeltime_accuracy()
## Warning: `bindFillRole()` only works on htmltools::tag() objects (e.g., div(),
## p(), etc.), not objects of type 'shiny.tag.list'.

Assim, podemos utilizar o modelo de Regressão com Erros ARIMA para realizar previsões para o futuro.