Modelo ARIMA no R

Demonstração prática da aplicação de modelos ARIMA em séries macroeconômicas brasileiras usando a metodologia Box-Jenkins, desde a identificação da ordem do modelo até a previsão de valores futuros.
data-science
time-series
econometria
Author

Vinicius Oike

Published

July 1, 2025

O modelo

A teoria dos modelos ARMA(p,q) serve apenas para séries estiacionárias; na prática, contudo, é bastante frequente encontrar séries não-estacionárias. Séries com tendência de crescimento ou cuja variância aumenta ao longo do tempo são exemplos típicos de não-estacionaridade que encontra-se nos dados.

A ideia da modelagem ARIMA é de tirar diferenças na série original de forma a torná-la estacionária. Isto é, se tivermos uma série \(y_{t}\) que não é estacionária, vamos transformá-la em \(x_{t} \equiv y_{t} - y_{t-1}\) para torná-la estacionária.

Neste caso, a série \(x_{t}\) é chamada a primeira diferença de \(y_{t}\). Eventualmente será necessário tirar mais do que uma diferença para tornar a série estacionária.

Contexto

A análise de séries temporais ocupa um lugar central na macroeconomia moderna, fornecendo ferramentas essenciais para compreender a dinâmica das principais variáveis econômicas. Os modelos ARIMA (AutoRegressive Integrated Moving Average) representam uma das abordagens mais consolidadas e amplamente utilizadas para modelagem de séries temporais univariadas.

O ARIMA combina simplicidade e boa capacidade preditiva. Neste sentido, ele serve tanto como uma ferramenta central, dentro de uma análise macroeconômica, como também como benchmark para modelos mais sofisticados.

R

A modelagem de séries temporais no R é um tanto caótica. Isto acontece pela diversidade e complexidade de dados em séries de tempo. Via de regra, as funções-base do R servem muito bem para séries mensais, trimestrais e anuais. Quando se lida, por exemplo, com dados financeiros diários ou intra-diários costuma ser necessário usar pacotes especializados.

Com o tempo, as funções base do R se tornaram um tanto defasadas. Os gráficos do R foram revolucionados pelo pacote ggplot2 mas a sua sintaxe não funciona tão diretamente com séries de tempo; como resultado, alguns ajustes são necessários para compatibilizar os dois.

Neste post, tentei manter o código o mais simples possível usando somente o pacote forecast1, que é excelente para o tratamento de séries de tempo univariadas2. Abaixo segue a lista dos pacotes necessários para acompanhar este post3.

library(forecast)
library(urca)
library(qs)

Metodologia

Este post tem como objetivo principal demonstrar a aplicação prática dos modelos ARIMA em séries macroeconômicas brasileiras, explorando tanto a metodologia tradicional de Box-Jenkins quanto abordagens automatizadas mais recentes.

Box-Jenkins

A metodologia Box-Jenkins parte do princípio da parcimônia. Esta abordagem iterativa é composta por três etapas fundamentais: identificação, estimação e verificação.

A identificação da ordem do modelo é feita visualmente usando gráficos de autocorrelação e autocorrelação parcial.

Vários modelos concorrentes são estimados e escolhe-se o mais apropriado usando uma mistura de parcimônia, critérios de informação e testes estatísticos. Resumidamente,

  1. Usando um teste de raiz unitária, identifica se a série é não-estacionária. Se a série for estacionária pode-se usar um modelo ARMA(p, q).
  2. Se a série for não-estacionária, tira a primeira diferença da série e repete o processo até que a série se torne estacionária.
  3. Avalia-se a FAC e FACP para propor um modelo de “ordem máxima”.
  4. Estima-se vários modelos ARIMA(p,q) de ordens menores.
  5. Seleciona-se o melhor modelo segundo algum critério. A abordagem mais simples é escolher algum critério de informação como o BIC.

Aplicaremos esta metodologia completa à série de consumo das famílias brasileiras.

Macroeconomia: consumo no Brasil

Para acompanhar estes exemplos vamos usar as séries das contas nacionais brasileiras. Estas séries são trimestrais e dessazonalizadas. Além disso, os seus valores são indexados, em base-100, em 1995, o primeiro ano da série completa.

