Estimação de Modelos ARMA

Este material tem como objetivo contribuir para o entendimento sobre o processo de estimação de modelos AR(p), MA(q) e ARMA(p,q). Para tanto, vamos mostrar como escrever uma função de verossimilhança geral para uma série temporal (aqui, os retornos) e alterá-la para cada tipo de modelo (AR, MA e ARMA).

INTRODUÇÃO

Assim como proposto por Box, Jenkins, and Reinsel (1994), o processo de modelagem de séries temporais univariadas por meio de modelos AR(p), MA(q) e ARMA(p,q) compreende os seguintes passos:

1. Identificar as ordens \(p\) e \(q\) do modelo
2. Estimar parâmetros
3. Verificar se os resíduos estimados se comportam como um ruído branco. Se sim, passa-se ao próximo passo. Caso contrário, retorna-se ao primeiro passo.
4. Fazer previsões

Como vimos no estudo dos modelos autorregressivos de médias móveis (ARMA), nem sempre é fácil identificar as ordens do modelo apenas usando as funções de autocorrelação (FAC) e autocorrelação parcial (FACP). Em função disso, precisamos também fazer uso de algum critério de informação (por exemplo, AIC ou BIC) que nos ajude com a correta decisão das ordens \(p\) e \(q\) (como mostramos para o o modelo ARMA neste link).

A escolha incorreta das defasagens pode ocasionar em resíduos indesejados (que não se comportam como um ruído branco). Lembre-se que no momento de derivação das propriedades dos modelos AR(p), MA(q) e ARMA(p,q) nós sempre tratamos \(a_t\) como um ruído branco com média \(0\), variância \(\sigma_{a}^{2}\), \(E\left[a_t-E(a)\right]\left[a_{t-l}-E(a)\right]=E[a_{t}a_{t-l}]=0\) e independente e identicamente distribuído (iid).

Neste documento, vamos mostrar como escrever a função de verossimilhança bem como seu logaritmo para cada modelo e entender quais testes devem ser aplicados para verificar os resultados obtidos.

IDENTIFICAR AS ORDENS \(p\) e \(q\) DO MODELO

Já abordamos anteriormente como fazer uso das funções de autocorrelação e autocorrelação parcial bem como dos critérios de informação. Você pode verificar nestes links (AR, MA e ARMA).

ESTIMAÇÃO

Nesta seção vamos mostrar como estimar os modelos AR, MA e ARMA usando-se o método de máxima verossimilhança. Optamos pelo método de máxima verossimilhança condicional que é assintoticamente equivalente ao da máxima verossimilhança exata. A vantagem do méotodo de máxima verossimilhança condicional é a facilidade de especificar e, consequentemente, estimar a função de verossimilhança. Porém, em caso de pequenas amostras o método não é tão eficiente quanto o método exato.

Vamos adotar a estratégia de derivar uma função de verossimilhança geral para modelos univariados estacionários. Após isso, podemos alterar tal função para considerar cada modelo com suas características.

• VEROSSIMILHANÇA GERAL

Em séries temporais as observações são dependentes. Logo, teremos uma função de probabilidade conjunta que dependerá da função de probabilidade condicional e da função de probabilidade marginal.

Assim, para um ativo financeiro, assuma que você tem dois retornos consecutivos \(r_{1}\) e \(r_{2}\). Neste caso temos que:

\[ f(r_{2},r_{1}) = f(r_{2}|r_{1})f(r_{1}) \] Para três retornos, \(r_{1}\), \(r_{2}\) e \(r_{3}\), a expressão se torna:

\[ f(r_{3},r_{2},r_{1}) = f(r_{3}|r_{2},r_{1})f(r_{2},r_{1}) \]

\[ f(r_{3},r_{2},r_{1}) = f(r_{3}|r_{2},r_{1})f(r_{2}|r_{1})f(r_{1}) \] Generalizando para uma série temporal de \(T\) retornos, teremos;

