Previsão econômica com modelos ARIMA no Python

A família de modelos ARIMA constitui uma abordagem de previsão para séries temporais buscando descrever as autocorrelações dos dados, que é uma característica comum em dados econômicos. Esses modelos se popularizaram muito em 1970 com o trabalho de Box e Jenkins, que é a referência base sobre o assunto, e evoluíram para abordagens modernas de automatização com o trabalho de Hyndman e Khandakar em 2008, que contribuíram com a popularização de software open source para uso e aplicações destes modelos.

Existem muitos modelos ARIMA e aqui apresentaremos uma visão geral focando em 4 variações, nesta ordem:

  • Modelo AR
  • Modelo MA
  • Modelo ARIMA
  • Modelo SARIMA

O foco será desenvolver a intuição sobre os modelos para entender a mecânica de funcionamento, as representações em equações e a aplicação aos dados, com finalidade de previsão de séries temporais econômicas do Brasil. Exploraremos recursos visuais (gráficos) para facilitar a compreensão e demonstraremos aplicações dos modelos com exemplos práticos em Python.

A referência técnica utilizada é o livro-texto de previsão de séries temporais intitulado “Forecasting: Principles and Practice” de Hyndman e Athanasopoulos (2021).

Para aprender mais e ter acesso a códigos confira o curso de Modelagem e Previsão usando Python ou comece do zero em análise de dados com a formação Do Zero à Análise de Dados com Python.

Modelos auto-regressivos (AR)

Em regressão linear  buscamos modelar a relação entre uma variável dependente y_t e outras variáveis independentes x_t, para prever y_t usando uma combinação linear de x_t. De forma similar, em modelos auto-regressivos (AR) estamos interessados em prever y_t com base em uma combinação linear de valores passados de y_t. Em outras palavras, o modelo AR é uma regressão da variável contra ela mesmo.

intuição por trás desse modelo é que as observações passadas da série temporal possuem uma influência ou são preditoras dos valores futuros da série.

Por exemplo:

  • A taxa de juros Selic para o mês que vem será, provavelmente, um valor muito próximo ou igual ao valor do mês corrente.
  • A taxa de inflação, medida pelo IPCA, apresenta inércia (um termo do econômes que significa, de forma simplificada, o aumento dos preços correntes com base nos preços passados) devido a vários fatores, como a indexação de contratos de aluguel.

A forma comumente utilizada para representação desse modelo é:

    \[y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \varepsilon_{t}\]

onde:

  • y_t é a série temporal alvo de previsão.
  • c é a constante ou intercepto do modelo, relacionado com a média de y_t porém não é a média em si.
  • \phi_p é o p coeficiente auto-regressivo a ser estimado, geralmente restrito em intervalos que geram modelos estacionários.
  • \epsilon_t é o termo de erro do modelo, um ruído branco.

É comum se referir a esse modelo usando o acrônimo AR(p), ou seja, um modelo auto-regressivo de ordem p, onde p é a ordem de defasagens no modelo.

A ilustração abaixo mostra a linha de ajuste de um modelo AR(1) para diferentes valores de \phi_1, no intervalo entre -1 e +1, usando como exemplo a série temporal da taxa de inflação medida pelo IPCA (IBGE):

A linha preta representa os dados observados e a linha azul representa o ajuste do modelo.

estimação dos coeficientes c, \phi_1, \dots, \phi_p do modelo pode ser feita pela máxima verossimilhança (MLE), que é uma técnica para encontrar os valores para os coeficientes que maximizam a probabilidade de obter os dados que foram observados. A técnica é bastante similar à minimização da soma dos quadrados dos resíduos na regressão linear.

Exemplo de modelo AR(p) com código: taxa de inflação - IPCA

Agora vamos mostrar como utilizar o modelo AR(p) para previsão, usando como exemplo os dados exibidos acima da taxa de inflação medida pelo IPCA (IBGE). Abaixo reportamos (a) o sumário estatístico do modelo, (b) as métricas de acurácia de treino/teste e (c) o gráfico de valores observados com ajuste dentro da amostra e previsão fora da amostra.

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  211
Model:               SARIMAX(1, 0, 0)   Log Likelihood                 -24.881
Date:                Mon, 05 Feb 2024   AIC                             55.763
Time:                        08:06:41   BIC                             65.818
Sample:                    01-01-2005   HQIC                            59.828
                         - 07-01-2022                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.1906      0.037      5.217      0.000       0.119       0.262
ar.L1          0.5864      0.052     11.185      0.000       0.484       0.689
sigma2         0.0740      0.005     14.893      0.000       0.064       0.084
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                50.09
Prob(Q):                              0.90   Prob(JB):                         0.00
Heteroskedasticity (H):               4.20   Skew:                            -0.30
Prob(H) (two-sided):                  0.00   Kurtosis:                         5.31
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
RMSE de treino: 0.2720257157158568
RMSE de teste: 0.2990384885047852