data pib consumo governo fbkf export import
1996-03-31 100.84 98.63 99.14 96.84 98.21 91.34
1996-07-01 100.36 100.88 100.74 98.94 95.94 99.91
1996-10-01 104.40 103.93 104.24 102.12 95.63 107.19
1996-12-31 103.28 109.17 88.58 106.81 100.45 118.61
1997-03-31 104.07 106.42 100.31 108.11 105.16 118.78
1997-07-01 105.10 106.92 99.91 109.03 111.29 123.62

Inicialmente, vamos usar somente a série de consumo.

# Série de consumo
cons = ts(macro$consumo, start = c(1996, 1), frequency = 4)

1. Estacionaridade

Grosso modo, uma série estacionária possui propriedades estatísticas que não variam no tempo. Isto garante que métricas usuais como média e variância continuem fazendo sentido.

Propriedades de séries estacionárias.

  1. Média constante.
  2. Variância finita e constante.
  3. Autocovariância que depende somente do tamanho do intervalo de tempo considerado.

Costuma ser fácil enxergar se uma série é estacionária ou não-estacionária. Visualmente, uma série não-estacionária cresce ou diminui ao longo do tempo; em alguns casos, ela possui períodos de oscilação (i.e. variância) mais/menos intensa.

Além disso, é muito mais comum encontrar séries que são não-estacionárias do que o contrário. Séries estacionárias oscilam de maneira regular em torno de uma média constante e são visualmente similares a um ruído branco. Já a maioria das séries de tempo que se encontra na prática são ou uma linha que sobe ou uma linha que desce.

Análise visual:

  1. Série aumenta ou diminui muito ao longo do tempo (i.e. parece ter uma tendência).
  2. Volatailidade da série cresce ou diminui ao longo do tempo.

O gráfico abaixo, feito com a função autoplot, mostra a série de consumo dessazonalizada. Como se vê, ela apresenta claros indícios de ser não-estacionária, pois a média dela não é constante no tempo.

# echo: false
autoplot(cons) +
  ggtitle("Consumo das famílias (dessazonalizado)") +
  xlab(NULL) +
  ylab("Index (100 = 1995)") +
  theme_ts

Autocorrelação

Pode-se verificar o gráfico de autocorrelação e autocorrelação parcial da série como um critério informal para avaliar a estacionaridade4. Em geral, se a FAC demora muito para decair, a série é não-estacionária. O gráfico abaixo mostra a FAC e FACP de cinco “ciclos” da série.

ggAcf(cons, lag.max = frequency(cons) * 5)
ggPacf(cons, lag.max = frequency(cons) * 5)

Teste de raiz unitária

Testes de raiz unitária podem ser bastante complexos. Tipicamente, aplica-se um teste ADF (Augmented Dickey-Fuller) mais geral (com constante e tendência) e vai-se afunilando até o ADF mais simples (sem constante e sem tendência). Enders (1995) apresenta uma boa exposição do assunto e oferece um fluxograma de como melhor aplicar os testes.

O grande problema do teste ADF é que ele tem baixo poder e queremos não-rejeitar (aceitar) a hipótese nula. O fluxograma de Enders visa minimizar o risco de se aceitar a hipótese de raiz unitária devido ao baixo poder do teste ADF. Alternativamente, pode-se usar o teste Philips-Perron (ur.pp).

Abaixo mostro como fazer o teste ADF, nas suas três versões, usando o critério BIC para seleção automática do número de defasagens. Neste caso, fica claro que a série possui uma ruiz unitária sem constante e sem tendência. A rigor, seria necessário aplicar os mesmos testes sobre a primeira diferença da série para verificar se há somente uma raiz unitária.

Os resultados podem ser melhor detalhados usando summary , i.e., summary(adf1).

ur.df(cons, type = "trend", lags = 12, selectlags = "BIC")
adf1 = ur.df(cons, type = "trend", lags = 12, selectlags = "BIC")
adf2 = ur.df(cons, type = "drift", lags = 12, selectlags = "BIC")
adf3 = ur.df(cons, type = "none", lags = 12, selectlags = "BIC")