\[ f(r_{t},r_{t-1},...,r_{2},r_{1}) = f(r_{t}|r_{t-1},...r_{1}) f(r_{t-1},...,r_{1}) \] \[ f(r_{t},r_{t-1},...,r_{2},r_{1}) = f(r_{t}|r_{t-1},...r_{1})f(r_{t-1}|r_{t-2},...,r_{1})f(r_{t-2},...,r_{1}) \] \[ f(r_{t},r_{t-1},...,r_{2},r_{1}) = f(r_{t}|r_{t-1},...r_{1})f(r_{t-1}|r_{t-2},...,r_{1})f(r_{t-2}|r_{t-3},...,r_{1}) f(r_{t-3},...,r_{1}) \] \[ f(r_{t},r_{t-1},...,r_{2},r_{1}) = \left[\prod_{t=2}^{T}{f(r_{t}|r_{t-1},...r_{1})}\right]f(r_{1}) \]

Se \(r_{t}|r_{t-1},...r_{1}\) segue uma distribuição normal com média \(\mu\) e variância \(\sigma^{2}\) a função de verossimilhança condicional se torna, para \(R = \left\{r_{t}\right\}_{t=1}^{T}\):

\[ L(R|\boldsymbol{\theta}) = f(r_{t},r_{t-1},...,r_{1} | \boldsymbol{\theta}) = f(r_{1}|\boldsymbol{\theta}) \prod_{t=2}^{T}{ \frac{1}{\sqrt{2\pi\sigma^{2}}}exp\left\{\frac{-(r_{t}-\mu)^2}{2\sigma^{2}} \right\}} \]

Se ignoramos \(f(r_{1})\), a função de verossimilhança condicional dos retornos sob hipótese de normalidade será:

\[ L(R|\boldsymbol{\theta}) = f(r_{t},r_{t-1},...,r_{1}|\boldsymbol{\theta}) = \prod_{t=2}^{T}{ \frac{1}{\sqrt{2\pi\sigma^{2}}}exp\left\{\frac{-(r_{t}-\mu)^2}{2\sigma^{2}} \right\}} \]

Outras funções de distribuição de probabilidades (t-Student, log-Normal, Cauchy) podem ser usadas para lidar com os fatos estilizados discutidos neste material. A função de verossimilhança mudará para cada função de distribuição de probabilidades.

• VEROSSIMILHANÇA AR(p)

Suponha que temos um modelo AR(p) da seguinte forma:

\[ r_{t} = {\phi}_{0} + {\phi}_{1}r_{t-1} + {\phi}_{2}r_{t-2} + ... + {\phi}_{p}r_{t-p} + a_{t} \] Logo, \(E[r_{t}|r_{t-1},...,r_{t-p}] = \phi_{0}+\phi_{1}r_{t-1}+...+\phi_{p}r_{t-p}\) e \(Var(r_{t}|r_{t-1},...,r_{t-p}) = \sigma_{a}^{2}\) que significa que dado os retornos passados de \(r_{t}\) o retorno corrente é centralizado em torno de \(\phi_{0}+\phi_{1}r_{t-1}+...+\phi_{p}r_{t-p}\) com desvio padrão igual a \(\sigma_{a}\).

Assim, para um modelo AR(p) temos que a função de verossimilhança condicional será:

\[ L(R|\boldsymbol{\theta}) = f(r_{t},r_{t-1},...,r_{1};\boldsymbol{\theta}) = \prod_{t=p+1}^{T}{ \frac{1}{\sqrt{2\pi\sigma_{a}^{2}}}exp\left\{\frac{-(r_{t}-\phi_{0}-\phi_{1}r_{t-1}-...-\phi_{p}r_{t-p})^2}{2\sigma_{a}^{2}} \right\}} \] que no momento que observamos os valores de \(R = \left\{r_{t}\right\}_{t=1}^{T}\) temos como objetivo encontrar quais os parâmetros de \(\theta\) (neste caso, \(\phi_{0}\), \(\phi_{1}\), …, \(\phi_{p}\) e \(\sigma_{a}^{2}\)) que maximizam a função de verossimilhança ou seu logarítmo.

