Capítulo 7 Modelos ARIMA
Modelos ARIMA são outra classe de modelos populares em forecasting. A diferença entre modelos de suavização exponencial e os ARIMA é que os primeiros são baseados em descrever a tendência e a sazonalidade na série, enquanto os segundos se baseam em autocorrelações presentes nos dados.
Para este capítulo, continuaremos usando a série temporal do consumo de energia pela indústria no Brasil.
library(forecast)
library(tidyverse)
energia <- readRDS("data/ts_energia.Rda")Antes de explicar o que é ARIMA, precisamos introduzir conceitos básicos como estacionariedade.
7.1 Estacionariedade, diferenciação e autocorrelação
Uma série temporal é dita estacionária se suas propriedades (média, variância, etc.) não dependem do tempo da observação. Portanto, séries que apresentam tendência ou sazonalidade não são estacionárias. Por outro lado, uma série composta por valores gerados aleatoriamente (ex. pela função rnorm) são estacionárias, visto que a “aparência” da série é basicamente a mesma para qualquer período \(t\).
data("lynx")
head(lynx)## Time Series:
## Start = 1821
## End = 1826
## Frequency = 1
## [1] 269 321 585 871 1475 2821autoplot(lynx)
Por curiosidade, linces são isto:
Aparentemente, a série apresenta ciclos ou sazonalidade, certo? Vamos ver seus componentes por dentro:
lynx %>% decompose() %>% autoplot()## Error in decompose(.): time series has no or less than 2 periodslynx %>% stl(s.window = "periodic") %>% autoplot()## Error in stl(., s.window = "periodic"): series is not periodic or has less than two periodsPor que a função retorna um erro? Acontece que esta série não possui sazonalidade. Os algoritmos de decomposição não foram capazes de detectar períodos em que padrões se repetem, visto que os ciclos, além de possuírem níveis não constantes, acontecem em diferentes períodos, tornando-os imprevisíveis.
Voltando a nossa série de exemplo. Ela é estacionária?
energia %>% autoplot()
Claramente, ela não é. Contudo, existe um método muito simples para a tornar estacionária, que é a diferenciação. Ele se baseia nada mais em calcular a diferença absoluta (ou percentual) entre uma observação e a outra. O intervalo entre as observações pode ser escolhido pelo usuário:
head(energia, 10)## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1979 4215 4221 4329 4405 4563 4672 4716 4779 4890 4776head(diff(energia), 10)## Feb Mar Apr May Jun Jul Aug Sep Oct Nov
## 1979 6 108 76 158 109 44 63 111 -114 164head(diff(energia, 2), 10)## Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1979 114 184 234 267 153 107 174 -3 50 120# plotando as series
par(mfrow=c(3,1))
plot(energia, main = "Série original", xlab = "", ylab = "")
plot(diff(energia), main = "Série diferenciada", xlab = "", ylab = "")
# plotando uma serie aleatoria:
plot(rnorm(length(energia)), type = "l", main = "Dados aleatórios", xlab = "", ylab = "")
Vemos pelo gráfico que a série diferenciada se aproxima muito de dados gerados aleatoriamente, comprovando que o método da diferenciação tornou a série estacionária.
Uma outra maneira de analisar a estacionariedade de uma série é pelo conceito de autocorrelação, que é uma medida do relacionamento linear entre os valores dentro de uma mesma série (como se fosse uma correlação interna). A autocorrelação depende do lag que se deseja analisar. Por exemplo, \(r1\) se refere à relação entre \(y_t\) e \(y_{t-1}\).
A melhor maneira de se analisar a autocorrelação de uma série é por meio de gráficos:
lag.plot(energia, lags = 9)
No gráfico de dispersão acima, vimos que em todos os 9 lags analisados, existe uma correlação forte entre as observações entre si. Numericamente, temos:
Acf(energia, plot = FALSE)$acf %>% as.vector()## [1] 1.0000000 0.9893734 0.9790275 0.9695393 0.9612011 0.9538097 0.9470524
## [8] 0.9408976 0.9358154 0.9313553 0.9275797 0.9250635 0.9215740 0.9119906
## [15] 0.9027513 0.8942274 0.8870129 0.8810469 0.8752439 0.8699125 0.8661511
## [22] 0.8616840 0.8584167 0.8556130 0.8510220 0.8416753 0.8318787Outro gráfico muito útil para analisar a autocorrelação é o próprio gráfico de autocorrelação, também conhecido como correlograma:
Acf(energia)
A linha horizonal tracejada azul define o nível de significância mínimo. Isto é, se uma autocorrelação estiver acima dessa linha, ela é estatisticamente significativa.
Contudo, a situação é diferente para a série diferenciada, como mostram os gráficos abaixo:
lag.plot(diff(energia), lags = 9)
Acf(diff(energia))
O fato de ainda haver três valores podem indicar duas coisas:
- Que a série diferenciada ainda não é estacionária;
- Que a série diferenciada é sim estacionária, pois esses lags com forte autocorrelação podem ser devidos a uma aleatoriedade. De fato, espera-se que, para uma série aleatória, 5% dos valores na ACF encontrem-se fora do intervalo \(\pm 2/ \sqrt{T}\).
7.2 Modelos autoregressivos
A diferença entre modelos regressivos e autoregressivos é que os primeiros prevêem o valor de uma variável de interesse usando uma combinação linear (equação) das variáveis explanatórias. Já os segundos usam uma combinação linear de valores passados da própria variável. Matematicamente, um modelo autoregressivo é descrito como:
\(y_t = c + \phi_1y_{t-1} + \phi_2y_{t-2} + ... + \phi_py_{t-p} + e_t\)
Onde \(c\) é uma constante e \(e_t\) é um erro aleatório (ruído branco). Esse tipo de modelo é chamado de modelo AR(p) e são normalmente restritos a séries estacionárias.
No R, modelos autoregressivos podem ser ajustados por meio da função arima:
# explicaremos o que é o argumento order da função abaixo daqui a pouco
ar.p1 <- energia %>% diff() %>% arima(order = c(1, 0, 0))
ar.p2 <- energia %>% diff() %>% arima(order = c(2, 0, 0))
coefficients(ar.p1) %>% round(4)## ar1 intercept
## -0.0142 21.5167coefficients(ar.p2) %>% round(4)## ar1 ar2 intercept
## -0.0151 -0.0621 21.3925Portanto, as equações resultantes do ajuste de modelos autoregressivos de ordem 1 e 2 são, respectivamente:
\(y_t = 21.1135 - 0.0135y_{t-1} + e_t\)
\(y_t = 20.9643 - 0.01435y_{t-1} - 0.0143y_{t-2} + e_t\)
7.3 Modelos de média móvel
Modelos de média móvel utilizam valores passados de erro de previsão de maneira semelhante a um modelo de regressão:
\(y_t = c + e_t + \theta_1e_{t-1} + \theta_2e_{t-2} + ... + \theta_qe_{t-q}\)
Um modelo acima é chamado de modelo MA(q) e pode ser interpretado como um modelo onde \(y_t\) é uma média ponderada dos erros de previsão passados.
7.4 Modelos ARIMA
A combinação entre os métodos de diferenciação e os modelos de autoregressão e média móvel resultam em um modelo ARIMA (AutoRegressive Integrated Moving Average model) não-sazonal, que pode ser descrito matematicamente como:
\(\acute{y_t} = c + \phi_1\acute{y}_{t-1} + ... + \phi_p\acute{y}_{t-p} + \theta_1 + \theta_1e_{t-1} + ... + \theta_qe_{t-q} + e_t\)
Onde \(\acute{y}_t\) é a série diferenciada. A equação acima é o que descreve o modelo ARIMA(p, d, q), onde:
- \(p\) é a ordem do modelo autoregressivo;
- \(d\) é o grau de diferenciação;
- \(q\) é a ordem do modelo de média móvel.
Como seria complicado selecionar valores apropriada para cada um desses três parâmetros, a função forecast::auto.arima() faz isso automaticamente:
mod.arima <- auto.arima(energia, seasonal = FALSE)
summary(mod.arima)## Series: energia
## ARIMA(2,1,1) with drift
##
## Coefficients:
## ar1 ar2 ma1 drift
## 0.8145 -0.0965 -0.8672 20.8063
## s.e. 0.0645 0.0488 0.0468 7.3280
##
## sigma^2 estimated as 111001: log likelihood=-3344.47
## AIC=6698.93 AICc=6699.07 BIC=6719.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.6031951 331.3683 248.7715 -0.01133265 2.525879 0.4451406
## ACF1
## Training set -0.005477991O modelo resultante é um ARIMA(2, 1, 1), o que corresponde a uma combinação a uma série com um grau de diferenciação (p = 1) aplicada em um modelo AR(2) e em um MA(1).
É possível saber o comportamento da previsão de um modelo ARIMA apenas baseado nos valores dos coeficientes c e d:
- Se \(c = 0\) e \(d = 0\), previsões em longo prazo serão iguais a zero;
- Se \(c = 0\) e \(d = 1\), serão iguais a uma constante maior que zero;
- Se \(c \neq 0\) e \(d = 2\), seguirão uma linha reta;
- Se \(c \neq 0\) e \(d = 0\), convergirão para a média da série;
- Se \(c \neq 0\) e \(d = 1\), seguirão uma linha reta;
- Se \(c \neq 0\) e \(d = 2\), seguirão uma tendência quadrática;
No caso do nosso exemplo, vamos ver se isso se aplica:
mod.arima %>% forecast(h = 36) %>% autoplot()
Conforme esperado, as previsões seguem uma linha reta, o que, em palavras, indicam que a série crescerá indefinitivamente.
Como ficaram as previsões de outras combinações ARIMA? Vamos fazer o teste:
# criar funcao para extrair previsao de um modelo arima
prever.arima <- function(p = 0, d, q = 0, con) {
# argumento c = binario TRUE ou FALSE
# argumento d = 0, 1 ou 2
x <- Arima(energia, order = c(p, d, q), include.constant = con, seasonal = c(0, 0, 0))
x <- forecast(x)
as.numeric(x$mean)
}
# plotar previsoes
prever.arima(d = 1, con = FALSE) %>% plot(type = "l")
prever.arima(d = 0, con = FALSE) %>% lines(col = "red")
prever.arima(d = 2, con = FALSE) %>% lines(col = "orange")
prever.arima(d = 0, con = TRUE) %>% lines(col = "green")
prever.arima(d = 1, con = TRUE) %>% lines(col = "pink")
prever.arima(d = 2, con = TRUE) %>% lines(col = "blue")
7.5 Modelos SARIMA
Modelos ARIMA são capazes também de modelar séries que apresentam um componente sazonal, sendo descrito como:
ARIMA \((p, d, q)(P, D, Q)_m\)
Onde o primeiro parênteses se refere à parte não-sazonal do modelo e o segundo à parte sazonal. \(m\) corresponde ao número de períodos sazonais.
Ajustar um modelo SARIMA é semelhante ao processo de ajustar um ARIMA:
mod.sarima <- auto.arima(energia, seasonal = TRUE)
mod.sarima## Series: energia
## ARIMA(1,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 sar1 sar2
## -0.1114 0.4288 0.2269
## s.e. 0.0471 0.0446 0.0461
##
## sigma^2 estimated as 75746: log likelihood=-3259.25
## AIC=6526.49 AICc=6526.58 BIC=6543.04mod.sarima %>% forecast(h = 36) %>% autoplot()
Para a série de análise, qual modelo possui o melhor ajuste: o sazonal ou o não-sazonal?
# calculando a qualidade de ajuste:
list(mod.arima, mod.sarima) %>% map(accuracy)## [[1]]
## ME RMSE MAE MPE MAPE MASE
## Training set 0.6031951 331.3683 248.7715 -0.01133265 2.525879 0.4451406
## ACF1
## Training set -0.005477991
##
## [[2]]
## ME RMSE MAE MPE MAPE MASE
## Training set 7.740118 274.0315 201.3429 0.06620215 2.088898 0.3602739
## ACF1
## Training set -0.002166188# calculando a qualidade preditiva
treino <- window(energia, end = c(2015, 12))
teste <- window(energia, start = c(2016, 1))
mod.arima2 <- auto.arima(treino, seasonal = FALSE) %>% forecast(teste)## Warning in 1:h: numerical expression has 20 elements: only the first usedmod.sarima2 <- auto.arima(treino, seasonal = TRUE) %>% forecast(teste)
list(mod.arima2, mod.sarima2) %>% map(forecast) %>% map(accuracy, teste)## [[1]]
## ME RMSE MAE MPE MAPE MASE
## Training set 0.6579 330.8885 247.8696 -0.01504425 2.546858 0.4341097
## Test set -432.8558 566.6038 440.5747 -3.22377046 3.278950 0.7716063
## ACF1 Theil's U
## Training set -0.005512616 NA
## Test set 0.626147776 1.602787
##
## [[2]]
## ME RMSE MAE MPE MAPE MASE
## Training set 4.994897 276.0395 202.3100 0.04766251 2.123846 0.3543183
## Test set 926.070363 1027.7004 946.0789 6.70087421 6.860102 1.6569277
## ACF1 Theil's U
## Training set -0.001806213 NA
## Test set 0.624322671 3.055618Apesar de o modelo ARIMA possuir o pior qualidade de ajuste, seu desempenho para prever valores futuros é superior ao SARIMA.
