INTRODUÇÃO

A necessidade de se realizar previsões a partir de dados coletados ao longo de intervalos de tempo é cada vez maior, ainda mais tratando-se de temas de interesse público, como aqueles relativos à economia, hidrologia/climatologia e, até mesmo, à saúde pública. Para a realização de tais análises, muitos e diversificados são os métodos que podem ser empregados. Neste tutorial, o método utilizado será o ARIMA que, de forma geral, pode ser definido como "modelo autorregressivo integrado de móveis móveis.

Para a aplicação desse modelo, faz-se necessário o entendimento dos termos “p”, “d” e “q”, pois o ajuste a ser realizado será em função de tais variáveis. As definições são as seguintes:

Dependendo do objeto com o qual se está trabalhando, ou seja, a depender da série temporal, o modelo ARIMA também poderá ser denominado como SARIMA. A adição da letra “S” advém da presença de sazonalidade na série.

Abaixo encontram-se as fórmulas:

\[ ARIMA(p,d,q) \]

\[ SARIMA(p,d,q)(P,D,Q)_x \]

A diferenciação entre letras maiúsculas e minúsculas se deve unicamente à sazonalidade do modelo. Enquanto as minúsculas referem-se às componentes não-sazonais, as maiúsculas referem-se às sazonais. O entendimento da terminologia é útil para a execução de análises e, principalmente, de previsões temporais.

IMPLEMENTAÇÃO EM R

Para o ajuste do ARIMA em ambiente R, deve-se recorrer ao pacote forecast.

library(forecast)

Após a ativação do pacote acima, pode-se ter acesso à função Arima(), cuja descrição dos termos que o compõe se encontra abaixo e pode ser acessada da seguinte forma: ?Arima.

Arima(y, # Objeto que será analisado; nesse caso, uma série temporal
  order = c(0, 0, 0),  # Componente não-sazonal
  seasonal = c(0, 0, 0), # Componente sazonal
  xreg = NULL, 
  include.mean = TRUE,
  include.drift = FALSE,
  include.constant,
  lambda = model$lambda,
  biasadj = FALSE,
  method = c("CSS-ML", "ML", "CSS"),  
  model = NULL,
  x = y,
)

Para que os valores dos termos “p”, “d” e “q” sejam encontrados, e o ajuste seja viabilizado, deve-se recorrer às funções de autocorrelação e de autocorrelação parcial. Suponha que a série com a qual se esteja trabalhando é a seguinte:

plot.ts(AirPassengers) # Série temporal dafault

A série acima, denominada AirPassengers, descreve o número de passageiros contabilizados mensalmente ao longo do período de 1949 - 1960. O modelo já se encontra implementado dentre os dataset do R.

Para a verificação das funções de autocorrelação e autocorrelação parcial, pode-se proceder da seguinte forma:

ggtsdisplay(AirPassengers)

A partir dos gráficos:

Para a determinação da ordem de autoregressão (p), do grau de diferenciação (d) e da componente de médias móveis (q), deve-se realizar a observação do comportamento dos “lags” manifestados por cada gráfico de autocorrelação (ACF e ÁCF). “Nos gráficos gerados, as linhas tracejadas azuis são os limites de significância ou intervalos de confiança; sempre que houver uma ultrapassagem, diz-se que há, ali, um lag com significância”

Para que o número de diferenciações seja encontrado com maior facilidade, as funções ndiffs() e nsdiffs() podem ser empregadas. Ambas retornam a quantidade de diferenciações necessárias para a estabilização da estacionariedade e da variânica da série.

A determinação dos valores dos termos pode não ser exata, uma vez que o autor da análise pode incorrer em diversos erros e gerar diversos tipos de modelos. Devido a isso, é necessária a utilizaçaõ de uma métrica com o o intuito de se aferir a qualidade do modelo gerado. Usualmente, se utiliza os valors de AIC e BIC, por exemplo, gerados após o ajuste do modelo.

Pois bem, entendendo essa dificuldade, dentre as ferramentas disponibilizadas pelo pacote forecast(), está a função auto.arima().

AUTO-ARIMA

Como dito anteriormente, essa função realiza a verificação dos possíveis modelos gerados a partir da série temporal em questão, visando ao ajuste ideal. Veja abaixo:

?auto.arima()

# Description
# Returns best ARIMA model according to either AIC, AICc or BIC value. 
# The function conducts a search over possible model within the order constraints provided.

