Modelos ARIMA e extensões (SARIMA, ARIMAX, ARIMAX com sazonalidade)

Metodologia de Box & Jenkins para Identificação: FAC e FACP

(a) Serve para verificar se existe necessidade de uma transformação na série original com o objetivo de estabilizar a sua média, variância e autocovariâncias.

(b) Tomar diferenças da série tantas vezes quantas necessárias para obter uma série estacionária (atendendo assim as hipóteses e permitindo que façamos interpretações/inferências com propriedades conhecidas), de

modo que o processo seja reduzido a um ARMA(p,q). O número de diferenças, d, necessárias para que o processo se torne estacionário é alcançado quando as funções de autocorrelação e autocorrelação parcial

(FAC e FACP) amostral de Wt= decrescem rapidamente para zero. Neste estágio, a utilização de um teste para verificar a existência de raízes unitárias pode ser útil (o teste de Dickey-Fuller aumentado (ADF) é o mais utilizado).

(c) Identificar o processo ARMA(p,q) resultante por meio da análise das autocorrelações e autocorrelações parciais estimadas, cujos componentes devem imitar os comportamentos das respectivas quantidades teóricas. Uma forma alternativa de identificação é combinar as funções FAC e FACP com um penalizador (AIC, BIC ou HQ).

Usualmente, utilizamos um modelo “máximo” (considerando os lags significativos das funções FAC e FACP) e testamos combinações de ordem aninhada (menores) até encontrar aquele que minimiza os critérios de informação.

AIC seleciona melhor quando a amostra é pequena (100 observações em um processo estocástico é considerado muitas vezes como pequeno, como regra de bolso). BIC seleciona o verdadeiro modelo assintoticamente

(quando ) e costuma-se utilizar como regra de bolso . HQ é um penalizador intermediário.

Diagnósticos: assim como em qualquer análise estatística, após estimar e selecionar o modelo precisamos verificar se não conseguimos rejeitar (refutar) as hipóteses nos quais eles se assentam.

Teste de autocorrelação serial (estatística Q de Ljung-Box).

Modelos sazonais:

É comum em séries econômicas termos sazonalidade. Quando a sazonalidade é determinística (exemplo: consumo de sorvete aumenta no verão) então podemos captar este efeito utilizando uma variável dummy.

A sazonalidade, porém, pode não ser determinística. Os modelos SARIMA servem para captar este efeito quando ele é multiplicativo (a amplitude “sazonal” vai aumentando ao longo do tempo). Quando ele for aditivo podemos utilizar filtros pré-ARIMA (Holt Winters com sazonalidade, por exemplo). Os filtros mais conhecidos foram criados na década de 60, antes dos procedimentos de Box & Jenkins.

Cabe ressaltar que a modelagem de Box & Jenkins sofre muitas críticas por não ter um modelo teórico por trás (costuma-se chamar o procedimento de ateórico). O procedimento envolve uma decomposição numérica dos dados visando encontrar padrões de correlação teóricos (os mais similares dentro da classe de modelos utilizados) com base nos dados históricos (passados) sem a inclusão de outras variáveis (a série se auto-explica – mais adiante comentamos que isto não é um absurdo como é demonstrado na decomposição de Wolf). O método é popular, pois apresenta bons resultados quando validamos as suas hipóteses (se a série for estacionária o procedimento sempre nos trará resultados razoáveis, mesmo que não sejam os melhores possíveis).

Os modelos (S)ARIMAX são aqueles que incorporam regressores exógenos ao modelo (por exemplo, PIB e Inflação) constituindo assim uma combinação entre um modelo teórico e um ateórico.

Estacionariedade fraca

Um processo estocástico (ou uma série temporal) é dito fracamente

estacionário se a sua variância for finita ( ) e :

(1) A função valor médio, (uma constante) , para

todo