O código estima automaticamente a ordem p do modelo com base no menor critério de informação AICc. O modelo estimado é um AR(1):

    \[y_t = 0.1906 + 0.5864 y_{t-1} + \varepsilon_{t}\]

onde \varepsilon_t é um ruído branco com desvio padrão de 0.272= \sqrt{0.0740}. Note que as previsões geradas pelo modelo são relativamente boas pela simplicidade do modelo e captam, a olho nu, o retorno à média dos dados.

A linha preta no gráfico mostra os dados observados, enquanto que as linhas vermelha e azul mostram o ajuste dentro da amostra e a previsão fora da amostra, respectivamente.

Note que avaliar pontos de previsão sem levar em consideração o intervalo de confiança pode esconder uma grande incerteza e levar a conclusões errôneas.

Modelos de média móvel (MA)

Em modelos de média móvel (MA) busca-se prever os valores de uma série temporal  y_t com base em uma combinação linear de erros de previsão \varepsilon_t passados de y_t. Em outras palavras, é um modelo de regressão múltipla que usa os erros passados como variáveis independentes, em vez de outras variáveis observáveis.

intuição por trás deste modelo é a suposição de que a variável y_t reage a choques aleatórios com uma defasagem (delay), dessa forma, o choque defasado \varepsilon_t é naturalmente uma variável a se colocar no modelo.

Por exemplo:

  • Imagine que o retorno diário de uma ação negociada na bolsa possa ser explicado por uma tendência e pelos efeitos ponderados de choques que aconteceram nos dias anteriores (acontecimentos políticos, desempenho econômico, decisões de investimento, etc.). É razoável supor que o choque do dia t-1 terá um efeito importante no retorno da ação no dia t, além de efeitos de choques que podem acontecer no próprio dia t. Extrapolando essa intuição, trata-se de um processo que pode ser modelado por um MA.

A forma comumente utilizada para representação desse modelo é:

    \[y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q}\]

onde

  • \varepsilon_t é o termo de erro do modelo (choque), um ruído branco.
  • \theta_q é o q coeficiente de média móvel a ser estimado, geralmente restrito em intervalos que geram modelos invertíveis.

É comum se referir a esse modelo usando o acrônimo MA(q), ou seja, um modelo de média móvel de ordem q, onde q é a ordem de defasagens no modelo.

A ilustração abaixo mostra a linha de ajuste de um modelo MA(1) para diferentes valores de \theta_1, no intervalo entre -0.99 e +0.99, usando como exemplo a série temporal dos retornos diários simples do Ibovespa (Yahoo Finance):

A linha preta representa os dados observados e a linha azul representa o ajuste do modelo.

estimação dos coeficientes c, \theta_1, \dots, \theta_p do modelo pode ser feita pela máxima verossimilhança (MLE), que é uma técnica para encontrar os valores para os coeficientes que maximizam a probabilidade de obter os dados que foram observados. A técnica é bastante similar à minimização da soma dos quadrados dos resíduos na regressão linear.

Exemplo de modelo MA(q) com código: taxa de inflação - IPCA

Agora vamos mostrar como utilizar o modelo MA(q) para previsão, usando como exemplo os dados exibidos acima da taxa de inflação medida pelo IPCA (IBGE). Abaixo reportamos (a) o sumário estatístico do modelo, (b) as métricas de acurácia de treino/teste e (c) o gráfico de valores observados com ajuste dentro da amostra e previsão fora da amostra.

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  211
Model:               SARIMAX(0, 0, 4)   Log Likelihood                 -21.730
Date:                Mon, 05 Feb 2024   AIC                             55.459
Time:                        08:08:39   BIC                             75.570
Sample:                    01-01-2005   HQIC                            63.588
                         - 07-01-2022                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.4596      0.046      9.952      0.000       0.369       0.550
ma.L1          0.5910      0.056     10.560      0.000       0.481       0.701
ma.L2          0.3337      0.068      4.872      0.000       0.199       0.468
ma.L3          0.3297      0.067      4.925      0.000       0.199       0.461
ma.L4          0.1147      0.063      1.812      0.070      -0.009       0.239
sigma2         0.0718      0.005     13.745      0.000       0.062       0.082
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):                37.19
Prob(Q):                              0.98   Prob(JB):                         0.00
Heteroskedasticity (H):               4.11   Skew:                            -0.13
Prob(H) (two-sided):                  0.00   Kurtosis:                         5.04
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
RMSE de treino: 0.2679476899168494
RMSE de teste: 0.30324435848274084