Modelo ARIMA

Muito utilizado na modelagem e previsões de séries temporais, o modelo ARIMA foi sistematizado em 1976 pelos estatísticos George Box e Gwilym Jenkins. Seu nome deriva do inglês autoregressive integrated moving average, que significa modelo autorregressivo integrado de média móvel. Para ver alguns conceitos mais básicos sobre Estatística, clique no link.

Gráfico de série temporal e modelo ARIMA (em azul) com sua previsão e intervalos de confiança - veja como plotar um gráfico desses no final do post.

Gráfico de série temporal e modelo ARIMA (em azul) com sua previsão e intervalos de confiança - veja como plotar um gráfico desses no final do post.

Em estatística, regressão é uma técnica que permite explorar e inferir a relação de uma variável dependente Y (variável de resposta) com uma (X) ou mais variáveis independentes específicas (variáveis explicatórias ou regressoras). Pode ser usada como um método descritivo da análise de dados (como, por exemplo, o ajustamento de curvas) sem serem necessárias quaisquer suposições acerca dos processos que permitiram gerar os dados.

O método de estimação mais amplamente utilizado é o método dos mínimos quadrados ordinários (MMQ ou MQO) - técnica de otimização matemática que procura encontrar o melhor ajuste para um conjunto de dados tentando minimizar a soma dos quadrados das diferenças entre o valor estimado e os dados observados (tais diferenças são chamadas resíduos). A curva resultante é denominada de regressão de Y para X, sendo Y avaliado a partir de X. Uma boa planilha para praticar esse ajuste está disponível no STOA da matéria de Física Experimental: ajuste_de_reta.xls.

Regressão múltipla é uma coleção de técnicas estatísticas para construir modelos que descrevem de maneira razoável relações entre duas ou mais variáveis explicativas de um determinado processo. Admite-se que "Y" é a variável resposta (dependente) e "X1", "X2", etc, são as variáveis preditoras (independentes), onde os valores esperados de Y são dados por uma combinação linear das variáveis preditoras. Uma relação não-linear (como uma relação polinomial, por exemplo), mesmo envolvendo apenas uma variável preditora e a variável resposta, pode ser facilmente tratada no âmbito duma regressão linear múltipla.

Para calcular previsões de séries temporais com um grande número de observações, é razoável aceitar o ajuste de modelos autorregressivos. Se o número de dados for pequeno, essa análise não é apropriada, pois a hipótese básica de independência dos resíduos é quase sempre violada, gerando estimadores inconsistentes. Nesse caso, os métodos de alisamento são mais utilizados. Alternativamente, pode-se utilizar a técnica de filtragem adaptativa, que não faz discriminação entre ruído e padrão básico - sua hipótese básica é que o padrão de comportamento da série possa ser representado por uma soma ponderada de observações passadas.

Suponha que a variável aleatória y é linearmente relacionada com seus próprios valores defasados. Na forma mais geral da análise de regressão, a variável dependente transformada pode ser explicada por variáveis independentes também transformadas:

math formula

A hipótese de trabalho por trás da transformação é equivalente à afirmativa de que existe algum valor de λ para o qual yλ tem distribuição (aproximadamente) normal com média dada pelo produto vetorial β'X e variância constante. Para λ = 1, a transformação deixa a variável praticamente inalterada, e para λ = 0 implica na transformação logarítmica da variável.

Como o parâmetro da transformação (λ) é desconhecido, deve-se estimá-lo a partir dos dados da amostra. Isso torna o modelo não-linear nesse parâmetro e exige a utilização do método de máxima verosimilhança para se obter estimadores com propriedades desejáveis. Após realizar a padronização da variável de resposta y é possível transformá-la, segundo Box–Cox, com diferentes valores de λ, ajustar as várias regressões por MQO, e escolher aquele valor do parâmetro para o qual a soma dos desvios quadráticos seja mínimo. Existem alguns pacotes computacionais que realizam essas estimativas escolhendo o λ ótimo de maneira automática e fornecendo a estatística necessária para realizar testes de significância da estimativa desse parâmetro de transformação.