(2) (a autocovariância não depende de t, apenas do lag (distância) h entre dois pontos t e t+h.

Ljung-Box é um teste do tipo Portmanteau (teste no qual a hipótese nula está bem definida e a alternativa é mais “vagamente” especificada) Ao invés de testar a autocorrelação em cada lag distinto, ele testa a aleatoriedade "global" com base em um número de defasagens. Por esta razão, é muitas vezes referido como um teste de "Portmanteau”.

A hipótese nula para este teste é que as primeiras m autocorrelações são

conjuntamente zero, isto é, . A alternativa é que pelo uma das autocorrelações até m seja diferente de zero (não especifica qual). Em outras palavras, a hipótese nula é que as autocorrelações na população de onde a amostra foi retirada é zero, logo quaisquer correlações amostrais observadas são fruto da aleatoriedade contida na amostra. A alternativa é de que os dados exibem correlação serial.

A escolha de m afeta a performance do teste. Se T é o número de observações da série temporal series, Tsay recomenda utilizar ln(T) (devido ao poder do teste).

[2] Tsay, R. S. Analysis of Financial Time Series. 3rd ed. Hoboken, NJ: John Wiley & Sons, Inc., 2010.

Se aceitarmos a hipótese de que as autocorrelações dos resíduos são iguais a zero, então consideramos que o modelo é razoável (não refutamos a hipótese de que ele esteja bem ajustado com base na amostra)

Na prática, podemos adotar os seguintes passos em busca de um modelo mais parcimonioso:

1) Carregue os dados;

2) Fazer o gráfico da série temporal para checar anomalias (observações não usuais) e observar padrões como em qualquer análise de dados e o correlograma amostral (gráfico da FAC e da FACP);

3) Diferencie e/ou transforme os dados (se for necessário; use uma transformação da classe de Box-Cox se precisar estabilizar a variância). Caso não seja necessário fazer qualquer alteração, passe para (5).

4) Faça o gráfico (correlograma) das séries diferenciadas/transformadas;

5) Identifique as ordens de dependência (usando FAC e FACP ou então utilizando um algoritmo automático (auto.arima() no R, por exemplo). Caso você não utilize um algoritmo automático, siga os seguintes procedimentos: (a) se necessário diferencie novamente série até que ela pareça ser estacionária, (b) faça um teste de raiz unitária se você não estiver certo disto;

6) Reveja as funções (FAC e FACP) da série transformada e tente novamente determinar as ordens (possíveis candidatos “máximos”). Especifique o modelo “máximo” com base na FAC e FACP;

7) Escolha o modelo usando AIC, BIC ou AICc , isto é, a melhor combinação de ordem igual ou menor a do modelo “máximo” que minimize os critérios de informação desejados;

8) Faça um correlograma dos resíduos e diagnósticos (em especial, um teste de não-autocorrelação usando a estatística Q de Ljung-Box, por exemplo). Caso os resíduos não passem por este teste, passamos para o segundo melhor modelo encontrado em (7) e testamos novamente (e assim por diante até encontrarmos um modelo razoável para os dados que temos em mão – que não conseguimos refutar que tenham sido gerados pela nossa população estacionária hipotética).

9) Se desconfiarmos que o processo apresenta sinais de heterocedasticidade condicional, podemos fazer um teste ARCH. O modelos da classe ARIMA servem para modelar a média condicional.

10)Os resíduos parecem um ruído branco? Se não volte para (7) e escolha o segundo melhor modelo de acordo com o(s) critério(s) desejado(s),

11) Caso o modelo escolhido passe pelos diagnósticos, fazemos previsões e validações do modelo (checar a performance preditiva do modelo).

Graficamente, o procedimento pode ser descrito assim:

A previsão costuma ser avaliada através de uma gama de medidas (como comentamos na parte de regressão

linear).

Vamos colocar um pouco de notação aqui para explicar melhor os conceitos. Seja yt o valor a variável de interesse no período t e ft a previsão de yt. Defina o erro de previsão por et= yt-ft. Dada uma série temporal com T observações e suas previsões nós podemos construir medidas globais de precisão de suas previsões. As mais conhecidas são erro médio (ME), erro quadrático médio (MSE), a raiz do erro quadrático médio (RMSE), o erro absoluto médio (MAE), o erro percentual médio (MPE) e o erro absoluto percentual médio (MAPE). Estas medidas são

assim definidas:

Outra estatística relevante é a U de Theil (já comentamos antes) definida por:

A U de Theil pode ser definida como a razão entre RMSE do modelo de previsão proposto e a RMSE do modelo de previsão ingênuo yt+1=yt para todo t. Além disso, Theil (1966) também propôs uma decomposição do erro quadrático médio (MSE) que pode ser útil para avaliar um conjunto de previsões. Ele mostrou que o erro quadrático médio pode ser particionado em três componentes não negativos:

onde e são as médias amostrais das previsões e das observações, sf

e sy são seus respectivos erros-padrão (usando T no denominador), e r é a correlação amostral entre y e f. Dividindo os termos acima por MSE, obtemos:

Theil chamou os três termos de proporção do viés (UM), proporção da regressão (UR) e proporção residual (UD). Se y e f representam observações dentro da amostra da variável dependente e suas previsões, então esperaríamos, se o modelo está prevendo bem, que UM e UR fossem próximos de zero e que todo o erro quadrático médio correspondesse ao termo não sistemático residual.

Explicando de outra maneira: se planejamos ajustar um modelo para previsão, é uma boa prática avaliar a performance preditiva do modelo. Modelos que se ajustam bem dentro da amostra, não necessariamente fazem boas previsões. Por exemplo, um sobreajustamento (overfitting) , fará com que o ajuste seja muito bom dentro da amostra, mas provavelmente nos levará a uma pobre performance preditiva.

Ao avaliar a performance preditiva, é importante não utilizar os dados ”duas vezes”, isto é, os dados que você utiliza para ajustar o seu modelo devem ser diferentes daqueles que você utiliza para avaliar a performance de suas previsões. Os passos usuais a sere seguidos são:

(1) Divida a série temporal em duas partes: conjunto de treinamento e conjunto de teste (ou validação);

(2) Ajuste um modelo aos seus dados de treinamento;

(3) Faça previsões com base no modelo ajustado para o período de teste/validação;

(4) Verifique a performance através de gráficos e de medidas (além das já citadas, inclua o erro quadrático médio de previsão (PMSE)).

PMSE mede a discrepância entre as previsões do modelo e os dados observados. Suponha que nós temos uma série temporal de tamanho N e que reservemos M destas observações para validação, digamos

Após o ajuste do modelo com base na amostra de treinamento (com as primeiras N-M observações), nós estimamos (prevemos)

. O PMSE do modelo é então estimado através de,

Algumas vezes, é interessante calcular PMSE para várias opções de M a fim de verificar a robustez dos nossos resultados (a escolha de M é ad hoc, não costuma envolver um procedimento objetivo/automático).

As previsões podem ser (a) estáticas ou (b) dinâmicas. A previsão um passo adiante é igual para ambos os tipos. A diferença surge quando queremos prever mais adiante. A previsão dinâmica irá considerar, dentro do seu conjunto de informação, as previsões passadas (ela leva em conta o padrão passado da série) enquanto a previsão estática só considera os valores existentes da série no seu

conjunto de informação.

Em comparação com regressão, o número de hipóteses em um modelo (S)ARIMA é menor (precisamos garantir a estacionariedade (fraca) para conseguirmos fazer previsões razoáveis).

Vamos agora explorar um pouco mais alguns dos (11) pontos elencados acima.

Se, por exemplo, o processo cresce ao longo do tempo, será necessário tirar uma diferença (ou uma diferença fracionária se utilizarmos processos ARFIMA – F significa fracionário).

Se a variabilidade, por exemplo, cresce ao longo do tempo, será necessário transformar os dados (com alguma função côncava como logaritmo) para estabilizar a variância. Podemos utilizar para tal uma transformação da classe Box-Cox, isto é,

(estimamos o valor de lambda). Ás vezes, uma aplicação particular pode sugerir uma transformação apropriada. Por exemplo, suponha que um processo, digamos um investimento, pode ser descrito por:

, onde é o valor do investimento no tempo t e é uma mudança percentual do período t-1 para t (pode ser negativo). Aplicando logs, nós temos

ou onde é o operador primeira diferença, isto é, .

Se a percentagem é pequena em magnitude, então podemos

mostrar que e, portanto, que é um processo relativamente estável ( é chamado de retorno).

Após transformar adequadamente a série, o próximo passo envolve identificar preliminarmente as ordens autoregressivas, p, a ordem de diferenciação, d, e a ordem média móvel, q.