2. Identificação

A etapa de identificação constitui o ponto de partida da metodologia Box-Jenkins, onde determinamos a estrutura mais apropriada do modelo ARIMA(p,d,q) através de análise gráfica.

Train e test

Uma prática comum em séries de tempo é remover um percentual das observações finais afim de testar a capacidade preditiva do modelo5. Esta abordagem também é bastente comum na literatura de Machine Learning, onde o objetivo costuma ser gerar previsões.

O código abaixo cria uma série train com cerca de 85% das observações. Vale notar que há maneiras mais sofisticadas de fazer esta divisão e veremos como fazer uma delas mais adiante.

a = 0.15
h = floor(a * length(cons))
train = head(cons, length(cons) - h)
test = tail(cons, h)

Autocorrelação

O gráfico abaixo mostra a FAC e FACP da primeira diferença. Agora, temos somente a primeira defasagem significativa em cada um dos gráficos. Assim, podemos supor que o modelo de ordem máximo é um ARIMA (1, 1, 1).

dcons = diff(train)

ggAcf(train)
ggPacf(train)

3. Estimação

Vamos estimar quatro modelos: ARIMA(1, 1, 1), ARIMA (1, 1, 0), ARIMA (0, 1, 1) e ARIMA (0, 1, 0). A interpretação dos coeficientes requer cuidado especial:

  • Coeficiente AR(1): mede a persistência da série; valores próximos a 1 indicam alta persistência.

  • Coeficiente MA(1): captura o efeito de choques aleatórios passados sobre o valor atual.

  • Significância estatística: avaliada através de testes t individuais.

Vale notar que algumas estatísticas comuns como o R² ou R² ajustado não fazem sentido nos modelos ARMA e ARIMA devido às características assintóticas dos estimadores. Até por isso, é comum omitir ambos.

m1 = Arima(train, order = c(0, 1, 0), include.drift = TRUE)
m2 = Arima(train, order = c(1, 1, 1))
m3 = Arima(train, order = c(1, 1, 0))
m4 = Arima(train, order = c(0, 1, 1))
m5 = Arima(train, order = c(1, 1, 1), include.drift = TRUE)

A saída do modelo ARIMA (1, 1, 1) indica a estimativa dos coeficientes junto com a estimativa da variância do erro e o cálculo dos critérios de informação. O modelo estimado tem a forma:

\[ \Delta y_{t} = 0.9449 \Delta y_{t-1} - 0.7277 \varepsilon_{t-1} \]

m2
Series: train 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.9469  -0.7277
s.e.  0.0490   0.1162

sigma^2 = 2.907:  log likelihood = -153.58
AIC=313.17   AICc=313.49   BIC=320.28

4. Diagnóstico de resíduos

Para verificar o ajuste do modelo aos dados fazemos um diagnóstico sobre seus resíduos. Idealmente, os resíduos $e_{t} = y_{t} - $ ou os resíduos normalizados devem se comportar como ruído branco. Os resíduos normalizados são definidos como:

\[ e_{t} = \frac{y_{t} - \hat{y_{t}}}{\sqrt{\sigma^2}} \]

onde \(\sigma^2\) é o valor estimado da variância dos resíduos. A distribuição de \(e_{t}\) deve ser i.i.d. com média zero e variância unitária.

Queremos verificar as seguintes propriedades nos resíduos:

  • Ausência de autocorrelação.
  • Média igual a zero.
  • Variância constante (homocedasticidade).
  • Normalidade.

Em termos de previsão, o mais importante é verificar se a média do resíduo é zero. Qualquer valor diferente de zero indica que a previsão do modelo será enviesada. Em termos econométricos, o mais importante é verificar a ausência de autocorrelação. A hipótese de homocedasticidade é importante para a validade dos intervalos de confiança das previsões, embora violações moderadas costumem ser toleradas na prática

Para selecionar os resíduos de um modelo usa-se residuals. O valor de \(\sigma^2\) já foi calculado pela função Arima então precisamos somente selecionar o valor.