Na análise de modelos paramétricos, existe o método de Box & Jenkins, que consiste em ajustar modelos autorregressivos integrados médias móveis (ARIMA) a um conjunto de dados. Sua estratégia para a construção do modelo é baseada em um ciclo iterativo (a escolha do modelo é baseada nos próprios dados), cujos estágios são:

  1. Especificação - uma classe geral de modelos é considerada para análise (ARIMA, no caso);
  2. Identificação de um modelo com base na análise de autocorrelações e autocorrelações parciais;
  3. Estimação dos parâmetros do modelo identificado (as estimativas preliminares encontradas na fase de identificação serão usadas como valores iniciais nesse procedimento);
  4. Verificação do modelo ajustado através de uma análise de resíduos (menor erro quadrático médio) para saber se este é adequado para os fins (no caso, de previsão).

Caso o modelo não seja adequado, o ciclo é repetido.

O principal argumento utilizado para a modelagem é uma série temporal univariada ("univariate time series"). Uma série temporal consiste de observações únicas (escalares) registradas sequencialmente durante incrementos iguais de tempo, enquanto que a palavra "univariada" implica que há somente uma variável dependente - ou seja, existe uma função y = f(t) tal que t é a variável independente e y a variável dependente (de t, que no caso representa o tempo). Quando os valores da série podem ser escritos através de uma função desse modo, diz-se que a série é determinística - caso envolva também um termo aleatório, a série é chamada estocástica.

Uma série estacionária (ou convergente) flutua em torno de uma mesma média ao longo do tempo (a sequência de suas somas parciais convergem), diferentemente da série não estacionária (ou divergente), onde ocorrem tendências ascendentes e descendentes. Suas propriedades não dependem do momento em que foi observada. Para tornar a série estacionária deve-se tomar diferenças quantas vezes for necessário, até atingir estacionariedade - são chamadas de séries não estacionárias homogêneas. Um modelo é estacionário se, para todo t, tem-se os valores de média, variância e covariâncias constantes. Por exemplo, quando a série oscila ao redor de um nível e salta para outro, basta uma diferença para torná-la estacionária; se a série oscila em uma direção e depois muda, são necessárias duas diferenças para ficar estacionária.

O modelo autorregressivo de ordem p é usado quando há autocorrelação entre as observações, ou seja, quando o valor de uma variável no período t depende do seu valor no período anterior (t-1) e de um termo aleatório. O modelo de média móvel de ordem q é usado quando há autocorrelação entre os resíduos - de modo geral, quanto maior o número de termos que compõe a média móvel, mais suavizado fica. Já o modelo autorregressivo de média móvel (ARMA) é usado quando há autocorrelação entre as observações e autocorrelação entre os resíduos. Se o modelo ARMA(p,q) é uma diferença da série temporal Zt, então essa série temporal é uma integral do modelo, daí temos que Zt segue um modelo autorregressivo-integrado-médias móveis (ARIMA).

Considere uma série temporal escrita da seguinte forma:

math formula, t = 1,...,N

math formula (se S_t independe de T_t) ou math formula (caso contrário)

A função f(t), chamada "sinal", é composta da adição/multiplicação da componente ciclo-tendência (Tt), que inclui as flutuações cíclicas de longo período que não podem ser detectadas com os dados disponíveis, e da componente sazonal (St), enquanto que at é o "ruído aleatório". Pode-se considerar duas classes de modelos: modelos de erro ou de regressão, onde o sinal f(t) é uma função do tempo completamente determinada de observações independentes e at é uma sequência aleatória independente de f(t), e modelos ARIMA. A hipótese de erros não correlacionados introduz limitações na validade dos modelos que seguem equações do tipo apresentado acima. Para descrever o comportamento de séries onde os erros observados são correlacionados e influenciam a evolução do processo, são utilizados os modelos ARIMA.

O modelo ARIMA é dado em função de p, q e d:

math formula

Onde:

  • p é o número de termos autorregressivos;
  • d é o número de diferenças;
  • q é o número de termos da média móvel;
  • L é o operador defasagem;
  • Φ é o polinômio ligado ao operador autorregressivo de ordem p;
  • θ é o polinômio ligado ao operador de média móveis de ordem q;
  • εt (ou at) é o ruído branco de média zero.