O gráfico da série temporal tipicamente sugere se devemos ou não diferenciar os dados. Crescimentos lineares (tendências) indicam que devemos tirar a primeira diferença e tendências quadráticas duas diferenças (e assim por diante).

Após tirar a diferença inspecionamos o gráfico da série transformada .

ARMAtoMA # prints psi-weights (para calcular as funções de resposta ao impulso no R).

O próximo passo são os diagnósticos (após a estimação).

Estatística Q de Ljung-Box:

Tipicamente H=20, isto é, para que consideremos que o modelo esteja bem ajustado aos dados (nossa informação disponível) consideramos, tipicamente, aqueles cujas primeiras vinte autocorrelações tenham um p-valor associado grande (tipicamente maior que 0,05 ou 0,10).

Diagnósticos no R: tsdiag(serie ajustada, lags) # R Checar padrões nos gráficos. Usualmente, podemos considerar como

outliers aqueles distanciados por mais de 3 desvios-padrão em magnitude.

Testes de Normalidade (Shapiro-Wilk, por exemplo): usualmente as distribuições dos processos econômicos apresentam caudas mais pesadas do que a distribuição normal prevê.

histqqnorm (serie ajustada$resid)shapiro.test(serie ajustada$resid)

Obs: O enfoque de Box & Jenkins para modelar uma série temporal é geralmente implicado pela suposição de que a dependência entre valores adjacentes no tempo é melhor explicada em termos de uma regressão dos valores correntes sobre os valores passados. Esta suposição é parcialmente justificada, na teoria, pela decomposição de Wolf.

A decomposição de Wolf nos diz que cada processo estacionário pode ser representado como a soma de um e uma sequência, digamos determinística.

A representação de Wolf não nos diz que esta decomposição é a melhor possível para descrever o processo, mas nos dá uma confiança de que não estamos utilizando um modelo “completamente fora do mercado” quando ajustamos um ARMA a uma série temporal.

Tabela: Comportamento da FAC e FACP nos modelos AR, MA e ARMA

SARIMA: o padrão de identificação da sazonalidade é também feita utilizando a FAC e a FACP. Neste caso, nós precisamos analisar os lags sazonais. Por exemplo, imagine que um processo coletado mensalmente tenha uma sazonalidade de período s=12 (meses). Neste caso, esperamos que existam autocorrelações na série de ordem mais elevada em 12, 24 e assim por diante.

Pense em 12 como o primeiro elemento da série (como o primeiro lag para o caso não sazonal), 24 como o segundo e assim por diante. Para que o processo seja estacionário ele não pode depender do tempo t, isto é, para que ele seja estacionário as autocorrelações devem decair (exponencialmente) para zero se olharmos lags de qualquer ordem (12, por exemplo).

Vejamos um exemplo. Os gráficos abaixo representam a FAC e a FACP de um processo estocástico (uma série temporal).

A série tem claramente uma tendência (dado o decaimento lento na FAC e o fato de que a FACP apresenta um valor perto de 1 – ambas as características indicam comportamento não estacionário). Seguindo o procedimento recomendado, a primeira diferença foi aplicada e os correlogramas resultantes da primeira

diferença, , são:

As figuras acima representam a FAC e a FACP da série transformada, isto é, (1-L)xt, onde L é o operador primeira diferença.

Note agora os picos nos lags sazonais, h= 1s, 2s, 3s, 4s, onde s=12 (logo, h= 12, 24, 36, 48) . O decaimento na FAC continua lento nos lags sazonais, o que sugere que tiremos uma diferença sazonal da série diferenciada, isto é,

onde e aqui.

Os correlogramas associados as novas transformações seguem:

Agora a série parece estacionária, pois não há nenhum decaimento lento para nenhum lag (nenhum múltiplo) .

Ambas as séries parecem ter um pico forte em h=1s (isto é, em h= 12). Na FAC vemos picos menores aparecendo em 2s e 3s. Na FACP temos um pico forte em 2s e picos menores em 3s e 4s.

Antes de explicar como encontrar o modelo “máximo”, vamos definir o modelo multiplicativo sazonal autoregressivo integrado média móvel (SARIMA).

Definição: O modelo multiplicativo sazonal autoregressivo integrado média móvel