u1 = residuals(m1) / sqrt(m1$sigma2)
u2 = residuals(m2) / sqrt(m2$sigma2)
p1 = autoplot(u1) + labs(subtitle = "ARIMA (0, 1, 0)")
p2 = autoplot(u2) + labs(subtitle = "ARIMA (1, 1, 1)")
panel = p1 / p2

panel &
  plot_annotation(
    title = "Resíduo dos modelos ARIMA (0, 1, 0) e ARIMA (1, 1, 1)"
  ) &
  theme_ts +
    theme(axis.title = element_blank())

Novamente, faz-se uma verificação tanto visual como estatística.

Autocorrelação serial

Visualmente, temos duas opções para verificar a autocorrelação serial nos resíduos: a FAC e um lag plot.

Da mesma forma como fizemos antes, podemos fazer a FAC da série. Nota-se que a primeira defasagem é significativa (i.e. não é zero). Já a FAC do modelo ARIMA (1, 1, 1) não apresenta termos significativos, sugerindo que o resíduo não apresenta autocorrelação.

ggAcf(u1, lag.max = 20)
ggAcf(u2, lag.max = 20)

Alternativamente, podemos fazer um gráfico de defasagens (lag plot). Em cada um dos quadrantes abaixo o resíduo é plotado contra seu valor defasado. Se não houver autocorrelação, os gráficos devem ser nuvens “aleatórias” de pontos.

gglagplot(u1, do.lines = FALSE, colour = FALSE, lags = 4) + theme_ts

gglagplot(u2, do.lines = FALSE, colour = FALSE, lags = 4) + theme_ts

Testes Estatísticos

Formalmente, há pelo menos duas opções para testar a presença de autocorrelação: o teste Ljung-Box (ou Portmanteau) e o teste Breusch-Godfrey (BG test). O primeiro é o mais comum e o código abaixo mostra como calculá-lo.

Ljung-Box

Note que, neste caso fitdf = 0 pois não estimamos termos p e q no modelo ARIMA. A hipótese nula do teste Ljung-Box é de que a autocorrelação conjunta até defasagem k é igual a zero. Neste sentido, queremos não-rejeitar (aceitar) a hipótese nula.

No caso abaixo, vemos que o p-valor do teste com oito defasagens é pequeno e menor que 0.05, indicando que há alguma autocorrelação no resíduo do modelo.

Box.test(u1, type = "Ljung-Box", fitdf = 0, lag = 8)

    Box-Ljung test

data:  u1
X-squared = 16.311, df = 8, p-value = 0.03814

Pode ser interessante testar vários valores diferentes no teste Ljung-Box. O código abaixo mostra como fazer um loop simples para calcular vários valores desta estatística.

lags = c(1, 4, 8, 12, 16, 20, 24)
tabela = matrix(nrow = length(lags), ncol = 3)

for (i in seq_along(lags)) {
  lag = lags[[i]]
  lbtest = Box.test(u1, type = "Ljung-Box", fitdf = 0, lag = lag)
  tabela[i, ] = c(lag, as.numeric(lbtest$statistic), as.numeric(lbtest$p.value))
}

A tabela abaixo mostra o resultado do teste Ljung-Box para diversas defasagens.

Resultados do teste Ljung-Box sobre os resíduos do modelo
Lag Estatística P-valor
1 7.53578 0.00605
4 11.46158 0.02184
8 16.31102 0.03814
12 17.96947 0.11662
16 21.18401 0.17155
20 24.20865 0.23340
24 30.11548 0.18095

Breusch-Godfrey

O teste Breusch-Godfrey é, essencialmente, uma regressão linear do resíudo contra seus valores defasados. Na sua forma mais simples o teste é dado como:

\[ \hat{u}_{t} = \alpha + \rho_{1}\hat{u}_{t-1} + \rho_{2}\hat{u}_{t-2} + \dots + \rho_{p}\hat{u}_{t-p} + v_{t} \]

onde \(\hat{u}_{t-1}\) é o valor estimado do resíduo, \(\alpha\) é uma constante, \(\rho_{i}\) é um coeficiente que mede a autocorrelação entre o resíduo e a sua i-ésima defasagem e \(v_{t}\) é um termo de erro. Intuitivamente, queremos que todos os valores de \(\rho_{i}\) sejam iguais a zero pois o resíduo não deve ter uma relação linear com seus valores defasados.