Para d = 0, tem-se o modelo ARMA(p, q), e no caso de também q = 0, tem-se o modelo AR(p) - o modelo ARIMA(0, 1, 0) é o passeio aleatório (random walk). A parte autorregressiva do modelo prediz o valor no momento t, considerando valores anteriores nas séries no momento t-1, t-2, etc. A média móvel usa valores de resíduos passados (as diferenças entre o valor real e o valor predito baseado no modelo no momento t).

O vídeo abaixo explica esse modelo com mais detalhes e exemplos:

Interpretação dos critérios

Os modelos ARIMA são mais eficientes em situações que envolvem um pequeno número de séries temporais - comprimento entre 50 (mínimo) e 100 observações (Montgomery, 1990). Na escolha do modelo mais apropriado, busca-se o que for mais parcimonioso, ou seja, que envolva o mínimo de parâmetros possíveis a serem estimados e que explique bem o comportamento da variável resposta.

A verossimilhança (ou likelihood, em inglês) é uma característica de algo que é semelhante à verdade, provável. Em estatística, a noção de verossimilhança é uma função da probabilidade condicional (probabilidade de um evento A sabendo que ocorreu um outro evento B). Deve-se escolher o valor do parâmetro desconhecido que maximiza a probabilidade de obter a amostra particular observada, ou seja, o valor que torna aquela amostra a "mais provável". Assim, ao obter a função de verossimilhança, que depende dos parâmetros desconhecidos e dos valores amostrais, deve-se maximizar essa função (derivar a função e igualar a zero para encontrar o ponto de máximo) ou o logaritmo de base "e" dela, o que for mais conveniente (Morettin & Bussab, 2004).

Existem estimadores (também chamados 'critérios') definidos com base no máximo da função de verossimilhança (MFV ou MLE, em inglês) que são muito utilizados para avaliar os modelos. Em problemas regulares, a estatística da razão de verossimilhanças (ou "Likelihood Ratio" - LR) converge para a distribuição qui-quadrado com m graus de liberdade, onde m é a diferença entre as dimensões dos espaços paramétricos sob as duas hipóteses testadas (Terra, 2013).

A função de verossimilhança é definida como a função de densidade de probabilidade conjunta dos dados observados, em função dos parâmetros do modelo (ARIMA, no caso, apresentados na descrição da equação do modelo). Os estimadores de máxima verossimilhança são definidos como aqueles valores dos parâmetros, com relação aos dados observados, que são os mais verossímeis (prováveis), e maximizam a função de verossimilhança (MFV).

Existem critérios definidos para avaliar os modelos. De modo geral, quanto menor o "information criteria (IC)", melhor o modelo, já que são proprocionais ao valor do "log likelihood" e, consequentemente, à variância. Além disso, ambos os critérios penalizam modelos com muitas variáveis. O sinal positivo ou negativo não tem um significado importante. Dentre os principais critérios baseados no MFV, estão:

* AIC (Critério de Informação de Akaike) - distância relativa esperada entre dois modelos probabilísticos: admite a existência de um modelo "real" (com um alto número de dimensões) que descreve os dados e tenta escolher outro dentre um grupo de modelos avaliados, o que minimiza a distância (ou divergência) de Kullback-Leibler (veja mais no link). Essa distância é uma medida de quanta informação é "perdida" ao tentar representar um conjunto de medidas utilizando uma base conhecida. O AIC é definido por:

math formula

Onde Lp é a função de máxima verossimilhança do modelo e p é o número de variáveis explicativas consideradas no modelo.

* AICc - correção da AIC para amostras de tamanho finito, cujo termo somado ao AIC depende do modelo estatístico utilizado. Supondo-se que o modelo é univariado, linear e tem resíduos normalmente distribuídos, a fórmula é:

math formula

Onde n indica o tamanho da amostra e k indica o número de parâmetros. O AIC foi desenvolvido sob o conceito de que, assintoticamente (quando o tamanho da amostra tende ao infinito), ele converge para o valor exato da divergência de Kullback-Leibler. No entanto, quando há um número finito de amostras, este estimador se torna polarizado. Com isto, o AIC não pode falhar em escolher um modelo mais parcimonioso, podendo até escolher o modelo de maior ordem entre todos os modelos comparados. Daí surgiu a necessidade da correção.