(SARIMA) é dado por:

onde wt é o usual processo ruído branco. Notamos o modelo

por: ARIMA(p,d,q)x(P,D,Q)s. Os polinômios e representam as ordens p e q do processo ARMA comum,

respectivamente, enquanto os polinômios e representam as ordens P e Q da parte sazonal;

e .

Um modelo “máximo” poderia ser um ARIMA(2,1,2)x(4,1,2)12. Por quê? Colocamos novamente o gráfico abaixo para analisarmos melhor. Note que diferenciamos as séries (primeira diferença e diferença sazonal, logo d=1 e D=1). Existem dois lags não sazonais significativos (dentre os onze primeiros) na FACP, o que nos leva a considerar p=2. Na FAC, temos quatro lags não sazonais significativos, porém após decair até o terceiro o quarto volta a crescer em magnitude. Poderíamos considerar q= 4, mas na prática (do autor do texto), q=2 já seria uma boa escolha.

Para encontrar P e Q devemos olhar os lags sazonais. Na FACP temos 4 lags significativos (o terceiro e o quarto menores em magnitude, mas continua caindo). Por isso, considero P= 4 como uma alternativa para a ordem máxima. Na FAC temos dois lags significativos, logo, podemos Q =2.

Após a seleção do modelo máximo, temos que calcular os critérios de informação para ordens iguais ou menores (por exemplo, Q= 2 dificilmente será mantido com base na experiência do autor quando usarmos critérios de informação para selecionar o modelo, já que o valor ultrapassa em uma magnitude bem menor que o primeiro lag).

P.S. O melhor modelo para os dados deste exemplo, considerando AIC, foi ARIMA(1,1,1)x(2,1,1)12.

No R, podemos usar os seguintes comandos:

Segue uma tabela para explicação de como selecionar P e Q em um modelo SARIMA. Tails off = decai exponencialmente , cuts off = trunca.

Obs: Podemos calcular a função de resposta ao impulso (previsão do que acontece quando temos choques exógenos) representando as séries como um MA (infinito) no caso do processo conter um termo AR (se a série for estacionária). No R podemos usar a função ARMAtoMA para encontrar os valores da função de resposta ao impulso.

É possível fazer uma regressão com dados de séries temporais considerando que os resíduos são autocorrelacionados e seguem, por exemplo, um processo ARMA (ou qualquer outro dos citados neste texto – costuma-se chamar de Regressão Dinâmica neste contexto).

Neste caso, considere para explicação o seguinte modelo de regressão:

Os passos a serem seguidos são:

(1) Rode uma regressão (MQO) de yt sobre zt1,...,ztr (como se os erros não

fossem correlacionados). Guarde os resíduos ;

(2) Identifique os possíveis modelos ARMA para os resíduos (siga os passos já descritos);

(3) Rode MQG (mínimos quadrados generalizados, ou ainda MV, máxima verossimilhança) com os erros seguindo o processo que você encontrou (após todo o processo de seleção) em (2).

(4) Inspecione os resíduos, digamos, , do modelo encontrado em (4) para checar se ele parece ser um ruído branco (faça o teste de não autocorrelação usando a estatística Q de Ljung-Box, por exemplo). Ajuste o modelo se necessário.

Outras opções

- Redes Neurais

- Filtros: no entanto, eu descartaria esta parte, pois:

(i) Holt-Winters não sazonal: a função de previsão é um ARIMA(0,2,2);

(ii) EWMA (Exponentially Weighted Moving Average): a função de previsão é um ARIMA(0,1,1);

(iii) : as previsões de um modelo Holt-Winters (sazonal) são similares.

- Filtro de Hodrick-Prescott: a série temporal Yt é decomposta em tendência (componente de crescimento), digamos gt e um componente cíclico, isto é,

Yt=gt+ct,t=1,...,T. O filtro minimiza . O primeiro terno é a soma dos quadrados do componente cíclico ct =Yt - gt . O segundo termo é um múltiplo da soma dos quadrados das segundas diferenças dos componentes de tendência: quanto maior o valor de , maior a penalização e portanto mais suave será a série de tendência.

Hodrick & Prescott(1997) sugerem =1600 como razoável para dados trimestrais.

- Filtro de Butterworth