Infelizmente, a função lmtest::bgtest foi feita para funcionar com o output de modelos de regressão linear criados com a a função lm. Assim, a função não funciona diretamente com um output de Arima() e é preciso montar a regressão manualmente, o que pode ser um pouco chato.

No código abaixo uso a função ts.union para reunir todas as séries (funciona como um full_join). A função lmtest::bgtest aceita uma sintaxe de formula similar a da função lm do tipo y ~ x1 + x2.

A hipótese nula do teste é de que \(\rho_{i} = 0\), ou seja, novamente queremos não-rejeitar a hipótese nula. Como se abaixo, em ambos os casos temos evidência de autocorrelação.

# Monta um data.frame com o resíduo e suas defasagens
udata = ts.union(
  u = u1,
  ul1 = lag(u1, -1),
  ul2 = lag(u1, -2),
  ul3 = lag(u1, -3),
  ul4 = lag(u1, -4),
  dframe = TRUE
)
# Roda o BG test com uma defasagem
lmtest::bgtest(u ~ 1 + ul1, data = udata, fill = NA, type = "F")

    Breusch-Godfrey test for serial correlation of order up to 1

data:  u ~ 1 + ul1
LM test = 0.44503, df1 = 1, df2 = 75, p-value = 0.5068
# Roda o BG test com quatro defasagens
lmtest::bgtest(
  u ~ 1 + ul1 + ul2 + ul3 + ul4,
  data = udata,
  fill = NA,
  type = "F"
)

    Breusch-Godfrey test for serial correlation of order up to 1

data:  u ~ 1 + ul1 + ul2 + ul3 + ul4
LM test = 3.2901, df1 = 1, df2 = 69, p-value = 0.07405

Heterocedasticidade

A maneira mais simples de verificar heterocedasticidade é procurar pela presença de picos na série do quadrado do resíduo. O gráfico abaixo mostra os dois gráficos lado a lado. Vê-se que há picos de volatilidade em ambos os resíduos.

autoplot(u1^2)
autoplot(u2^2)
p1 = autoplot(u1^2) +
  labs(
    subtitle = "Quadrado dos resíduos do modelo ARIMA (0, 1, 0)",
    x = NULL,
    y = latex2exp::TeX("$u_{t}^2$")
  )

p2 = autoplot(u2^2) +
  labs(
    subtitle = "Quadrado dos resíduos do modelo ARIMA (1, 1, 1)",
    x = NULL,
    y = latex2exp::TeX("$u_{t}^2$")
  )

(p1 | p2) +
  plot_annotation(
    title = "Verificando homocedasticidade",
    subtitle = "Picos na série do quadrado dos resíduos indicam a presença de heterocedasticidade."
  ) &
  theme_ts

Arch test

O test ARCH (autoregressive conditional heteroskedasticity) é, conceitualmente, muito similar ao teste Breusch-Godfrey, visto anteriormente. A ideia é fazer uma regressão do quadrado do resíduo contra seus valores defasados para verificar a presença de heteroscedasticidade condicional. A hipótese nula do teste é de que não há efeitos ARCH presentes.

FinTS::ArchTest(u1, lags = 4)

    ARCH LM-test; Null hypothesis: no ARCH effects

data:  u1
Chi-squared = 11.216, df = 4, p-value = 0.02424

Normalidade dos resíduos

O teste Jarque-Bera avalia se os resíduos seguem uma distribuição normal. A normalidade dos resíduos é importante para a validade dos intervalos de confiança das previsões.

tseries::jarque.bera.test(u1)

    Jarque Bera Test

data:  u1
X-squared = 4.8301, df = 2, p-value = 0.08936

A função checkresiduals oferece uma visualização completa dos diagnósticos do modelo, combinando o gráfico da série dos resíduos, a FAC e um histograma com o resultado do teste Ljung-Box.

checkresiduals(m1) + theme_ts


    Ljung-Box test