* BIC (Critério de Informação Bayesiano) - definido como a estatística que maximiza a probabilidade de se identificar o modelo dentre os avaliados, considerando-se a existência de um "modelo verdadeiro" que descreve a relação entre a variável dependente e as diversas variáveis explanatórias dentre os diversos modelos sob seleção. A penalidade devido ao número de termos é maior do que no AIC.

math formula

Onde Lp é a função de máxima verossimilhança do modelo e p é o número de variáveis explicativas consideradas no modelo.

O AIC é o melhor para predição, pois é assintoticamente equivalente à validação cruzada, enquanto que o BIC é melhor para explicação, uma vez que permite uma estimativa consistente do processo de geração de dados subjacente.

Já adiantando um pouco do próximo tópico, a função auto.arima() do R pode estimar modelos através de um algoritmo de ajuste mais eficiente, usando somas condicionais de erros quadráticos, para evitar excesso de tempo computacional para séries longas ou para modelos sazonais complexos. Quando este processo de estimativa está em uso, a função aproxima os critérios de informação para um modelo (porque o "log likelihood" não foi calculado). A heurística simples é usada para determinar se a estimativa das somas condicionais quadráticas deve ser usada, se o usuário não indicar qual abordagem deve ser usada. Heurística é um método/processo criado com o objetivo de encontrar soluções para um problema muito difícil de se resolver, que envolve a substituição destas questões difíceis por outras de resolução mais fácil a fim de encontrar respostas viáveis, ainda que imperfeitas. Nesse caso, pode surgir o seguinte aviso (warning) em R: "Unable to fit final model using maximum likelihood. AIC value approximated".

Programação em R

O software livre R contém pacotes muito utilizados para trabalhar com estatística - em particular para elaborar modelos ARIMA. As rotinas devem envolver, basicamente:

  1. Ler arquivos CSV com as variáveis observadas e a previsão de regressores (se for o caso) em função do tempo
  2. Análise de auto correlação
  3. Construir modelo ARIMA
  4. Previsão usando o modelo
  5. Gráfico dos resultados (acertar limites e linhas exibidas)

Além do próprio programa R, também deve ser instalado o pacote "forecast" - veja mais sobre linguagem R clicando no link. Esse pacote possui funções de previsão para séries de tempo e modelos lineares, como:

  • BoxCox() - retorna o reultado de uma transformação "Box-Cox": f(x)=(x^lambda-1)/lambda
  • BoxCox.lambda() - calcula o parâmetro lambda da transformação Box-Cox
  • Arima() - ajusta o modelo ARIMA à série de dados univariada (função com uma só variável)
  • auto.arima() - retorna o melhor ajuste ao modelo ARIMA
  • forecast() - função genérica para uma série temporal ou modelo de série temporal.

A função ts() (do inglês time serie) considera a série temporal com os dados igualmente espaçados, a variável tempo não precisa ser dada explicitamente. Desse modo, é necessário somente informal o início da série e sua frequência ("freq/frequency"). O formato é ts(vector, start=, end=, frequency=); veja alguns exemplos:

# Dados anuais - não se utiliza o parâmetro 'frequency'
ts(dados, start = 2001)
# Dados mensais - informar ano e mês
ts(dados, start = c(2001,06), frequency = 12)
#Dados trimestrais - informar ano e trimestre
ts(dados, start = c(2001,1), frequency = 4)
#Dados diários - informar ano e dia do ano
ts(dados, start = c(2001,152), frequency = 365)

A função auto.arima() do R usa uma variação do algoritmo desenvolvido por Hyndman e Khandakar que combina testes de raiz unitária, minimização de AICc e MFV para obter um modelo ARIMA. As etapas do algoritmo podem ser vistas no link. A função auto.arima retorna uma lista com: coef, sigma2, var.coef, mask, loglik, aic, arma, residuals, call, series, code, n.code, model (phi, theta, delta, Z, a, P, T, V, h, Pn), bic, aicc, xreg, x e lambda, respectivamente.