\[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=p+1}^{T}{\ln{f(r_{t},r_{t-1},...,r_{1};\boldsymbol{\theta})}} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=p+1}^{T}{\left[ \ln{1} - \ln{(2\pi\sigma_{a}^{2})^{\frac{1}{2}}} - \frac{1}{2\sigma_{a}^{2}}(r_{t}-\phi_{0}-\phi_{1}r_{t-1}-...-\phi_{p}r_{t-p})^{2} \right]} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=p+1}^{T}{\left[-\frac{1}{2}\ln{(2\pi\sigma_{a}^{2})} - \frac{1}{2\sigma_{a}^{2}}(r_{t}-\phi_{0}-\phi_{1}r_{t-1}-...-\phi_{p}r_{t-p})^{2} \right]} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = -\frac{(T-p)}{2}\ln{2\pi} -\frac{(T-p)}{2}\ln{\sigma_{a}^{2}} -\frac{1}{2\sigma_{a}^{2}}\sum_{t=p+1}^{T}{(r_{t}-\phi_{0}-\phi_{1}r_{t-1}-...-\phi_{p}r_{t-p})^{2}} \] onde \(\hat{\boldsymbol{\theta}} = [\phi_{0},\phi_{1},...,\phi_{p},\sigma_{a}^{2}]\) contem os parâmetros do modelo AR(p).

• VEROSSIMILHANÇA MA(1)

Suponha que temos o seguinte modelo de médias móveis de primeira ordem:

\[ r_{t} = \mu + a_{t} + {\theta}_{1}a_{t-1} \] onde \(a_{t}\) é um ruído branco com média \(0\), variância \(\sigma_{2}^{a}\), \(E[a_{t}a_{t-l}] = 0\) e independente e identicamente distribuído (iid). Supondo \(a_{0}=0\), para \(t=1\), temos:

\[ r_{1}=\mu + a_{1} + \theta_{1}a_{0} = \mu + a_{1} \] Isso implica que:

\[ (R_{1}=r_1|a_0=0) \sim N\left(\mu,{\sigma}_{a}^{2}\right) \]

em função das propriedades estatísticas de \(a_t\). Uma vez que \(r_1\) é observado, temos que \(a_{1} = r_{1} - \mu\) o que implica que:

\[ (R_{2}=r_{2}|R_{1}=r_{1},a_{0}=0) \sim N\left(\mu + \theta_{1}a_{1},{\sigma}_{a}^{2}\right) \] e a partir que \(a_{1}\) é conhecido, podemos obter \(a_{2}\) por meio de \(a_{2}=r_{2}-\mu-\theta_{1}a_{1}\). Assim, podemos observar que dado que conhecemos \(a_{0}=0\), a sequência \(\left\{a_{t}\right\}_{t=1}^{T}\) pode ser calculada a partir de \(\left\{r_{t}\right\}_{t=1}^{T}\) iterando sobre:

\[ a_{t} = r_{t} - \mu -\theta_{1}a_{t-1} \] para \(t = 1, 2, ..., T\), iniciando a partir de \(a_{0}=0\). A verossimilhança condicional da amostra completa pode ser calculada pelo produto das densidades individuais:

\[ L(R|\boldsymbol{\theta}) = f(r_{t},r_{t-1},...,r_{1},a_{0}=0|\boldsymbol{\theta}) = \prod_{t=2}^{T}{ \frac{1}{\sqrt{2\pi\sigma_{a}^{2}}}exp\left\{\frac{-a_{t}^2}{2\sigma_{a}^{2}} \right\}} \] que no momento que observamos os valores de \(R = \left\{r_{t}\right\}_{t=1}^{T}\) temos como objetivo encontrar quais os parâmetros de \(\theta\) (neste caso, \(\mu\), \(\theta_{1}\) e \(\sigma_{a}^{2}\)) que maximizam a função de verossimilhança ou seu logarítmo.