data:  Residuals from ARIMA(0,1,0) with drift
Q* = 16.311, df = 8, p-value = 0.03814

Model df: 0.   Total lags used: 8
NULL

5. Comparação de Modelos

Após estimar os quatro modelos e verificar seus diagnósticos, o próximo passo é compará-los para selecionar o mais apropriado. A tabela abaixo mostra os critérios de informação para cada modelo estimado.

# Extrai informações dos modelos
tab_models = data.frame(
  Modelo = c(
    "ARIMA(0,1,0)",
    "ARIMA(1,1,1)",
    "ARIMA(1,1,0)",
    "ARIMA(0,1,1)",
    "ARIMA(1,1,1) com drift"
  ),
  AIC = c(m1$aic, m2$aic, m3$aic, m4$aic, m5$aic),
  AICc = c(m1$aicc, m2$aicc, m3$aicc, m4$aicc, m5$aicc),
  BIC = c(m1$bic, m2$bic, m3$bic, m4$bic, m5$bic),
  sigma2 = c(m1$sigma2, m2$sigma2, m3$sigma2, m4$sigma2, m5$sigma2)
)
Modelo AIC AICc BIC sigma2
ARIMA(0,1,0) 318.6831 318.8410 323.4220 3.18430
ARIMA(1,1,1) 313.1667 313.4867 320.2750 2.90749
ARIMA(1,1,0) 320.2546 320.4125 324.9935 3.23828
ARIMA(0,1,1) 325.5904 325.7483 330.3293 3.46862
ARIMA(1,1,1) com drift 313.3759 313.9165 322.8537 2.89641

Os critérios de informação penalizam modelos mais complexos, buscando encontrar o melhor equilíbrio entre ajuste aos dados e parcimônia. Menores valores indicam melhores modelos. Tanto o AIC quanto o BIC indicam que o modelo ARIMA(1,1,1) apresenta o melhor ajuste.

Vale notar que o modelo ARIMA(0,1,0), também conhecido como random walk, apresenta os piores critérios. Este modelo assume que a melhor previsão para o próximo período é simplesmente o valor atual, sem considerar a estrutura de dependência temporal da série.

Modelo final

Com base nos diagnósticos e critérios de informação, selecionamos o modelo ARIMA(1,1,1) como o mais apropriado para a série de consumo. O modelo final estimado é:

\[ \Delta y_{t} = 0.9449 \Delta y_{t-1} - 0.7277 \varepsilon_{t-1} \]

Este modelo indica que a primeira diferença do consumo possui alta persistência (coeficiente AR próximo de 1) e também incorpora choques passados através do componente MA.

6. Previsão

Com o modelo final estimado, podemos gerar previsões para períodos futuros. O código abaixo gera previsões para os próximos 8 trimestres (2 anos) usando o modelo ARIMA(1,1,1).

# Gera previsões
h = length(test)
fcst = forecast(m2, h = h)

O gráfico abaixo mostra a série histórica junto com as previsões e seus respectivos intervalos de confiança de 80% e 95%.

Os intervalos de confiança aumentam à medida que avançamos no horizonte de previsão, refletindo a crescente incerteza sobre valores futuros. Esta é uma característica comum de modelos de séries temporais.

Avaliando a precisão

Para avaliar a qualidade das previsões, podemos compará-las com os valores observados no conjunto de teste que separamos anteriormente.

# Calcula métricas de erro
accuracy(fcst, test)
                    ME     RMSE      MAE        MPE      MAPE      MASE
Training set 0.0923341 1.672860 1.253162 0.08911013 0.9504304 0.2512341
Test set     6.0885365 8.376132 6.584560 3.43589199 3.7244126 1.3200731
                  ACF1 Theil's U
Training set 0.0485215        NA
Test set     0.8172789  8.097087

A função accuracy calcula diversas métricas de erro de previsão. O RMSE (Root Mean Squared Error) e o MAE (Mean Absolute Error) medem o tamanho típico dos erros de previsão. O MAPE (Mean Absolute Percentage Error) expressa o erro em termos percentuais.

Note que as métricas calculadas sobre o conjunto de treino (Training set) são naturalmente melhores do que as calculadas sobre o conjunto de teste (Test set), pois o modelo foi estimado usando os dados de treino.