Dentre os parâmetros (opcionais) da função auto.arima, estão:

  • stationary: TRUE para restringir a pesquisa para modelos estacionários.
  • xreg: vetor/matriz de regressores ("regressors") externa. São pesos dados às observações da série temporal.
  • stepwise: TRUE fazer a seleção stepwise (default). Caso contrário, ele procura sobre todos os modelos - pode ser muito lenta, especialmente para os modelos sazonais.
  • trace: TRUE para imprimir os modelos considerados
  • approx/approximation: TRUE para estimativa através de somas condicionais quadráticas e os critérios de informação utilizados para seleção de modelos são aproximadas. O modelo final ainda é calculado usando estimativa máxima verossimilhança. A aproximação deve ser usado para a série de longa data ou um período muito sazonal para evitar tempos excessivos de computação.
  • seasonal: Se FALSE, restringe a pesquisa para modelos não- sazonais.
  • ic: informa o critério de seleção a ser utilizado (por default, ic="aicc")

Depois de uma série temporal ser "estacionarizada" por diferenciação, o passo seguinte na montagem de um modelo ARIMA é determinar se os termos autorregressivos (AR) ou de médias móveis (MA) são necessários para corrigir qualquer autocorrelação que permaneça na série diferenciada. O R tenta algumas diferentes combinações de termos e vê o que funciona melhor (com menor AICc). Também é possível olhar para a função de autocorrelação (ACF) e autocorrelação parcial (FACP) da série diferenciada, e assim identificar o número de termos AR e/ou MA que são necessários. No eixo das abcissas estão os "lags" (ou intervalos de tempo). Os gráficos podem ser vistos através do comando "tsdisplay(diff(x,difference))", onde "x" é a série temporal e "difference" é a ordem de diferenciação. Uma interpretação prática dessas funções pode ser vista no meio do artigo ARIMA: Como evitar a mentalidade de rebanho ao analisar os dados de séries temporais. Algumas dicas para escolher o número de termos (regras para identificar modelo ARIMA) e para comparar os resultados de diferentes modelos pode ser visto nos links.

Saída da função diff com gráficos de diferença e de autocorrelação (mesma série temporal exibida na primeira figura do post)

Saída da função diff com gráficos de diferença e de autocorrelação (mesma série temporal exibida na primeira figura do post)

Uma análise de resíduos pode ser realizada imprimindo-se os gráficos de resíduos padronizados, ACF dos resíduos e "p-value" de Ljung-Box através da função "tsdiag(modelo)", onde o parâmetro de entrada deve ser a saída da função auto.arima(). O modelo deve apresentar os resíduos estacionários, com média zero e variância constante - se o gráfico dos resíduos mostra uma tendência sistemática positiva ou negativa, isso significa que uma outra função (não linear) deve ser escolhida. O "valor-p" (ou p-value) é a probabilidade de se obter uma estatística de teste igual ou mais extrema que aquela observada em uma amostra, ou seja, um valor-p pequeno significa que a probabilidade de obter um valor da estatística de teste como o observado é muito improvável, levando assim à rejeição da hipótese nula (de que não existe relação entre dois fenômenos medidos). Veja mais sobre o p-valor no post sobre Testes de Hipóteses.

Análise de resíduos de modelo ARIMA (mesma série temporal exibida na primeira figura do post)

Análise de resíduos de modelo ARIMA (mesma série temporal exibida na primeira figura do post)

A função auto.arima retorna o modelo escolhido (conforme os valores de p, d e q), os termos autorregressivos (arX) e de média móvel (maX), de "drift" (termo constante em modelo que tem d > 0 devido à sazonalidade), "intercept" (média do termo (1-L)dXt, próximo à média da série histórica) e dos regressores (se for o caso), assim como o erro/desvio padrão (s.e. - "standard error") de cada valor apresentado. Como o desvio padrão mostra o quanto de variação ou "dispersão" existe em relação à média (ou valor esperado), se todos os pontos estiverem exatamente sobre a linha de regressão, então regressão explicaria toda a variação, ou seja, s.e. seria igual a zero. A função também imprime os critérios de informação AIC, AICc e BIC, assim como o log do MFV (log likelihood) e a variância (sigma^2), que é definida como o quadrado do desvio padrão. Veja um exemplo (diferente dos dados usados nos gráficos):

Series: exemplo
Series: resid 
ARIMA(1,1,1)(1,0,1)[12]

Coefficients:
         ar1      ma1    sar1     sma1    prec
      0.5488  -0.9031  0.8937  -0.7234  0.0029