\[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=1}^{T}{\ln{f(r_{t},r_{t-1},...,r_{1},a_{0}=0;\boldsymbol{\theta})}} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=1}^{T}{\left[ \ln{1} - \ln{(2\pi\sigma_{a}^{2})^{\frac{1}{2}}} - \frac{a_{t}^{2}}{2\sigma_{a}^{2}} \right]} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=1}^{T}{\left[-\frac{1}{2}\ln{(2\pi\sigma_{a}^{2})} - \frac{a_{t}^{2}}{2\sigma_{a}^{2}} \right]} \] \[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = -\frac{T}{2}\ln{2\pi} -\frac{T}{2}\ln{\sigma_{a}^{2}} - \sum_{t=1}^{T}{-\frac{a_{t}^{2}}{2\sigma_{a}^{2}}} \] Esta é a verossimilhança condicional usada para estimar tanto modelos MA(1) e MA(q). Sendo que a única diferença entre eles é que enquanto assumimos no MA(1) que \(a_{0}=0\), no MA(q) é preciso assumir que os valores iniciais para todos os termos de erro são nulos, ou seja, \(a_{0}=a_{-1}=...=a_{q-1}=0\)

• VEROSSIMILHANÇA ARMA(p,q)

Assuma que observamos os valores iniciais para \(R = \left\{r_{t}\right\}_{t=1}^{T}\). Neste caso, inicia-se a estimação em \(t=p+1\) e fixa-se:

\[ a_{p} = a_{p-1} = ... = a_{p-1+1} = 0 \]

e a função de verossimilhança se torna:

\[ L(R|\boldsymbol{\theta}) = f(r_{t},r_{t-1},...,r_{1}|r_p,r_{p-1},...,r_1,a_{p} = a_{p-1} = ... = a_{p-1+1} = 0|\boldsymbol{\theta}) = \prod_{t=2}^{T}{ \frac{1}{\sqrt{2\pi\sigma_{a}^{2}}}exp\left\{\frac{-a_{t}^2}{2\sigma_{a}^{2}} \right\}} \]

que no momento que observamos os valores de \(R = \left\{r_{t}\right\}_{t=1}^{T}\) temos como objetivo encontrar quais os parâmetros de \(\theta\) (neste caso, \(\mu\), \(\phi_1\), …, \(\phi_p\), \(\theta_{1}\), …., \(\theta_{p}\) e \(\sigma_{a}^{2}\)) que maximizam a função de verossimilhança ou seu logarítmo.

\[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = \sum_{t=1}^{T}{\left[-\frac{1}{2}\ln{(2\pi\sigma_{a}^{2})} - \frac{a_{t}^{2}}{2\sigma_{a}^{2}} \right]} \]

\[ \ln{L(\hat{\boldsymbol{\theta}}|R)} = -\frac{T-p}{2}\ln{2\pi} -\frac{T-p}{2}\ln{\sigma_{a}^{2}} - \sum_{t=p+1}^{T}{-\frac{a_{t}^{2}}{2\sigma_{a}^{2}}} \]

É importante salientar que todas as funções de verossimilhanças mostradas assumem a hipótese de que a série temporal dos retornos segue uma distribuição Normal. Caso contrário, por exemplo uma t-Student, as funções devem ser escritas assumindo uma nova hipótese de distribuição.

DIAGNÓSTICO DOS RESÍDUOS

Uma hipótese nos modelos ARMA estacionários é que o termo de erro \(a_{t}\) é um ruído branco. Isso implica que \(a_{t}\) tem média \(0\), variância \(\sigma_{a}^{2}\), \(E[a_{t}a_{t-l}] = 0\) e indepente e identicamente distribuído (iid).

Desta forma, espera-se que os resíduos estimados pelo modelo se comportem como um ruído branco se o modelo está corretamente especificado.

Outra hipótese assumida na derivação da função de verossimilhança dos modelos AR(p), MA(q) e ARMA(p,q) é que a série temporal dos retornos segue uma distribuição Normal o que, por conseguinte, define que os resíduos também têm a mesma distribuição de probabilidade. Assim, devemos fazer um teste estatístico que confirme se os resíduos são normalmente distribuídos.

Resumidamente, temos o seguinte conjunto de testes sobre os resíduos:

• \(E[a] = 0\): o próprio processo de estimação dos parâmetros garante esta hipótese
• \(\sigma_{a}^{2}\): podemos usar a FAC dos resíduos ao quadrado para verificar se os resíduos são homocedásticos ou se apresentem heterocedasticidade condicional. O ideal é que nenhuma defasagem seja estatisticamente significante. Caso contrário, talvez seja preciso estimar a variância condicional para a série tem análise (aprenderemos isso com os modelos ARCH e GARCH)
• \(E[a_{t}a_{t-l}] = 0\): a FAC dos resíduos permite avaliar se os mesmos são autocorrelacionados. O ideal é que nenhuma defasagem seja estatisticamente significante.
• Resíduos seguem uma distribuição Normal: podemos fazer uso do teste Jarque-Bera que verifica se os momentos da série estimada são iguais aos da Normal. Esperamos não rejeitar a hipótese nula de normalidade dos resíduos.