Podemos fazer os mesmos cálculos para todos os modelos usando lapply e definindo uma função que extrai as métricas de erro somente para o conjunto de teste.

models = list(
  "ARIMA(0,1,0)" = m1,
  "ARIMA (1,1,1)" = m2,
  "ARIMA (1,1,0)" = m3,
  "ARIMA (0,1,1)" = m4,
  "ARIMA (1,1,1) com drift" = m5
)

collect_measures = function(m) {
  fcast = forecast(m, h = length(test))
  measures = accuracy(fcast, test)
  measures = measures[2, c("ME", "RMSE", "MAE", "MAPE")]
  return(measures)
}

measures = lapply(models, collect_measures)
tab_measures = do.call(rbind, measures)
ME RMSE MAE MAPE
ARIMA(0,1,0) -7.38142 7.54294 7.38142 4.22250
ARIMA (1,1,1) 6.08854 8.37613 6.58456 3.72441
ARIMA (1,1,0) 0.35612 2.86907 2.63088 1.50421
ARIMA (0,1,1) -0.50112 2.83434 2.49317 1.43254
ARIMA (1,1,1) com drift -0.72171 1.71544 1.38625 0.80201

Note como, apesar do ARIMA (1,1,1) ter os melhores critérios de informação, os modelos ARIMA (0,1,1) e ARIMA (1,1,0) apresentaram melhores métricas de erro. O modelo ARIMA (1,1,1), sobretudo, as melhores métricas de erro.

Visualmente, podemos ver estes modelos tem uma melhor performance pois eles não seguem a trajetória de queda recente da série. Na verdade, ambos os modelos ARIMA (0,1,1) e ARIMA (1,1,0) geram uma previsão “naive” onde o último valor da série é repetido diversas vezes.

Vale notar, que nosso período de train termina no último trimestre de 2015, em meio à crise econômica. Nosso modelo ARIMA (1,1,1) tem ótima performance durante os primeiros 4 trimestres de 2016: ele continua prevendo a queda no consumo acuradamente. Contudo, dada a sua estrutura matemática, ele continua a gerar previsões de que o consumo vai continuar diminuindo. Já no modelo ARIMA (1,1,1) com drift as previsões voltam a aumentar em valor à medida que o tempo passa.

Code 
f1 = forecast(m2, h = length(test))
f2 = forecast(m3, h = length(test))
f3 = forecast(m5, h = length(test))

autoplot(train) +
  autolayer(test, series = "Obsevado") +
  autolayer(f1$mean, series = "ARIMA (1,1,1)") +
  autolayer(f2$mean, series = "ARIMA (1,1,0)") +
  autolayer(f3$mean, series = "ARIMA (1,1,1) com drift") +
  labs(
    title = "Previsão do Consumo das Famílias",
    subtitle = "Modelo ARIMA(1,1,1) com intervalos de confiança de 80% e 95%",
    x = NULL,
    y = "Index (100 = 1995)"
  ) +
  scale_color_manual(
    name = "",
    values = c("#2C6BB3", "#1abc9c", "#f39c12", "#34495e")
  ) +
  theme_ts

O gráfico abaixo faz uma previsão 12 anos à frente para enfatizar o comportamento de longo prazo destes modelos. Os modelos ARIMA (1, 1, 0) rapidamente converge para um valor constante. As previsões de longo prazo do ARIMA (1,1,1) com drift seguem uma tendência linear simples, enquanto o modelo sem drift converge lentamente para uma constante.

Code 
f1 = forecast(m2, h = 48)
f2 = forecast(m3, h = 48)
f3 = forecast(m5, h = 48)


autoplot(train) +
  autolayer(f1$mean, series = "ARIMA (1,1,1)") +
  autolayer(f2$mean, series = "ARIMA (1,1,0)") +
  autolayer(f3$mean, series = "ARIMA (1,1,1) com drift") +
  labs(
    title = "Previsão do Consumo das Famílias",
    subtitle = "Modelo ARIMA(1,1,1) com intervalos de confiança de 80% e 95%",
    x = NULL,
    y = "Index (100 = 1995)"
  ) +
  scale_color_manual(
    name = "",
    values = c("#2C6BB3", "#1abc9c", "#f39c12", "#34495e")
  ) +
  theme_ts

