Introdução

O objetivo dessa publicação é explicar a utilização de um modelo SARIMA, de forma conceitual e prática. Não tem o objetivo de gerar uma previsão com alta acurácia, mas sim mostrar o passo a passo da tomada de decisão e construição de um modelo preditivo univariado.

Conceitos

Introdução a alguns conceitos que serão úteis no decorrer da análise.

ARIMA

Primeiramente, se faz necessário, o entendimento conceitual referente ao modelo ARIMA, que é um modelo muito popular quando o assunto é previsão de séries temporais univariadas. Podemos dividi-lo em três partes, sendo elas:

Autoregressivo (AR), que significa que a variável será analisada baseada em seus próprios valores anteriores. Integrado (I), mostra a diferenciação dos dados com os valores atuais e anteriores. Média móvel (MA), que significa modelar os termos de erros baseado em uma combinação linear, ou seja, multiplicar os termos por uma constante.

Parâmetros ARIMA:

Baseado na fórmula, podemos ver os parâmetros do modelo, onde:

p é o número de termos autoregressivos; d é o número de diferenças; q é o número de termos da média-móvel;

Assim, podemos dizer que em seus parâmetros (p,d,q), cada letra representa uma parte do modelo e são representadas por um número inteiro não negativo. Por exemplo, um ARIMA (1,0,1), se trata de um modelo ARMA, pois o “d” mostra que o modelo não é integrado ou um ARIMA(1,0,0) se trata de um modelo AR(1), pois “d” e “q”, estão zerados, dessa forma o modelo é apenas autoregressivo.

Séries temporais - Sazonalidade

Modelos em séries temporais, são estudos que procuram analisar variáveis no passado com o objetivo de entendê-las no presente ou prevê-las no futuro. Assim, dentro de séries temporais é importante entender o conceito de estacionariedade:

Podemos dizer que uma série é estacionária, quando suas propriedades não tem relação com o tempo. Ou seja, matematicamente, ela flutua em torno da mesma média ao longo do período da série. Na prática, os “choques” tendem a não afetar suas propriedades (média, variância, etc) e se estabilizar com o passar do tempo.

E, ao contrário do que temos para séries estacionárias, podemos ter algo não estacionário. Um exemplo desse caso, são séries com sazonalidade, caso peguemos uma série temporal que mostre as vendas de blusas, podemos supor que a quantidade de produtos vendidos aumente no inverno. Esse tipo de produto é conhecido como produto sazonal, que tem maior procura em determinado período do ano (O mercado de vestuário costuma ser um bom exemplo de sazonalidade).

Lidando com a sazonalidade - Modelo SARIMA

Cada negócio terá sua própria sazonalidade, o de vestuário, pode ser nas estações do ano, mas um vendedor de sorvete na praia, pode vender mais de acordo com a temperatura do dia, nesse caso uma sazonalidade diária. Ou no caso da inflação brasileira, neste caso o IPCA, que demonstra altas mais expressivas no início e no fim do ano.

Para esse tipo de análise usamos uma variação do modelo ARIMA, o chamado SARIMA. Que inclui a letra “S”, que representa a sazonalidade.

Pacotes

library(Quandl)

library(tidyverse)

library(gridExtra)

library(forecast)

library(scales)

library(tsibble)

library(feasts)

library(caret)

library(fpp3)

Pegando e visualizando os dados

Usaremos a tabela 433, do Banco Central, para pegar a série histórica do IPCA.

data = Quandl('BCB/433', type = "raw")

ggplot(data, aes(x=Date, y=Value))+
  geom_line()

Tratamento e visualização da amostra a ser utilizada

ipca <- 
  data %>%
  filter(Date >= '2007-01-01' & Date <= '2021-12-31') %>%
  mutate(Date = yearmonth(Date)) %>%
  as_tsibble(index = Date)

ipca %>%
  ggplot(aes(x=Date, y=Value))+
  geom_line()

Nos ateremos a analisar o período de 01/2007 até 12/2021. Assim, não pegaremos o período de hiperinflação (1994), nem o choque inflacionário após a vitória de Lula nas eleições (2003), pulando também a tendência de baixa no final de seu primeiro mandato (2006).

Visualização da sazonalidade

Podemos verificar a sazonalidade da série, como estamos falando de inflação, normalmente temos um aumento mais expressivo no começo e no fim do ano, devido o aumento do consumo.

ipca %>%
  gg_subseries(Value)

Análise de autocorrelação da amostra

acf1 = ggAcf(ipca$Value)
pacf1 = ggPacf(ipca$Value)
grid.arrange(acf1, pacf1)

Podemos observar pelo gráfico do ACF (autocorrelation function), que há uma correlação significativa no lag 1, que descreve no decorrer de alguns lags. Padrão característico de um termo auto-regressivo.

Podemos então olhar para o gráfico de PACF (partial autocorrelation function), para entender a ordem desse termo. No caso, temos uma correlação significativa apenas no primeiro lag, seguido por outras não significativas, habitual de um termo de ordem 1.

Criando o modelo e separando a amostra de treino e teste

ts_ipca = ts(ipca$Value, start = c(2007,01), frequency=12)

auto.arima(ts_ipca)
## Series: ts_ipca 
## ARIMA(1,0,0)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##         ar1    sma1    mean
##       0.574  0.1492  0.4770
## s.e.  0.061  0.0762  0.0502
## 
## sigma^2 = 0.06532:  log likelihood = -8.67
## AIC=25.33   AICc=25.56   BIC=38.1
model = Arima(ts_ipca, order = c(1,0,0), 
              seasonal = c(0,0,1))

intrain = createDataPartition(ts_ipca, p=0.9, list=F)
training = ts_ipca[intrain]
testing = ts_ipca[-intrain]

Temos a transformação dos dados em time series e depois o uso da função auto.arima para validar a nossa análise anterior, que aparenta estar correta.

Usamos essa coleta para o nosso modelo, e logo depois fazemos a separação da amostra em treino (90%) e teste (10%).

Análise de autocorrelação dos resíduos

fit_sarima <- auto.arima(training, seasonal = T)

ggtsdisplay(residuals(fit_sarima))

ggAcf(residuals(fit_sarima))

Box.test(residuals(fit_sarima), lag=24, 
         fitdf=length(coef(fit_sarima)),
         type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(fit_sarima)
## X-squared = 30.602, df = 20, p-value = 0.06067

A análise dos residuos, não mostra uma autocorrelação significativa visualmente para os resíduos. Demonstra um p-value relativamente alto, portanto, não rejeitamos a hipótese nula.

Diagnóstico

sarima_treino = Arima(training, order = c(1,0,0), seasonal=c(0,0,1))
sarima_teste = forecast::forecast(sarima_treino, h = length(testing))

diag = forecast::accuracy(sarima_teste,testing)
diag
##                         ME      RMSE       MAE       MPE     MAPE     MASE
## Training set  0.0001775496 0.2634840 0.2000025      -Inf      Inf 0.894878
## Test set     -0.0583417244 0.2853361 0.2400191 -20.18367 121.9462 1.073926
##                     ACF1
## Training set -0.04665686
## Test set              NA

Visualização da previsão

fsarima = forecast::forecast(model, h=12)

autoplot(fsarima)

Conclusão

O modelo tem um mean error razoável, comparado com o RMSE. Porém, provavelmente é um modelo pouco preciso fora da amostra, até devido ao aumento de variância do IPCA nos últimos anos, que dificulta a previsão dela baseado apenas em seus valores históricos.