# Usage
auto.arima(
  y, # Série temporal
  d = NA, # Valor de "d" não sazonal estipulado
  D = NA, # Valor de "D" sazonal estipulado
  max.p = 5, 
  max.q = 5,
  max.P = 2,
  max.Q = 2,
  max.order = 5,
  max.d = 2,
  max.D = 1,
  start.p = 2,
  start.q = 2,
  start.P = 1,
  start.Q = 1,
  stationary = FALSE,
  seasonal = TRUE,
  ic = c("aicc", "aic", "bic"), # Critério de seleção do modelo
  stepwise = TRUE, # Se TRUE, buscará a seleção mais rápida
  nmodels = 94,
  trace = FALSE,
  approximation = (length(x) > 150 | frequency(x) > 12),
  method = NULL,
  truncate = NULL,
  xreg = NULL,
  test = c("kpss", "adf", "pp"),
  test.args = list(),
  seasonal.test = c("seas", "ocsb", "hegy", "ch"),
  seasonal.test.args = list(),
  allowdrift = TRUE,
  allowmean = TRUE,
  lambda = NULL,
  biasadj = FALSE,
  parallel = FALSE,
  num.cores = 2,
  x = y,
  ...
)

Para que uma verificação seja mais satisfatória e completa, coloque a função stepwise = FALSE. O processo, a depender do tamanho da série analisada, poderá ser demorado.

Realizando o ajuste da série temporal AirPassengers com parte dos dados, para que se possa comparar os dados gerados.

modelo_treino <- ts(AirPassengers[1:120]) # Amostra treino
modelo_teste <- ts(AirPassengers[121:144]) # Amostra para confirmação
modelo <- auto.arima(modelo_treino, stepwise = FALSE) # Criação do modelo

A partir da geração do modelo, os seus atributos podem ser exibidos.

modelo
## Series: modelo_treino 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2   drift
##       1.6431  -0.9066  -1.8743  0.9688  2.1761
## s.e.  0.0371   0.0361   0.0572  0.0595  0.7107
## 
## sigma^2 estimated as 488.2:  log likelihood=-537.51
## AIC=1087.03   AICc=1087.78   BIC=1103.7

Como pode-se perceber, o modelo ajustado foi um ARIMA, devido à ausência de componentes sazonais. Verificando a aproximação dos dados.

accuracy(modelo$fitted, modelo_treino)
##                 ME     RMSE      MAE        MPE     MAPE       ACF1 Theil's U
## Test set 0.3246097 21.53481 17.21691 -0.5203958 7.005998 0.08174572 0.8001771

Os termos resultantes da acurácia acima são os seguintes:

Agora, se realizará a previsão para os próximos 24 períodos, o que permitirá a visualização da acurácia do modelo para a realização de previsões.

previsao <- forecast(modelo, h = 24)

Os dados gerados são os seguintes:

head(previsao$lower) # Limite inferior
## Time Series:
## Start = 121 
## End = 126 
## Frequency = 1 
##          80%      95%
## 121 314.4380 299.4384
## 122 324.9450 306.0139
## 123 347.5203 327.4080
## 124 372.3897 352.1691
## 125 391.5983 371.3388
## 126 400.8089 380.3839
head(previsao$upper) # Limite superior
## Time Series:
## Start = 121 
## End = 126 
## Frequency = 1 
##          80%      95%
## 121 371.1076 386.1071
## 122 396.4682 415.3992
## 123 423.5063 443.6186
## 124 448.7850 469.0056
## 125 468.1406 488.4001
## 126 477.9761 498.4010

Visualização:

autoplot(previsao, predict.colour = "red") + labs(x = "Tempo", y = "Passageiros", 
                                                  title = "Previsao usando o forecast") +
  theme_test() + scale_x_continuous(breaks = seq(0,150,15)) +
  scale_y_continuous(breaks = seq(100, 600, 50)) + geom_vline(xintercept = 120, lty = "dashed")

Verificação da acurácia da previsão com a amostra teste:

accuracy(modelo_teste, as.numeric(previsao$mean))
##                 ME     RMSE      MAE       MPE     MAPE
## Test set -37.34932 67.58919 46.91646 -8.559127 10.94617

De acordo com a acurácia, podemos perceber que o modelo arima não é infalível, mas dependendo da qualidade do ajuste de dados, pode ser uma ferramenta poderosa para a análise de séries temporais. Cabe ressaltar que há outras ferramentas com as quais se pode realizar previsões, como as redes neurais, por exemplo.