EXEMPLO:

Vamos simular um modelo ARMA e fazer uso da abordagem proposta por Box, Jenkins, and Reinsel (1994) da seguinte forma:

1. Identificar diversas ordens \(p\) e \(q\) por meio funções de autocorrelação (FAC) e autocorrelação parcial (FACP) para uma série temporal de interesse
2. Estimar os modelos por meio do princípio da máxima verossimilhança usando pacotes do R;
3. Escolher dentre os modelos estimados o que apresenta menor AIC ou BIC;
4. Verificar se os resíduos do modelo se comportam como um ruído branco.

ARMA(1,1)

Suponha que temos o modelo ARMA(1,1) no seguinte formato:

\[ r_{t} = 10 + 0.4r_{t-1} + a_{t} + 0.3a_{t-1} \] onde \(a_{t}\) é um ruído branco com média \(0\), variância \(\sigma_{a}^{2}\), \(E\left[ \left(a_t-\bar{a}\right)\left(a_{t-l}-\bar{a}\right)\right]=E[a_{t}a_{t-l}]=0\) e independente e identicamente distribuído (iid).

O gráfico abaixo apresenta a série temporal que representa tal modelo.

É possível observar que realmente há estacionariedade fraca dado que a série oscila em torno de uma média (\(10\)) e com uma variância constante.

Agora, imagine que temos a série temporal acima e que não sabemos que ela foi originada por um processo ARMA(1,1), mas gostaríamos de usar a abordagem proposta por Box, Jenkins, and Reinsel (1994) para estimar um modelo para a série. Para tanto, vamos seguir os passos propostos anteriormente.

• IDENTIFICAÇÃO DAS ORDENS \(p\) e \(q\):

É preciso avaliar as funções de autocorrelação e autocorrelação parcial para verificar as ordens que devemos especificar para o modelo. Lembre-se que vimos anteriormente que a FACP é ideal para definir a ordem \(p\) do AR enquanto a FAC deve ser usada para encontrar a ordem \(q\) do MA. Abaixo, os gráficos das duas funções.

Aparentemente parece que se trata de um modelo ARMA(2,2), pois para a FAC e FACP há significância até a defasagem 2 [^1]. Observe que, como esperado, o uso da FAC e FACP nos direcionaria a uma especificação que não é verdadeira, pois sabemos que a série temporal foi gerada por um modelo ARMA(1,1). Assim, como proposto no processo de estimação, vamos fazer uso de algum critério de informação que nos ajude na seleção da melhor especificação para a série que queremos estimar.

Uma abordagem útil é fazer todas as combinações possíveis de modelos ARMA para \(p=0\) até \(p=2\) e \(q=0\) até \(q=2\), como segue:

• ARMA(0,0)
• ARMA(1,0)
• ARMA(2,0)
• ARMA(0,1)
• ARMA(1,1)
• ARMA(2,1)
• ARMA(0,2)
• ARMA(1,2)
• ARMA(2,2)

e a partir disso estimar todos os modelos e escolher aquele com menor AIC e/ou BIC.

• ESTIMAÇÃO:

Uma vez que temos os modelos de interesse para avaliação do AIC ou BIC, vamos fazer a estimação usando o princípio da máxima verossimilhança. Desta forma, uma alternativa é escrever a função de verossimilhança no R e por processo de otimização buscar os parâmetros que maximizam o seu logaritmo. Como optamos na simulação dos dados que \(a_{t}\) seguisse uma distribuição Normal, em caso de criar uma função de verossimilhança no R, tal distribuição deveria ser usada.

Porém, aqui usamos a função arima do pacote stats para estimar os modelos selecionados que usa como default a distribuição Normal para o termo de erro. Caso você tenha interesse em entender outras opções da função, faça help(arima) no console do R para visualizar a documentação da função.