Esta discussão mostra a complexidade da previsão de séries de tempo e como seguir roteiros automatizados podem levar a erros. No nosso exemplo, a escolha do corte train/test levou a série ao meio de uma crise econômica, o que tem impactos grandes sobre a estimação dos modelos. Há maneiras de contornar isto usando, por exemplo, cross-validation para selecionar janelas de train/test variáveis, ou incorporando outros dados para ajudar a melhor lidar com o período de recessão.

Também fica evidente como os modelos ARIMA geram previsões de longo prazo questionáveis. Dada a estrutura matemática destes modelos, eles geram previsões que podem ser bastante irrealistas em horizontes grandes de tempo. Assim, tipicamente, deve-se usar modelos ARIMA como uma ferramenta de previsão de curto prazo.

7. Resumo

Neste post, demonstramos a aplicação completa da metodologia Box-Jenkins para estimar modelos ARIMA em séries temporais macroeconômicas brasileiras. Seguindo as etapas de identificação, estimação, verificação e previsão, modelamos a série de consumo das famílias brasileiras.

Principais Conclusões

  1. Estacionaridade: A série de consumo apresentou raiz unitária, exigindo uma diferenciação (d=1) para torná-la estacionária.

  2. Identificação: A análise das funções de autocorrelação (FAC e FACP) sugeriu um modelo ARIMA(1,1,1) como candidato de ordem máxima.

  3. Diagnóstico: A verificação dos resíduos através de testes estatísticos e análise gráfica confirmou que o modelo ARIMA(1,1,1) apresenta resíduos bem comportados, sem autocorrelação significativa.

  4. Comparação: Os critérios AIC e BIC confirmaram que o modelo ARIMA(1,1,1) oferece o melhor equilíbrio entre ajuste e parcimônia.

  5. Previsão: O modelo final foi usado para gerar previsões fora da amostra, com intervalos de confiança apropriados que refletem a incerteza crescente ao longo do horizonte de previsão. A previsão do modelo escolhido foi inicialmente boa, mas se deteriorou ao longo do intervalo. Modelos alternativos tiveram melhor performance preditiva.

Quando usar modelos ARIMA

Modelos ARIMA são apropriados quando:

  • Trabalha-se com séries temporais univariadas.
  • Há evidência de dependência temporal nos dados.
  • O objetivo é fazer previsões de curto a médio prazo.
  • Busca-se um modelo parcimonioso e fácil de estimar.

Limitações importantes:

  • Não incorpora informação de outras variáveis (para isso, considere modelos ARIMAX ou VAR).
  • A capacidade preditiva tende a deteriorar rapidamente em horizontes longos.
  • Quebras estruturais podem invalidar o modelo.

Alternativas e extensões

  • auto.arima(): O pacote forecast oferece a função auto.arima() que seleciona automaticamente a ordem do modelo usando critérios de informação
  • Modelos sazonais: Para séries com sazonalidade, considere modelos SARIMA
  • Modelos multivariados: Para incorporar múltiplas séries, considere modelos VAR ou VECM
  • Não-linearidade: Para capturar assimetrias e não-linearidades, considere modelos GARCH ou TAR

Posts Relacionados

Footnotes

  1. Para mais sobre o pacote forecast vale consultar o livro Forecasting: Principles and Practice.↩︎

  2. O pacote forecast começa a ficar um pouco limitado quando se lida com muitas séries de tempo (e.g. mais de 20).↩︎

  3. Alguns pacotes adicionais foram utilizados para melhorar as visualizações apresentadas no post. Para ver o código completo consulte meu GitHub.↩︎

  4. A rigor, a FAC e FACP são apenas definidas para séries estacionárias, então não faz sentido plotar a FAC de uma série não-estacionária.↩︎

  5. Caso seu interesse seja encontrar o melhor ajuste aos dados, pode fazer sentido usar a série inteira e não remover observações.↩︎