s.e.  0.1212   0.0694  0.1037   0.1601  0.0006

sigma^2 estimated as 0.1234:  log likelihood=-45.31
AIC=102.62   AICc=103.37   BIC=119.3

Nesse exemplo, houve uma derivação, um de autorregressão e um de médias móveis. O segundo conjunto de parênteses refere-se ao modelo desenvolvido para a parte sazonal, cujos coeficientes recebem um "s" na frente. O valor dentro dos colchetes é o número de períodos por temporada. Se a soma dos coeficientes de AR é próxima de 1, isso que mostra que os parâmetros estão perto da borda da região de estacionaridade.

O número de linhas dos regressores deve ser o mesmo da série temporal. O modelo criado por arima() e auto.arima() com um argumento xreg é uma regressão linear com erros ARMA:

math formula math formula

Onde:

  • yt é o valor da projeção no instante t
  • xt é a covariância no instante t e β é seu coeficiente;
  • Φ é o polinômio ligado ao operador autorregressivo de ordem p;
  • θ é o polinômio ligado ao operador de média móveis de ordem q;
  • zt é o ruído branco no instante t.

Mais informações nesse link.

Outro pacote em R bom para esse tipo de estatística é o ASTSA. Nele, é possível plotar gráficos de ACF, PACF, p-value, possui filtros, técnicas e várias séries de dados para serem usadas como exemplo.

Previsão

Do mesmo pacote do R, existe a função "forecast". Quando é feita a previsão com um modelo ARIMA sem regressores, ele simplesmente usa os valores passados de sua série histórica para prever valores futuros. Nesse caso, você pode especificar o seu horizonte "h", que limita os períodos de previsão.

Ao usar regressores para construir um modelo ARIMA, é preciso incluir valores futuros dos regressores - o parâmetro "h" é ignorado. Por exemplo, se você usou temperatura como regressor para prever a incidência de uma doença, então você precisaria de valores futuros de temperatura para prever a incidência da doença.

A função forecast retorna uma lista composta de: method, model, level, mean, lower, upper, x, xname, fitted e residuals, respectivamente.

A variável de entrada level é opcional e representa os níveis/intervalo de confiança -  intervalo estimado onde a média tem uma dada probabilidade de ocorrer, que no caso do exemplo abaixo são as opções 65% (cinza claro) e 80% (cinza escuro). Sua interpretação seria algo como "existe 80% de probabilidade de que o cálculo do intervalo de confiança de alguns experimentos futuros envolvam o valor verdadeiro do parâmetro da população”. Veja mais no blog do Rob Hyndman.

O gráfico da série temporal de valores observados junto com os valores modelados para esse período e a previsão com os respectivos intervalos de confiança pode ser plotado através do R com as seguintes linhas:

# Ler dados observados de arquivo CSV
dados.obs = read.csv(filepath, header = TRUE)
# Criar série temporal com uma coluna dos dados observados
dados.ts = ts(dados.obs$valor, frequency = 12, start = c(start.year.hist, start.month.hist))
# Fazer modelo ARIMA
model = auto.arima(dados.ts)
# Projeção 
forec = forecast(model, level = c(65,80))
# Impressão de gráfico com série temporal e projeção com intervalos de confiança
plot(forec)
# Incluir no gráfico os valores do modelo para dados observados em azul
lines(fitted(model),col="blue")

O código para fazer o gráfico em si corresponde às duas últimas linhas não comentadas, sendo que as linhas anteriores correspondem a leitura de dados, modelagem e previsão. Note que "forec" é a saída da função "forecast".

Fontes

ps: Diferençar ou diferenciar? As duas palavras existem na língua portuguesa. Significam "distinguir, discriminar" ou também "aplicar cálculo diferencial (calcular as diferenças de uma função ou sequência)". O verbo "diferençar" é formado a partir de derivação sufixal (quando acrescenta-se um sufixo a uma palavra já existente, alterando o sentido da mesma). Já o verbo diferenciar poderá ter sua origem na palavra em francês différencier.
Compartilhe o link desse texto, mas se for copiar algum trecho, cite a fonte. Valorize nosso trabalho.
Mais informações na licença de uso do site.

5 Pingbacks/Trackbacks