Como era de se esperar, o modelo ARMA(1,1) apresenta menor AIC e BIC. Apesar deste modelo não ter o maior valor do logaritmo da verossimilhança, ele é o melhor modelo dentre os demais. Lembre-se que no cálculo do AIC ou BIC há penalização para modelos com maior quantidade de parâmetros.

Uma vez selecionado o modelo, vamos verificar quais os parâmetros estimados e a significância estatística de cada. Observe que os parâmetros estimados estão bem próximos dos verdadeiros parâmetros populacionais que usamos para simular a série temporal.

Resultado Estimação modelo ARMA(1,1)
=============================================
                      Dependent variable:    
                  ---------------------------
                            arma11           
---------------------------------------------
ar1                        0.349***          
                            (0.048)          
                                             
ma1                        0.318***          
                            (0.047)          
                                             
intercept                  10.029***         
                            (0.064)          
                                             
---------------------------------------------
Observations                 1,000           
Log Likelihood            -1,414.004         
sigma2                       0.990           
Akaike Inf. Crit.          2,836.007         
=============================================
Note:             *p<0.1; **p<0.05; ***p<0.01

• DIAGNÓSTICO DO MODELO SELECIONADO:

Uma hipótese nos modelos ARMA estacionários é que o termo de erro \(a_{t}\) é um ruído branco. Isso implica que \(a_{t}\) tem média \(0\), variância \(\sigma_{a}^{2}\), \(E[a_{t}a_{t-l}] = 0\) e indepente e identicamente distribuído (iid) (aqui, supomos que era Normalmente distribuído). Desta forma, espera-se que os resíduos estimados pelo modelo se comportem como um ruído branco se o modelo está corretamente especificado.

Para a hipótese de não autocorrelação dos resíduos (\(E[a_{t}a_{t-l}] = 0\)), podemos usar a função de autocorrelação (FAC) dos resíduos, onde esperamos que nenhuma defasagem tenha significância estatística. Como mostra o gráfico abaixo, não existe autocorrelação nos resíduos do modelo escolhido, pois as defasagens \(l>0\) não são estatisticamente significantes.

Já a hipótese de homocedasticidade dos resíduos (\(\sigma_{a}^{2}\)), podemos usar o gráfico da função de autocorrelação serial (FAC) dos resíduos ao quadrado. Como mostra o gráfico abaixo, não há heterocedasticidade condicional, ou seja, os resíduos se comportam conforme a hipótese de homocedasticidade.

Outra hipótese assumida na derivação da função de verossimilhança dos modelos AR(p), MA(q) e ARMA(p,q) é que a série temporal dos retornos segue uma distribuição Normal o que, por conseguinte, define que os resíduos também têm a mesma distribuição de probabilidade. Uma alternativa para isso é o teste Jarque Bera.

    Jarque Bera Test

data:  modelo[[5]]$residuals
X-squared = 0.88234, df = 2, p-value = 0.6433

A hipótese nula do teste é que os dados são normalmente distribuídos contra a hipótese alternativa que os dados não são normalmente distribuídos. Sabemos que o p-valor é o nível de significância mínimo para se rejeitar a hipótese nula. Assim, queremos um p-valor alto (\(>0,05\)) o que não nos permite rejeitá-la. Observe que é exatamente isso que o resultado do nosso teste mostra e, assim, não podemos rejeitar a hipótese nula de normalidade nos dados.

REFERÊNCIAS

Box, George EP, Gwilym M Jenkins, and Gregory C Reinsel. 1994. Time Series Analysis Forecasting and Control. Englewood Cliffs Prentice Hall.

Campbell, John Y, Andrew Wen-Chuan Lo, and Archie Craig MacKinlay. 1997. The Econometrics of Financial Markets. Princeton (NJ) Princeton University Press.

Morettin, Pedro Alberto. 2008. Econometria Financeira Um Curso Em Séries Temporais Financeiras. Edgard Blucher.

Tsay, Ruey S. 2010. Analysis of Financial Time Series. John Wiley & Sons.

———. 2014. An Introduction to Analysis of Financial Data with R. John Wiley & Sons.