Embed Size (px)
Citation preview
Universidade de Brasília - UnBInstituto de Exatas
Departamento de Estatística
Modelagem de Retornos Climatológicos
via Modelos GARCH
Gabriel Fonseca Sarmanho - 06/19370
Brasília
2009
GABRIEL FONSECA SARMANHO - 06/19370
Modelagem de Retornos Climatológicos via
modelos GARCH
Relatório apresentado à disciplina Estágio Supervisionado
II do curso de graduação em Estatística, Departamento de
Estatística, Instituto de Exatas, Universidade de Brasília,
como parte dos requisitos necessários para o grau de
Bacharel em Estatística.
Orientação: Prof. Afrânio Márcio Corrêa Vieira
Co-orientação: Prof. Paulo Sérgio Lucio - Depto de Estatística / UFRN
Brasília - DF
2009
Agradecimentos
À “Tia Cida”, minha querida mãe, por ser essa pessoa tão incrível e por sua compreensãoneste ano tão complicado; pela educação, exemplos, atitudes, perseverança, que foramdados; pelo amor incondicional e mútuo entre nós; por sua crença e confiança; enfim, porser a pessoa que mais amo na vida.
Aos meus irmãos, Izabela e Léo, por simplesmente existirem; pelo apoio em todas asetapas, admiração mútua e certeza da vossa presença.
À Tia Liliam, minha segunda mãezinha, por absolutamente tudo em minha vida; por terme acolhido em “nossa” casa, por ser engraçada, divertida, carinhosa, muito querida pormim.
Ao meu pai pela confiança e respeito aprendidos durante os últimos anos.
À toda família Fonseca e família Sarmanho: tios, primos e avós, em especial in memorianao meu avô materno Antônio, exemplo a ser seguido.
Ao meu orientador, Afrânio, pela confiança neste trabalho, enorme prestatividade em to-dos os momentos desde o estágio no INMET, pelos conselhos quanto à vida acadêmica eprincipalmente pelos conselhos de vida.
Ao Paulo, meu coorientador, que mesmo à distância, foi o idealizador deste e esteve sem-pre disposto a tirar dúvidas, acreditou em meu potencial e me fez acreditar no resultadodesta pesquisa
Aos professores do Departamento de Estatística - UnB, por todo o conhecimento difun-dido. Aos funcionários do departamento pela paciência durante esses quatro anos.
Aos membros da banca que contribuíram com sugestões para a melhoria deste trabalho.Agradeço também à professora Cibele por toda a atenção e disposição em me ajudarquando precisei.
Aos amigos do INMET: Dr. Lauro, Fábio, Andrea, Dani, André, Mozar, Dida, pelos doisanos de convivência, aprendizado e compreensão. Ao grande amigo Fabrício que con-tribui de forma substancial nessa monografia e em todos os trabalhos no INMET, sempre
com paciência, dedicação, humor e confiança.
Aos meus amigos de infância que nos últimos anos tenho pouco encontrado, mas queestão guardados em minha memória.
Aos colegas de ensino médio que, mesmo um pouco distantes, em todos os encontros sãosempre divertidos, companheiros, simplesmente amigos.
Aos “computadores amigos” do Tiago (“Presunto”), da Luna, do Vinícius, do Fabrício edo Wallace, que sem eles os resultados não sairiam tão “rapidamente”, ou melhor, nãoteriam saídos.
Aos amigos formados (“prováveis formandos”) durante a graduação: Gilson, Fred, Mair-con, Putão, Serenata, Perigo, Gago e tantos outros que ajudaram em momentos difíceis,sorriram em momentos felizes (e muitos foram); por terem feito esses quatro árduos anosbem mais prazerosos em minha vida.
E, final e novamente, à Luna, pelos sorrizos, pelos choros, pelas brincadeiras, pelas con-versas séries, por ser uma pessoa tão especial em minha vida, por trazer lembranças ale-gres nesses últimos anos, pela parceria no curso de graduação e por não ter somente feitoo curso de graduação comigo. Serei eternamente grato por todo o aprendizado.
4
Climate is what we expect, weather is what we get.
Mark Twain
5
Resumo
Neste trabalho serão apresentados os aspectos teóricos de uma classe de mod-elos heteroscedásticos condicionais, mais conhecidos como modelos da família ARCH(autoregressive conditional heteroscedasticity). A proposta é verificar se tais modelos,largamente utilizado em cenários financeiros, podem ser também utilizados no âmbitometeorológico, trazendo resultados satisfatórios, consistentes e representativos para pre-visões de volatilidade a curto prazo de variáveis que caracterizam o clima. Serão feitosdois enfoques de estimação: frequentista e bayesiano, apresentando suas aplicabilidadese verificando qual dentre as diferentes combinações de modelos melhor se ajusta à cadauma de seis séries de retornos de extremos climáticos fornecidas pelo Instituto Nacionalde Meteorologia (INMET). Salienta-se aqui o forte uso de recursos computacionais naconstrução dos algoritmos necessários para a modelagem estatística, todos programadosem R, devido à inexistência de ferramentas prontas em pacotes de tal linguagem.
Palavras-chave: Autoregressive Conditional Heteroskedasticity, Generalized Autoregres-sive Conditional Heteroskedasticity, Extremos Climáticos, Inferência Bayesiana, métodosMCMC.
Sumário
1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.2 Específicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
3 Aplicação na Área Climática . . . . . . . . . . . . . . . . . . . . . . . . . . 33.1 Conceitos de Climatologia . . . . . . . . . . . . . . . . . . . . . . . . . 33.2 Variáveis climatológicas . . . . . . . . . . . . . . . . . . . . . . . . . . 4
4 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Fundamentação Teórica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
5.1 Séries Temporais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85.1.1 Aspectos gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . 85.1.2 Séries Temporais Auto-Regressivas . . . . . . . . . . . . . . . . 115.1.3 Estrutura dos dados - Os Retornos . . . . . . . . . . . . . . . . . 12
5.2 Modelo ARCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135.2.1 Princípios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135.2.2 Estrutura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145.2.3 Propriedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
5.3 Modelo GARCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.3.1 Princípios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.3.2 Estrutura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.3.3 Propriedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175.3.4 Limitações do Modelo GARCH . . . . . . . . . . . . . . . . . . 185.3.5 Identificação de Modelos GARCH . . . . . . . . . . . . . . . . . 18
5.4 Inferência Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195.4.1 Conceitos importantes . . . . . . . . . . . . . . . . . . . . . . . 205.4.2 Inferência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.5 Simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245.5.1 Métodos de Integração de Monte Carlo . . . . . . . . . . . . . . 25
5.6 Cadeias de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 275.6.1 Conceitos Importantes . . . . . . . . . . . . . . . . . . . . . . . 275.6.2 Probabilidades Limites e Distribuições Estacionárias . . . . . . . 295.6.3 Espaço de Estados Contínuo . . . . . . . . . . . . . . . . . . . . 31
6 Estimação de Modelos GARCH . . . . . . . . . . . . . . . . . . . . . . . . . 326.1 Método Clássico: Otimização de Funções Não-Lineares . . . . . . . . . . 35
6.1.1 Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . 35
8
6.1.2 Método Scoring de Fisher . . . . . . . . . . . . . . . . . . . . . 366.1.3 Métodos Quase-Newton . . . . . . . . . . . . . . . . . . . . . . 36
6.2 Método Bayesiano - MCMC . . . . . . . . . . . . . . . . . . . . . . . . 376.2.1 Princípios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376.2.2 Amostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . 386.2.3 Amostrador de Griddy-Gibbs . . . . . . . . . . . . . . . . . . . 396.2.4 Algoritmo de Metropolis-Hastings . . . . . . . . . . . . . . . . . 40
7 Escolha e adequação do modelo . . . . . . . . . . . . . . . . . . . . . . . . 437.1 Critérios Utilizadas na escolha . . . . . . . . . . . . . . . . . . . . . . . 43
7.1.1 Akaike Information Criterion - AIC . . . . . . . . . . . . . . . . 437.1.2 Bayesian Information Criterion - BIC . . . . . . . . . . . . . . . 43
7.2 Testes de adequação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 447.2.1 Teste de Ljung-Box . . . . . . . . . . . . . . . . . . . . . . . . . 457.2.2 Teste de Jarque-Bera . . . . . . . . . . . . . . . . . . . . . . . . 45
7.3 Diagnósticos de Convergência . . . . . . . . . . . . . . . . . . . . . . . 457.4 Adequabilidade do modelo . . . . . . . . . . . . . . . . . . . . . . . . . 47
8 Previsão da Volatilidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . 489 Aplicação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
9.1 Um estudo de Simulação . . . . . . . . . . . . . . . . . . . . . . . . . . 509.1.1 Análise Descritiva . . . . . . . . . . . . . . . . . . . . . . . . . 509.1.2 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
9.2 Aplicação a dados reais . . . . . . . . . . . . . . . . . . . . . . . . . . . 599.2.1 Análise descritiva dos retornos . . . . . . . . . . . . . . . . . . . 609.2.2 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 659.2.3 Previsão de volatilidade . . . . . . . . . . . . . . . . . . . . . . 82
9.3 Aspectos Computacionais . . . . . . . . . . . . . . . . . . . . . . . . . . 8610 Conclusão e Trabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . . . 89Bibliografia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91APÊNDICE A - Demonstrações das propriedades dos modelos ARCH e GARCH 94APÊNDICE B - Métodos via simulação estocástica . . . . . . . . . . . . . . . . 97APÊNDICE C - Critérios de escolha AIC e BIC para as seis variáveis . . . . . . 99APÊNDICE D - Diagnóstico de Heidelberg para as seis variáveis . . . . . . . . 100APÊNDICE E - Diagnóstico de Geweke para as seis variáveis . . . . . . . . . . 101APÊNDICE F - Verossimilhança para a Distribuição T de Student . . . . . . . . 102APÊNDICE G - Cálculo do Vetor Score e da Matriz Hessiana . . . . . . . . . . 104
1 Introdução
Meteorologia é a ciência que estuda os fenômenos da atmosfera terrestre visando,a partir da dinâmica de processos físicos, a previsão de fenômenos em um período curto detempo (previsão de tempo), assim como explorar tendências de flutuações climáticas emperíodos longos pré-determinados como semanas, meses ou até anos (previsão climática).
Os Modelos Dinâmicos, que levam em consideração a complexidade de fenô-menos físicos e tentam simular a dinâmica atmosférica terrestre, mostram-se melhor ca-pacitados à realização de previsões de tempo (ver Moura, D. (1995)); enquanto os Mode-los Estocásticos, ou seja, aqueles que consideram no processo a existência de componentealeatória estocástica, apesar de menos complexos em comparação aos dinâmicos, são maissatisfatórios para a previsão climática, caracterizada por ser de longo prazo.
Vários modelos, como regressão linear simples e múltipla, regressão por compo-nentes principais (RCP), modelos de séries temporais (Box e Jenkins, Holt-Winters, etc),entre outros, têm sido experimentados com graus variados de sucesso.
Ainda que não identificado, na literatura das ciências meteorológicas, o conceitode volatilidade, isto é, de blocos de variação (variância não constantes e dependente tem-poralmente) em séries de quaisquer variáveis climatológicas, é presente. Assim sendo,deve-se recorrer aos modelos ditos heteroscedásticos condicionais, mais conhecidos comofamília ARCH (do inglês Autoregressive Conditional Heteroscedasticity), para tratar a járeferida anomalia.
Os modelos ARCH não foram ainda amplamente explorados, talvez devido à car-acterização tida das principais variáveis estudas na previsão climática (precipitação etemperatura) ou pelo desconhecimento das técnicas estatísticas mais sofisticadas nessecampo. Portanto, um estudo de séries climatológicas fazendo uso de tais modelos podeapontar alternativas poderosas aos modelos utilizados tradicionalmente.
2 Objetivos
2.1 Gerais
• Verificar se modelos da família ARCH, largamente utilizado para cenários finan-ceiros, pode ser também utilizado, trazendo resultados satisfatórios e consistentes,no âmbito meteorológico;
• Desenvolver e implementar uma metodologia para previsão de dados diárias noInstituto Nacional de Meteorologia (INMET);
2.2 Específicos
• Estudar aplicar a metodologia de séries temporais associada aos modelos da famíliaARCH à séries diárias de variáveis climatológicas;
• Estimar os parâmetros de modelos da família ARCH, de acordo com os métodos:clássico e bayesiano;
• Implementar métodos computacionais associados aos métodos de estimação dosparâmetros dos modelos e comparar o desempenho de ambos;
• Identificar, caracterizar e fazer previsões de volatilidade, verificando a qualidade doajuste
2
3 Aplicação na Área Climática
3.1 Conceitos de Climatologia
Meteorologia é uma ciência muito vasta e complexa, principalmente devido à suaextensão e a infinidade de tipos de ambientes que podem influenciar tanto na composição,como na dinâmica dos movimentos atmosféricos, gerando condições diversas de tempo eclima em todas as partes do globo terrestre.
Basicamente, algumas variáveis são suficientes para caracterizar o tempo e o climade uma região. Seriam elas: temperatura do ar, umidade relativa, pressão atmosférica,velocidade e direção dos ventos, precipitação em seus diversos tipos, radiação solar equantidade de horas com brilho solar (insolação).
Uma base de dados observados destas variáveis, entre outras, é de fundamentalimportância para aplicação não só na previsão do tempo, mas em outras áreas de im-portância econômica e social para um País, como agropecuária, defesa civil, indústria,turismo, transporte aéreo e marítimo, ecologia, recursos hídricos, energia e saúde pública.
Atualmente há muita curiosidade e bastante esforço em entender a questão defenômenos meteorológicos extremos. Eventos extremos deflagram problemas principal-mente em camadas sociais menos privilegiadas economicamente. Com a comprovação, apartir de dados observados, de que em muitas regiões já são detectadas modificações nassuas características climáticas, uma metodologia que dê suporte a previsões diárias dasvariáveis citadas é muito pertinente.
Com previsões precisas das variáveis supracitadas, a sociedade poderia, por ex-emplo, se prevenir com antecedência com relação a episódios de ondas de calor, a gerarpolíticas mais eficientes que auxiliasse a mitigação em áreas de desertificação, a um plane-jamento por parte do setor de saúde contra efeitos indesejáveis de doenças que se intensi-ficam sob certas condições de tempo.
3
4
3.2 Variáveis climatológicas
Para o estudo serão consideradas certas variáveis climatológicas que aparente-mente possuem uma característica desejável à modelagem: grande variabilidade e estaaparentemente não constante no tempo. Normalmente são variáveis conhecidas como“extremos climatológicos” que assim se comportam, tais como: temperatura máxima,temperatura mínima, umidade relativa máxima, umidade relativa mínima, média da ve-locidade do vento, insolação e radiação solar média.
Informações sobre tais eventos extremos para ainda são poucas na literatura es-pecializada, onde o enfoque principal é dado aos eventos extremos de precipitação. Estetrabalho visa desenvolver uma metodologia para previsão de dados diários e avaliar con-seqüentemente sua eficiência em simular bem casos extremos.
Os dados usados neste estudo serão a nível diário, de séries temporais do InstitutoNacional de Meteorologia, para o período de janeiro de 2007 a dezembro de 2008. Taisdados são medidos por dois tipos de estações meteorológicas:
• Estações automáticas: compostas de uma unidade de memória central (“data log-ger”), ligadas a vários sensores dos parâmetros meteorológicos, que integra os val-ores observados minuto a minuto e os atualiza automaticamente a cada hora. Afigura (3.1a) abaixo mostra uma foto de estação automática de Brasília.
• Estações Convencionais: compostas de vários sensores isolados que registramcontinuamente os parâmetros meteorológicos, que são lidos e anotados por um ob-servador a cada intervalo e este os envia a um centro coletor por um meio de comu-nicação qualquer. A figura (3.1b) abaixo representa uma estação convencional emBrasília.
Mais especificamente, a temperatura para fins meteorológicos é medida em graus,onde a escala Celsius (C) é a mais comum e adotada, cujos instrumentos são os ter-mômetros. Tais termômetros, para que indiquem a leitura representativa da temperaturado ar, são protegidos da radiação do sol, da terra e das radiações de quaisquer objetoscircunvizinhos, sendo acondicionados de maneira a serem protegidos destes fatores. Taisinstrumentos medem a temperatura instantânea do ar, temperaturas máximas e mínimas, atemperatura máxima geralmente é observada próxima as 15 horas e as mínimas próximasas 4 horas da manhã.
A medição da umidade relativa do ar segue padrão análogo a da temperatura, sendoo higrômetro o instrumento utilizado para medir tal variável. A unidade de medida é dadaem %. A umidade máxima registrada geralmente se dá nas primeiras horas do dia e amínima próxima às 15 horas da tarde.
5
(a) Estação Automática (b) Estação Convencional
Figura 3.1: Estações em Brasília (Fonte: INMET)
A insolação ou a duração do brilho solar é registrada pelo heliógrafo, em horas edécimos. A evaporação é medida pelo evaporímetro de Piche, em mililitro (ml) ou emmilímetros de água evaporada (mm), a partir de uma superfície porosa, mantida perma-nentemente umedecida por água. A evaporação também pode ser medida, em mm, deuma superfície livre de água, chamada tanque evaporimétrico classe A, como é o caso dasmedidas nas estações convencionais do INMET.
A velocidade instantânea do vento, em metros por segundo (m/s), pode ser obtidapor meio dos anemômetros ou dos anemógrafos.
4 Metodologia
A dependência temporal é uma característica notável das variáveis observadas nocampo da meteorologia e, nesse sentido, a análise de modelos baseados em séries tempo-rais é útil na identificação dos mecanismos (processos) que descrevem o comportamentode tais variáveis climatológicas.
Basicamente foi feito o uso da modelagem de séries temporais com o objetivo deestabelecer modelos que permitam obter previsões e elaborar cenários úteis na tomada dedecisões. Particularmente, os modelos da família ARCH foram explorados no trabalho,com o intuito de se prever a volatilidade tão presente nas séries diárias analisadas.
A estimação do modelo foi feita segundo o método da máxima verossimilhança,abordando-se a metodologia clássica, porém dado um enfoque na metodologia Bayesiana.Para o caso clássico, devido ao caráter não-linear do problema de estimação, este é re-solvido segundo métodos de estimação numérica (Newton-Raphson, Scoring de Fisher,Quase-Newton).
Quando utiliza-se a inferência Bayesiana, estimativas exatas podem ser feitas nosentido de não confrontar-se o problema da não-linearidade visto no caso clássico. Porém,a dificuldade da resolução da estimativas utilizando-se a distribuição à posteriori do mo-delo requer o uso de métodos de simulação estocástica via cadeias de Markov. Tais méto-dos permitem obter amostras das distribuições marginais a posteriori dos parâmetros.Baseado nestas amostras, serão feitas análises de diagnósticos para: avalição da con-vergência da cadeia, análise de adequação do modelo e interpretação das medidas clima-tológicas de interesse, no caso a volatilidade da série.
O corpo do texto desta monografia está organizado segundo 10 capítulos, alémde 3 apêndices. Para melhor entendimento da dinâmica do trabalho, abaixo há descriçãosucinta de cada um destes capítulos:
Capítulos 1 ao 4: Apresentação do problema proposto, os objetivos gerais e es-pecíficos do trabalho, além da metodologia a ser utilizada para a resolução do problema.Apresenta-se ainda uma breve descrição do tema usado na aplicação das técnicas;
6
7
No capítulo 5 encontra-se toda a fundamentação teórica do trabalho, isto é, a ap-resentação de todos os conceitos utilizados na literatura que são fundamentais para o en-tendimento das técnicas utilizadas. São conceitos de: séries temporais; propriedades dosModelos ARCH e GARCH; aspectos importantes de inferência bayesiana e os métodosde simulação de Monte Carlo e de Cadeias de Markov, para o entendimento dos métodosMCMC.
No capítulo 6 estão esmiuçados os passos para a estimação dos parâmetos dosmodelos GARCH segundo as metodologias frequentista e bayesiana, sendo detalhadosos métodos de otimização numérica para o primeiro caso e os métodos MCMC para osegundo.
No capítulo 7 encontram-se os argumentos utilizados para a seleção de modelos:critérios AIC e BIC. Além disso, testes de adequação do modelo, feitos posteriormenteà estimação dos parâmetros, são descritos, sendo eles: teste de Ljung-Box, e o teste deJarque-Bera.
No capítulo 8 estão as fórmulas para a previsão da volatilidade, medida diagnós-tico utilizada para a validação dos modelos propostos.
O capítulo 9 é inteiramente dedicado à aplicação e análise dos resultados do ajustedos modelos propostos às séries univariadas de retornos de extremos climatológicos.
Por fim, no capítulo 10 o trabalho é concluído destacando-se os pontos mais rele-vantes encontrados no capítulo anterior. É feita ainda uma análise crítica, concluindo-secom a proposta de trabalhos acadêmicos futuros, além da situação da implementação dométodo descrito no Instituto Nacional de Meteorologia (INMET).
Para os cálculos, foram utilizados alguns softwares estatísticos amplamente di-vulgados na área estatística, porém ainda pouco difundidos no campo metodológico dasciências atmosféricas:
• R: software livre disponibilizado em <http://www.r-project.org/>. Destaque espe-cial à ferramenta “Sweave” que possibilita a interface entre o R e o TEX, isto é,possibilita a migração de códigos e saídas decorrentes da programação, na escritado presente relatório.
• TEX( Miktex - <http://miktex.org/> e Winedt - <www.winedt.com/> ): Editor detexto usado para a construção deste relatório.
• Tinn-R: software livre disponibilizado em <http://sourceforge.net/projects/tinn-r/>:é um editor de texto para o ambiente R de fácil manuseio.
5 Fundamentação Teórica
5.1 Séries Temporais
Neste capítulo serão abordados alguns aspectos importantes na análise de sériestemporais. Tais aspectos são de grande utilidade para a compreensão dos capítulos subse-quentes. Os conceitos aqui usados foram retirados principalmente das referências ([14]) e([21]), respectivamente (Fuller, 1976) e (Brockell e Davis, 2002). Para melhor compreen-são, aconselha-se o uso de tais referências.
5.1.1 Aspectos gerais
Para se entender o conceito estatístico de série temporal, deve-se ter em mentealguns conceitos importantes. A idéia intuitiva que normalmente se tem de uma SérieTemporal (S.T.) é a de uma seqüência de observações que podem ser ordenadas no tempo.Entretanto, uma série temporal pode também ser vista como a realização de um processo
estocástico, conceito que será usado no presente trabalho.Um processo estocástico pode ser definido como uma coleção de variáveis aleatórias
Xt , t ∈ T , definidas em um espaço de probabilidades (Ω,A ,P). Aqui, T é o conjuntode índices que representa o espaçamento temporal entre duas variáveis aleatórias, geral-mente sendo um subconjunto dos números reais, isto é, T ⊆ R. O conjunto de todos osvalores S que cada uma das variáveis aleatórias pode assumir é chamado de espaço de
estados do processo.Para cada ω∈Ω, o conjunto formado pelos valores da coleção de variáveis aleatórias
do processo, X(ω),ω ∈Ω é chamado de realização ou função amostral do processo es-tocástico Xt , t ∈ T.
Uma realização pode ser vista como o “caminho por onde as observações reg-istradas percorreram”. O conjunto de todas as realizações possíveis de um processo es-tocástico é chamado de ensemble das realizações. Analogamente, o ensemble pode servisto como todos os possíveis caminhos por onde as observações podem ocorrer.
8
9
Agora, pode-se definir uma série temporal adequadamente: uma Série Temporal(S.T.) é um processo estocástico Xtt∈T , onde o T é caracterizado de duas formas:
• se T são intervalos: T = t : t1 < t < t2, diz-se que Xt é uma Série Temporal atempo contínuo (ou um processo a tempo contínuo).
• se T é um conjunto discreto: T = t1, t2, . . . , tn finito ou não, se diz que Xt é umaSérie Temporal a tempo discreto (ou um processo a tempo discreto).
Quanto ao tipo de variável aleatória, pode-se ter series temporais discretas: quandocada uma das variáveis aleatórias Xtt∈T é discreta, isto é, com espaço de estados Sdiscreto; e S.T. contínuas: quando S é contínua.
Dada essa definição, na prática quando se fala de uma série temporal, normalmenterefere-se à realização (ou parte dela) xt de uma S.T. Xt, pois pode ser bem complicadofazer várias medições de uma mesma variável a cada tempo.
No presente trabalho, serão consideradas Séries Temporais contínuas a tempo dis-creto, sendo T = Z
Na maioria das aplicações não há conhecimento da distribuição acumulada paraqualquer conjunto de variáveis,o que caracterizaria probabilisticamente as séries. Assim,é usual a utilização dos momentos das séries temporais, além da função de autocovariân-cia para tentar descrever a série.
Função de Autocovariância e Estacionariedade
A função de autocovariância extende o conceito de matriz de covariância paratratar de coleções infinitas de variáveis aleatórias, no caso um processo estocástico ouSérie Temporal com T infinito. Assim, a função de autocovariância de um processoXtt∈T onde Var(Xt) < ∞ é definida como:
γX(i, j) = Cov(Xi,X j) = E[(Xi−E(Xi))(X j−E(X j)
], ∀i, j ∈ T. (5.1)
onde a quantidade h = j− i é chamada de defasagem ou “lag”.Outro conceito importante no estudo de S.T. é o de estacionariedade. Em proces-
sos estacionários, propriedades estatísticas são dependentes apenas da defasagem da sériee não da variação temporal ocasionada pela defasagem. Diz-se que um processo “atingeo equilíbrio” quando o mesmo é estacionário.
Formalmente, uma série Xtt∈T é dita estacionária se:
1. E|Xt |2 < ∞ , ∀t ∈ T
10
2. E(Xt) = m , ∀t ∈ T
3. γX(i, j) = γX(i+ t, j + t) , ∀ i, j, t ∈ T .
Logo, para um processo estacionário, a função de autocovariância (5.1) pode serreescrita em função apenas da defasagem, isto é:
γX(i, j) = γX(i− i, j− i) = γX(0,h) = Cov(Xt ,Xt+h)≡ γX(h),∀ t,h ∈ T. (5.2)
Assim, em uma série estacionária a covariância entre valores defasados é depen-dente apenas da própria defasagem.
A função de autocovariância é importante para a comparação de propriedades bási-cas entre séries temporais, pois sumariza a estrutura de variância das séries. Porém, γX(h)é influenciada pela unidade de medida de cada um dos conjuntos de variáveis. Para re-solver isso, uma medida que não é influenciada por tal efeito é definida como função de
autocorrelação (FAC), sendo dada por:
ρX(i, j) =γX(i, j)√
γX(i, i) ·√
γX( j, j)=
Cov(Xi,X j)σi ·σ j
(5.3)
Sendo que, para séries estacionárias, é definida como:
ρX(h)≡ γX(h)γX(0)
,∀t,h ∈ T. (5.4)
Esta autocovariância padronizada serve para medir o comprimento e a memóriado processo, ou seja, a extensão para o qual o valor tomado no tempo t depende daqueletomado no tempo t−h.
A idéia da autocorrelação pode ser extendida tomando-se a correlação entre duasobservações da série, Xt e Xt+h, porém desconsiderando todos os termos intermediários,isto é, Xt+1,Xt+1, . . . ,Xt+h−1. A medida resultante é chamada de função de autocorre-
lação parcial (φhh), dada por:
φhh = ρX(t, t +h|t +1, t +2, . . . , t +h−1) (5.5)
Dada uma amostra aleatória de uma série temporal estacionária, isto é, o conjuntode observações x1,x2, . . . ,xn, pode-se estimar a função de autocovariância visando mel-hor conhecimento acerca da estrutura de dependência da série. Para isso, usa-se o esti-mador função de autocovariância amostral, definido por:
11
γX(h) =1n·
n−h
∑t=1
(xt+h− x)(xt− x) , 0≤ h < n (5.6)
onde x é a média amostral, dada por 1:
x =1n·
n
∑t=1
xi
Pode-se definir o estimador do coeficiente de correlação amostral (ρX(h)) em ter-mos do estimador acima, sendo dado por:
ρX(h) =γX(h)γX(0)
=
n−h
∑t=1
(xt+h− x)(xt− x)
n
∑t=1
(xt− x)2|h|< n (5.7)
Existem diversas maneiras para se estimar a função de autocorrelação parcial, paramaiores detalhes ver ([14]).
Um importante instrumento usado na identificação de modelos de séries temporaisé o Correlograma: gráfico dos coeficientes de autocorrelação ρX(h) versus a defasagemou “lag” h. Nos correlogramas amostrais são traçados limites paralelos ao eixo do númerode defasagens k (abcissa).
Tais limites representam os intervalos de confiança calculados de acordo com adistribuição assintótica do estimador ρX(h), e são usados para verificar se as autocorre-lações são significativamente diferentes de zero.
Há ainda o gráfico da função de autocorrelação parcial amostral (FACP), chamadode correlograma parcial.
A análise dos correlogramas tanto para a FAC quanto para FACP pode indicarmodelos potenciais que explicam determinados comportamentos característicos. Sãoainda muito utilizados na análise dos resíduos da série pois, dadas as suposições do mo-delo, espera-se que os resíduos se comportem de maneira aleatória e, dessa forma, ocorrelograma não pode aparentar padrões que indiquem relação temporal.
5.1.2 Séries Temporais Auto-Regressivas
Uma série temporal é dita autoregressiva quando o valor da variável observacionalé expressa em função de alguns de seus valores defasados no tempo. A equação abaixo,
1 a utilização de n ao invés de n-h se dá pelo caráter assintótico do estimador. Ver as referências citadas noinício do capítulo
12
chamada Equação de Diferenças, representa o modelo AR(p), um modelo de série autore-gresiva, com a variável observada dependente das p observações passadas.
p
∑i=0
αiXt−i = εt , t ∈ Z; (5.8)
com α0 6= 0, . . . ,αp 6= 0 e
εt ∼ N(0,σ2) , t ∈ Z ,σ2constante (5.9)
O conceito de autoregressão aqui explicitado é importante para o entendimento dametodologia relacionada aos modelos da família ARCH.
Um fato a se destacar é o comportamento característico do correlograma e do cor-relograma parcial de séries autoregressivas. Para tais séries, a função de autocorrelação(FAC) possui um decaimento para zero sob forma exponencial. Já a função de autocorre-lação parcial (FACP) possui um decaimento brusco para zero a partir de um “lag” k.
Determinar esse “lag” é, muitas vezes, uma tarefa difícil de ser realizada, porémé de extrema utilidade na medida que é um indicativo do número de defasagens p de ummodelo AR(p).
5.1.3 Estrutura dos dados - Os Retornos
Para os testes de ajustamento dos modelos serão trabalhados dados diários de al-gumas variáveis climatológicas. Em estudos práticos dos modelos GARCH em finançasé prática comum modelar diferenças entre termos defasados no tempo, também conheci-dos como “retornos” de tais variáveis. A preferência por retornos se dá por: além deserem livres de escala, suas propriedades estatísticas, como estacionariedade, são dese-jáveis para atender às suposições dos modelos, facilitando a manipulação estatística dosdados.
Dentre as diversas maneiras de se definir um retorno, será trabalhada a abordagemabaixo no presente trabalho. Sendo Xt o valor da variável climatológica no tempo t, osretornos de tais variáveis no tempo t são definidos por:
Rt = log(
Xt
Xt−1
)= log(Xt)− log(Xt−1) (5.10)
Espera-se, na modelagem dos retornos, que seja evidente a presença de “clustersde variação” (explicados mais a frente), sendo que tal característica, de variância variávelno tempo, é fundamental para a modelagem de séries segundo modelos GARCH.
13
Apesar das séries de retornos se mostrarem aparentemente simétricas, além dapresença de tal heterocedasticidade em clusters, observa-se ainda a presença de elevadacurtose (elevado grau de “achatamento” em relação à distribuição normal), e de leptocur-tose ou “cauda pesada” da distribuição dos retornos. Deste modo torna-se complicado asuposição de normalidade e independência da série de retornos.
Nota-se também que, se os retornos de fato são concisamente modelados porGARCH, os coeficientes de autocorrelação para os resíduos mostram-se baixos. Entre-tanto, valores absolutos dos retornos ou dos quadrados dos retornos mostra-se substan-cialmente elevados, isto é, indicando a presença de dependência em segunda ordem (nãouma dependência na média, mas sim na variância).
A utilização do retorno Rt (5.10) se dá devido ao caráter não linear de algumasséries de retornos. Apesar do uso do logaritmo, a razão é sustentada pelo fato de umasérie transformada preservar as características básicas da série original (ver ENDERS,1995).
Um cuidado a se ter é com relação à interpretação dos retornos na área climatólóg-ica, pois as análises não serão mais feitas com as variáveis originais. Quanto ao primeiroretorno, pode ser visto como o “ganho relativo” entre o tempo t e o tempo t +1. Para o pre-sente contexto, pode ser visto como a variação relativa de determinada variável climáticade um dia para o outro. Logo, no trabalho serão analisadas séries de ganhos relativos.
5.2 Modelo ARCH
5.2.1 Princípios
Modelos estatísticos são caracterizados pela existência de componentes aleatórias,sendo geralmente denominado um erro aleatório εt ao processo. E para a estimação deparâmetros de tais modelos, suposições são feitas em relação às diversas partes que ocompõe, sobretudo à distribuição desses erros.
A boa modelagem de tais erros é de extrema importância para a análise do modelo,pois dela são retiradas importantes informações sobre a distribuição dos dados, possibili-tando inferências de qualidade sobre os parâmetros.
Por características inerentes às variáveis em estudo, ou ainda para simplificaçãode cálculos, são feitas suposições sobre a distribuição de tais erros, como independênciae mesma distribuição, além de suposições sobre os momentos (ordinários), como médiaconstante e igual a zero, e ainda variância constante. Tais suposições conduzem à con-
14
strução de estimadores de tais parâmetros de forma robusta, com determinação analíticade características destes.
Porém, em muitas situações práticas como comumente em finanças, esta suposiçãode variância constante não atende à expectativas quando do ajustamento do modelo, pro-duzindo resultados não satisfatórios. Nesse contexto, modelos da família ARCH(AutoregressiveConditional Heteroscedaticity) são uma poderosa ferramenta, pois consideram a presençade grupos (“clusters”) de volatilidade 2.
5.2.2 Estrutura
O modelo ARCH foi inicialmente proposto por Engle ([11]) , com o intuito detratar sobre variações flutuantes em séries temporais, especificamente para modelar e pre-ver a inflação. Segundo ele, devido a tais flutuações a previsão da variabilidade torna-setão importante quanto a revisão dos valores da série original.
Precisamente, tal modelo é dito heterocedástico condicional, pois supõe que a var-iância ou volatilidade σ2
t (definida abaixo) de uma determinada variável em um instantet não se comporta de maneira constante. Supõe ainda que esse comportamento é depen-dente de informações desta e de outras variáveis em até certo instante passado 3. Portanto,o interesse não é modelar a série temporal em si, mas sim a volatilidade da mesma aolongo do tempo.
Deste modo, sendo ψt−1 o conjunto de todas as informações disponíveis e rele-vantes ao modelo até o momento t. e εt a série dos erros do modelo, define-se umavariância condicional σ2
t = Var(εt |ψt−1) , dependente de todas as informações até o in-stante imediatamente anterior t−1.
Assim, o modelo ARCH com q-defasagens ou ARCH(q) para a variável aleatóriaerro εt é dado por:
εt = σtzt , (5.11)
onde,
ztiid∼ N(0,1) , t ∈ Z , (5.12)
e
2 Por “clusters” de volatilidade, entende-se que há períodos de grandes erros e períodos de pequenos errosque ocorrem em grupos.
3 Do mesmo modo, é possível que a média varie com o tempo, ou outros momentos da distribuição domodelo variem com o tempo.
15
σ2t = Var(εt |ψt−1) = E[ε2
t |ψt−1] = α0 +q
∑i=1
αiε2t−i , (5.13)
sendo que α0 > 0 e αi ≥ 0 ∀i = 1,2, . . . ,q visando a não negatividade da var-iância condicional.
Antes de prosseguir, convém uma importante observação: embora Engle tenhaintroduzido o modelo utilizando a suposição de normalidade em (5.12), na prática outrasdistribuições podem ser utilizadas, de acordo com as características do modelo. Ressalta-se que, independente da distribuição utilizada, as suposições de média nula e variânciadependente temporalmente com expressão dada em (5.13) são o que caracteriza o modelo.
Com a estrutura acima, tem-se uma suposição sobre a distribuição dos erros condi-cionados ao passado, dada por:
εt |ψt−1 ∼ N(0,σ2t ) , t ∈ Z (5.14)
Observa-se ainda que a utilização de outras distribuições se deve ao fato de que asuposição de normalidade nem sempre é atendida em modelagens de séries financeiras,devido à caudas pesadas nas distribuições dos dados. Inúmeras outras distribuições podementão ser utilizadas para a modelagem. Aqui será dada importância ao caso da normal-idade, não descartando o uso futuro de outras distribuições, caso seja útil no ajuste dedados climatológicos.
5.2.3 Propriedades
Dada a construção do modelo, as propriedades dos modelos ARCH são mostradasabaixo. No apêndice encontra-se a demonstração de todas elas.
Tanto a esperança condicional, quanto a esperança incondicional são constantes eiguais a zero:
E[εt |ψt−1] = 0
e
E[εt ] = 0
Como propriedade fundamental, a variância condicional do modelo é dada por:
Var(εt |ψt−1) = σ2t
16
Já a variância incondicional, isto é, a variância não constante
Var(εt) =α0
1−∑qi=1 αi
, (5.15)
com 0 ≤ ∑qi=1 αi < 1 sendo a restrição de não-negatividade da variância incondi-
cional.Assumindo a estacionariedade do processo, mostra-se que a função de autoco-
variância é nula, isto é, a série dos erros é independente entre si. Tal característica éimportante na identificação e verificação do ajuste de modelos ARCH.
γε(h) = 0
As propriedades descritas acima caracterizam os modelos ARCH(q). Mostra-seainda que, sendo os erros εt pertencentes a um modelo ARCH(q), o quadrado dos mesmos,ε2
t , segue um modelo auto-regressivo com q defasagens no tempo AR(q), porém com errosnão-gaussianos ηt . Esses erros são variáveis aleatórias independentes, com média zero evariância não constante.
5.3 Modelo GARCH
5.3.1 Princípios
Evidências empíricas mostram que, para se estimar a volatilidade (σ2t ) de maneira
satisfatória e robusta usando o modelo ARCH(q), o número de defasagens (q) deve ser sat-isfatoriamente grande com pesos (parâmetros estimados) declinantes, acarretando assimum grande número de parâmetros a serem estimados pelo modelo.
Para solucionar esse problema, em 1986 Bollerslev criou uma nova versão do mo-delo que ficou conhecida como modelo GARCH (Generalized Autoregressive ConditionalHeteroscedaticity) ou GARCH (p,q). A idéia foi introduzir uma variável de média móvelpara deixar o modelo mais “parcimonioso”, no sentido de haver um número menor deparâmetros a se estimar.
5.3.2 Estrutura
17
A idéia é parecida com a extensão dos processos AR para processos ARMA,isto é, introduzir uma variável de média móvel no modelo. Nesse caso, a variânciacondicional do erro ou volatilidade (σ2
t ) depende não só dos quadrados dos erros pas-sados (ε2
t−1,ε2t−2, . . . ,ε
2t−q), mas também das próprias volatilidades defasadas no tempo
(σ2t−1,σ
2t−2, . . . ,σ
2t−p).
Assim, o modelo é descrito exatamente como no ARCH, isto é:
εt = σtzt , (5.16)
onde,
ztiid∼ N(0,1) , t ∈ Z , (5.17)
porém com uma pequena modificação em (5.13), sendo agora dada por:
σ2t = Var(εt |ψt−1) = E[ε2
t |ψt−1] = α0 +q
∑i=1
αiε2t−i +
p
∑j=1
β jσ2t− j , (5.18)
onde α0 > 0 , β j ≥ 0, ∀ j = 1,2, . . . , p e αi ≥ 0, ∀i = 1,2, . . . ,q visando a nãonegatividade da variância condicional.
Pela expressão acima, observa-se que modelos ARCH são casos especiais dosmodelos GARCH sendo que, especificamente, um modelo GARCH(0,q) é equivalente aum modelo ARCH(q).
GARCH, portanto, é um mecanismo que inclui as variâncias passadas na expli-cação das variâncias futuras. Mais especificamente, segundo Bollerslev (1986) ([5]),GARCH é uma técnica de séries temporais que permite utilizar o modelo de dependênciaserial da volatilidade.
O modelo, aparentemente mais complexo, é mais parcimonioso no sentido dehaver um número menor de parâmetros a ser estimado.
5.3.3 Propriedades
Na subseção de mesmo nome no capítulo anterior foram provadas as principaispropriedades dos modelos ARCH. Para os modelos GARCH mostra-se facilmente que osresultados encontrados são análogos, apenas mudando a expressão da variância condi-cional σ2
t pela suposta no modelo GARCH, dada em (5.18), isto é, considerando osquadrados das volatilidades anteriores.
18
Outra mudança se dá na expressão da variância incondicional do processo (σ2),dada por:
Var(εt) =α0
1−(
∑qi=1 αi +∑
pj=1 β j
) , (5.19)
sendo a restrição para a não-negatividade da variância incondicional dada por:k
∑i=1
(αi +βi) < 1 , onde k = maxp,q.
Analogamente ao mostrado nos modelos ARCH, modelos GARCH podem ser rep-resentados de forma conhecida, especificamente como um modelo ARMA (ver apêndice).
5.3.4 Limitações do Modelo GARCH
Abaixo são descritas algumas limitações ou problemas que são comumente en-contrados na utilização de modelos GARCH na modelagem de séries.
• Resultados satisfatórios são dependentes da própria capacidade de percepção (ex-periência) do usuário;
• GARCH é apenas parte de uma solução. As decisões financeiras são baseadas emoutros critérios além de retornos esperados e volatilidades;
• Freqüentemente falham por não capturar fenômenos altamente irregulares, incluindoflutuações ferozes do mercado e outros eventos altamente inesperados que podemconduzir a uma mudança estrutural significante;
• Comumente falham em não capturar o comportamento de curtose observadas emséries para pequenas amostras.
• A heteroscedasticidade consegue explicar alguns comportamentos de excesso decurtose, mas não todos;
5.3.5 Identificação de Modelos GARCH
Como visto, a variância condicional dos erros em modelos GARCH se comportacomo um processo autoregressivo e, desta forma, o correlograma é um instrumento bas-tante útil na identificação. Caso o modelo ajustado seja adequado, os correlogramas in-dicam a aleatoriedade o processo, enquanto o quadrado dos mesmos possui característicaautoregressivas.
Como visto, o correlograma de séries auto-regressivas possui um comportamentocaracterístico de decaimento exponencial para zero e, assim, na modelagem GARCHespera-se que o mesmo padrão seja encontrado nos quadrados dos resíduos.
5.4 Inferência Bayesiana
Em inferência estatística há o costume de se diferenciar os métodos de estimaçãode acordo com duas correntes de pensamentos: a “freqüentista” e a “bayesiana”. Aprimeira busca quantificar o valor do parâmetros populacionais (considerados fixos) deacordo com freqüências relativas ou proporções obtidas de uma amostra. Já a correntedita “bayesiana” baseia-se na informação obtida na amostra, porém combinada com umainformação adicional, subjetiva.
A subjetividade surge do fato da utilização de informações prévias acerca dosparâmetros populacionais, possivelmente relevantes, de um especialista naquela área doconhecimento tratada pelo problema. Tal medida “pessoal” da incerteza quantificadaprobabilisticamente, isto é, tem-se uma distribuição de probabilidade dos parâmetros con-siderados, portanto, aleatórios.
Para os freqüentistas tal informação não é válida, pois não pode ser observadadevido ao seu caráter subjetivo e, portanto, não é passível de verificação empírica.
Deste modo, a idéia básica da metodologia bayesiana é a tentativa da redução daincerteza sobre o parâmetro desconhecido, tal como sob a ótica freqüentistas, porém ointerpretando como uma variável aleatória. Para tanto, é utilizada a informação tantode um conhecimento prévio do parâmetro (denominado distribuição à priori) quanto doconhecimento obtido a partir de uma amostra (verossimilhança).
Para o presente contexto, será considerada ε = ε1,ε2, . . . ,εn uma amostra aleatóriade tamanho n, com função de verossimilhança dada mais a frente por: L(θ) = L(θ;ε) =f (ε|θ), onde f é a função de distribuição conjunta das variáveis aleatérias em questão.
Para os modelos GARCH, como visto anteriormente, denota-se por θ o conjuntode parâmetros do modelo, isto é:
θ = (θ0,θ1, . . . ,θp+q)T = (α0,α1, . . . ,αq,β1,β2, . . . ,βp)T ,
A quantidade θ assume valores em um conjunto denotado por Θ, conhecido comoespaço paramétrico. O tamanho de tal espaço é dependente do número de dimensões deθ, isto é, do número de parâmetros do modelo.
p(θ) = p(θ0,θ1, . . . ,θp+q) é dita a distribuição à priori conjunta dos parâmetrosdo modelo.
19
A regra de atualização que agrega as informações obtidas, tanto pelo conheci-mento do comportamento probabilístico dos parâmetros (priori) quanto da informaçãoobtida dos dados observados pela amostra (Verossimilhança), é denominada Teorema deBayes, dado por (Bayes, 1764):
p(θ|ε) =p(θ,ε)p(ε)
=f (ε|θ)p(θ)
p(ε)=
L(θ)p(θ)p(ε)
, (5.20)
onde π(θ) = p(θ|ε) é chamada distribuição a posteriori conjunta dos parâmetrosdo modelo. Pode-se observar que, como denominador da expressão que caracteriza adistribuição a posteriori, tem-se uma função que não traz informação alguma a respeitodo parâmetro θ, atuando como uma constante padronizadora.
p(ε) =
∫
p(ε|θ)p(θ)dθ, se θ é contínuo;
∑ p(ε|θ)p(θ), se θ é discreto;
Neste trabalho será utilizada a notação para o caso contínuo. Lançando mão daconstante normalizadora, obtém-se a forma usual do Teorema de Bayes, expressa por:
π(θ) ∝ L(θ)p(θ) , (5.21)
onde ∝ é o símbolo de proporcionalidade.
5.4.1 Conceitos importantes
Além dos conceitos fundamentais esboçados na introdução, outros termos impor-tantes para o entendimento do método estão descritos a seguir.
Distribuição à Priori
Como dito anteriormente, tal distribuição é a chave fundamental para a diferenci-ação entre a metodologia clássica e a bayesiana. Se refere à distribuição dos parâmetrospopulacionais de interesse, pois estes são considerados variáveis aleatórias.
Como quaisquer outras distribuições de probabilidades, as distribuições à prioripodem conter parâmetros, sendo estes chamados de “hiperparâmetros”. Ao contrário dosparâmetros de interesse, são considerados fixos e normalmente julga-se conhecidos osvalores dos mesmos para o estimativa da posteriori.
Para os modelos GARCH, devido às restrições de não-negatividade dos parâmet-
ros θi ≥ 0 e de quep+q
∑i=1
(θi) < 1, algumas prioris já foram exploradas na literatura (vide
[2] e [25] para exemplos).
20
21
Como desempenha um papel fundamental na metodologia bayesiana, deve-se tomarcuidado com a escolha de tal distribuição. Os “tipos” de distribuição a priori são basica-mente divididos em: Prioris informativa e não-informativa; Pioris Conjugadas; e PriorisHierárquicas.
• Priori informativa: quando há ganho substancial de informação com seu conhec-imento;
• Priori não-informativa: quando não há quaisquer conhecimentos a cerca da dis-tribuição dos parâmetros;
• Prioris Conjugadas: quando tanto a priori quanto a posteriori pertencem à mesmafamília de distribuições;
• Prioris Hierárquicas:
Distribuição à Posteriori
De acordo com Gelmam (2002), a distribuição à posteriori resume o estado atualde conhecimento sobre todas as quantidades incertas (incluindo parâmetros não observáveis,bem como dados potenciais não observáveis, latentes ou ausentes) na análise bayesiana.Isto é, é a “revisão” de conhecimento prévio (priori) de acordo com a informação trazidapela amostra (verossimilhança).
A diferença entre a probabilidade à priori e à posteriori caracteriza a informaçãoobtida pelo experimento.
Para o caso do modelo GARCH, a distribuição à posteriori não possui forma con-hecida quando se supõe as prioris sugeridas em estudos prévios. E devido à sua com-plexidade, métodos de simulação são exigidos para a estimação dos parâmetros, comodescritos no próximo capítulo.
Distribuição Marginal à Posteriori
A distribuição a posteriori π(θ) do conjunto de parâmetros θ =(θ0,θ1, . . . ,θq+p)T =(α0,α1, . . . ,αq,β1,β2, . . . ,βp)T , foi especificada acima, pela equação (5.20). Porém, emaplicações práticas pode se estar interessado na distribuição marginal do vetor de parâmet-ros, isto é, na distribuição de um parâmetro específico, digamos θi. Assim, a distribuiçãomarginal é derivada da conjunta, a partir da expressão:
π(θi) = p(θi|ε) =∫
π(θ0,θ1, . . . ,θq+p)dθ−i (5.22)
onde θ−i = (θ0, . . . ,θi−1,θi+1θq+p).
22
Distribuição Marginal ou Preditiva à Priori
É a função ou constante normalizadora utilizada na construção da posteriori, dadaacima em (5.20). É dita “preditiva”, pois utiliza uma quantidade observável e “à priori”,e não é condicional à nenhuma observação prévia dos dados, isto é, não depende deum conhecimento prévio, já sendo um. Expressa uma idéia da distribuição esperada davariável sendo modelada, dado o conhecimento do parâmetro, isto é:
p(ε) =∫
p(ε|θ)p(θ)dθ = E [p(ε|θ)] (5.23)
É importante na medida que fornece um modo de checagem da adequação dadistribuição à priori, ao fazer predições via distribuição preditiva à priori. Caso as ob-servações (amostra) observada tivesse pouca probabilidade preditiva, deve-se questionaro modelo. Além disso, tal distribuição é útil para a escolha de modelos bayesianos, nãosendo tão utilizada para as predições. (Congdon, 2003).
Distribuição Preditiva à Posteriori
É utilizada quando há o interesse em uma atualização, dada pela previsão de novosvalores da distribuição dos dados. Sendo y o novo valor a ser predito, deseja-se obter adistribuição de probabilidades de y|ε, onde ε é a toda a informação amostral disponívelaté então, isto é, o conjunto de dados observado.
Para tanto, utiliza-se a distribuição à posteriori da seguinte forma:
p(y|ε) =∫
p(y,θ|ε)dθ =∫
p(y|θ,ε)p(θ|ε)dθ (5.24)
*=∫
p(y|θ)p(θ|ε)dθ = E [p(y|θ)] (5.25)
A última igualdade ocorre quando há suposição de independência entre os dadosobservados e o novo valor a ser predito, dado o conhecimento do parâmetro.
Neste trabalho será feita análise da distribuição preditiva à posteriori das volatili-dades, objetivo de previsões.
Distribuição Condicional à Posteriori
π(θi|θ j) =π(θi,θ j)
π(θ j),∀i (5.26)
onde j ∈ 0,1, . . . , i−1, i+1, . . . , p+q
Como caso particular, tem-se a distribuição condicional completa de cada um dosparâmetros, designada por:
π(θi|θ−i) =π(θi,θ−i)
π(θ−i)=
π(θ)∫π(θ)dθi
,∀i (5.27)
onde θ−i = (θ0, . . . ,θi−1,θi+1θq+p).O conhecimento de tais distribuições é a chave para um poderoso algoritmo de
simulação: o amostrador de Gibbs, demonstrado adiante. Porém, em muitas ocasiões, aresolução analítica de tais integrais é impossível dada tamanha complexidade. Nesse casooutros métodos de simulação como: a amostragem por importância (Paulino et al., 2003),o slice sampling proposto por Neal (2003), a amostragem por aceitação e rejeição e oMetropolis-Hastings (Chib e Greenberg, 1995; Hastings, 1970), podem ser usados para asimulação de amostras da posteriori.
5.4.2 Inferência
Toda inferência será feita com base na distribuição à posteriori, de onde se obtémas estatísticas necessárias para resumir o comportamento dos parâmetros. Dentre as prin-cipais estatísticas a posteriori pode-se citar:
Média da posteriori (Estimador de Bayes)
É resultado da aplicação da teoria bayesiana de decisão para encontrar estimadorespontuais de interesse, quando há uma função quadrática de perda (ver Bolfarine, 2003em [4]). Os estimadores de Bayes não são invariantes como os estimadores de MáximaVerossimilhança, isto é, se θ é estimador de Bayes de θ, não necessariamente tem-se queg(θ) é estimador de Bayes de g(θ); Apesar de sempre ser função de estatística suficientepara o parâmetro, o estimador de Bayes não é necessariamente não viesado. Porém talviés diminui assintoticamente.
E(θ|ε) =∫
θ∈Θ
θp(θ|ε)dθ (5.28)
Quantil α à posteriori
É utilizado para a construção dos intervalos de confiança dos parâmetros do mo-delo. Um fato interessante a se destacar é a interpretação dos IC’s segundo a metodologiabayesiana.
23
24
Q(α) =
θ′∈Θ :
∫θ′
−∞
p(θ|ε)dθ = α
(5.29)
Intervalo 100(1−α)% de credibilidade a posteriori , (L,U)
Enquanto na linha frequentista se diz, após a retirada de inúmeras amostras dadistribuição, o verdadeiro valor do parâmetro encontra-se no interior do intervalos em1−α % das ocorrências, o IC’s bayesianos não fazem suposição a respeito de retiradashipotéticas de amostras, mas sim é calculado um intervalo “de credibilidade”, isto é, aondeespera-se que o parâmetro assuma valores em tais intervalos 1−α% das vezes.
(L,U) =
(L′,U ′)⊂Θ2 :∫ U ′
L′p(θ|ε)dθ = 1−α
,α ∈ (0,1); (5.30)
Quando o intervalo é simétrico, tem-se que L = Q(α/2) e U = Q(1−α/2). Quantomenor for o tamanho do intervalo mais concentrada é a distribuição do parâmetro, ou sejao tamanho do intervalo informa sobre a dispersão de θ.
Assim como em IC’s para inferência clássica, os intervalos de credibilidade bayesianossão invariantes a transformações 1 a 1 nos parâmetros, isto é, se o IC para α é (L,U), o ICpara g(α) é (g(L),g(U)).
Observa-se que em tais distribuições, que caracterizam a inferência bayesiana, háconstante presença de integrais que nem sempre são resolvidas analiticamente em prob-lemas reais. Nesses casos são requeridos métodos computacionalmente intensivos paraa resolução das mesmas. Alguns desses métodos estão no escopo deste trabalho e serãodiscutidos mais à frente.
5.5 Simulação
Usado em diversos ramos da ciência, os métodos de Monte Carlo são técnicasestocásticas baseadas em simulação. Tipicamente, envolve a geração aleatória de obser-vações de alguma distribuição de probabilidades, sendo tal amostra utilizada posterior-mente na obtenção de aproximações numéricas de funções complexas.
No presente contexto, há o interesse na aproximação numérica de integrais que nãopossuem resolução analítica, sendo então discutidos aqui alguns métodos de Integraçãode Monte Carlo.
25
5.5.1 Métodos de Integração de Monte Carlo
A idéia básica do Método é a escrita de integrais 4 de interesse, geralmente semresolução analítica, como valores esperados. Deste modo, sendo I a integral a ser aproxi-mada, geralmente sem solução analítica, escreve-se:
I =∫ b
ag(x)p(x)dx = Ep(x)[g(x)] (5.31)
Isto é, a integral é reescrita como sendo o valor esperado de uma função (g(X)),com respeito a uma variável aleatória (X) com densidade p(x) de onde se saiba simu-lar valores. Deste modo, simulando n valores (x1, . . . ,xn) de p(x), estima-se o valor daintegral pela expressão:
I = gn(x) =1n
n
∑i=1
g(xi) (5.32)
A equação acima define o chamado estimador de Monte Carlo Simples. Dentreas diversas aplicações, pode-se citar o cálculo das densidades marginais completas vis-tas anteriormente em (5.27). Observa-se que tal estimador nada mais é que uma médiaamostral de g(x) usado na aproximação de um valor esperado (média teórica).
Porém, nem sempre é possível gerar valores da distribuição p(x), tendo portantoproblemas com este estimador. Nesses casos, uma das possíveis alternativas é usar umafunção auxiliar, por exemplo q(x), da qual consegue-se simular valores mais facilmente.Assim, a integral de interesse (5.31) é reescrita como:
I =∫ b
ag(x)p(x)dx =
∫ b
a
g(x)p(x)q(x)
q(x)dx = Eq(x)
[g(x)p(x)
q(x)
](5.33)
Analogamente, após simular n valores da distribuição auxiliar q(x), calcula-se oestimador de Monte Carlo via Função de Importância, dado por:
I = gqn(x) =
1n
n
∑i=1
g(xi)p(xi)q(xi)
(5.34)
Para os dois estimadores ao definir σ2g como a variância de g(x), no caso do
primeiro estimador, e a variância de [g(x)p(x)/q(x)] para o segundo estimador 5, pode-secalcular uma estimativa consistente para tais variâncias nos dois casos. Tal estimativa édenominada Erro Padrão de Monte Carlo, dada por:
4 Para o caso discreto, não tratado no presente contexto, o interesse seria em somatórios5 Usando a desigualdade de Cauchy-schwarz, mostra-se que a variância deste estimador é minimizada
quando q(x) ∝ |g(x)|(Rubinstein,1981)
26
σg =
√1n2
n
∑i=1
(h(xi)− h)2 (5.35)
onde h(x) = g(x) para o estimador simples e h(x) = [g(x)p(x)/q(x)], para o esti-mador via função de importância.
Os dois estimadores, ambos gerados via método dos momentos, gozam das mes-mas propriedades assintóticos (Gamerman, 1997):
1. São não viesados para a estimação de I;
2. Por serem estimadores do tipo média amostrais, usando a Lei Forte dos GrandesNúmeros, assegura-se que a média amostral gn(x)egq
n(x) converge quase certamente(q.c) para a média teórica I, isto é, são estimadores fortemente consistentes para I:
Iq.c→ I , quando n→ ∞
3. Pelo Teorema do Limite Central(TLC) tem-se um importante resultado assintóticoda convergência em distribuição dos dois estimadores dada abaixo. Nota-se quecom tal resultado pode-se construir intervalos de confiança para a integral I;
I− IσI
TLC→ N(0,1) , quando n→ ∞
onde σI é Erro padrão de Monte Carlo para o determinado estimador de I;
No contexto Bayesiano, um problema já mencionado é o conhecimento da con-stante de proporcionalidade dada em (5.23). Quando a Integral de interesse é algumaesperança calculada em relação à distribuição à posteriori, é inevitável o cálculo da con-stante normalizadora e, nesse caso, usa-se outro método para a resolução de tais integrais.Na verdade, o método consiste do uso de um estimador composto pela razão de doisestimadores Monte Carlo de integrais.
Sendo π(θ) uma densidade à posteriori, e∫
g(θ)π(θ)θ uma função de interesse,tem-se a equação abaixo explicitando a constante normalizadora:
I =∫
g(θ)π(θ)dθ =∫
g(θ)l(θ)p(θ)dθ∫l(θ)p(θ)dθ
Assim, sendo q(θ) uma função de importância, aproxima-se a integral superiorsegundo:
I = ∑ni=1 g(θi)π(θi)/q(θi)∑
ni=1 π(θi)/q(θi)
27
O estimado acima, apesar de assintoticamente não-viesado, é um estimador forte-mente consistente para I.
Há ainda alguns métodos baseados em simulação estocásticas, de suma importân-cia na evolução do método MCMC, descritos melhor no apêndice deste trabalho.
5.6 Cadeias de Markov
5.6.1 Conceitos Importantes
Como definido anteriormente, um processo estocástico é uma sequência de var-iáveis aleatórias Xt , t ∈ T , definidas em um espaço de probabilidades (Ω,A ,P), ondeT é o conjunto de todos os índices do processo. Devido ao caráter estrutural das sériesem estudo, foi assumido T finito ou infinito enumerável. Foi chamado de S o conjunto detodos os valores possíveis do processo ou espaço de estados.
Uma Cadeia de Markov é um caso especial de processo estocástico, onde os es-tados são governados por uma probabilidade de transição que satisfaz a propriedade deMarkov, isto é:
P(Xt |Xt−1,Xt−2, . . . ,X0) = P(Xt |Xt−1)
Em outras palavras, dado o estado presente da cadeia Xt , os estados passados nãotem influência alguma sobre o futuro. Mais geralmente pode ser escrita, sendo t0 < t1 . . . <
tn−1 < tn < tn+1, como:
P(Xtn+1 = in+1|Xtn = in, . . . ,Xt1 = i1,Xt0 = i0) = P(Xtn+1 = in+1|Xtn = in) (5.36)
Assim, define-se a função ou Kernel de transição da cadeia como:
Pn(i, j) = P(Xtn+1 = j|Xtn = i) , ∀i, j ∈ S (5.37)
Serão estudadas apenas cadeias homogêneas, isto é, um caso particular onde afunção de transição da cadeia não depende de n. Nesse caso, Pn(i, j) = Pi, j,∀n é a proba-bilidade de transição a um passo da cadeia.
Quando o S é discreto, define-se a Matriz de Probabilidades de Transição de umpasso ou Matriz de Transição ou ainda Matriz Estocástica como sendo: P = [Pi, j]i, j∈S, isto
28
é, e a matriz formanda pelas probabilidades de transição de todos os estados da cadeia.Nota-se que:
• Pi, j ≥ 0 , ∀i, j ∈ S
• ∑ j∈S Pi, j = 1 , ∀i ∈ S
Será denotado por π(0) a distribuição inicial de uma cadeia, isto é:
π(0)i = P(X0 = i) , ∀i, j ∈ S (5.38)
As Probabilidades de Transição de n-passos, isto é, a probabilidade de, partindodo estado “i‘”, chegar ao estado “j” em exatas “n” transições, são definidas da seguinteforma:
P(n)i, j = P(Xm+n = j|Xn = i) = P(Xn = j|X0 = i) (5.39)
Sendo assim, Chapman e Kolmogorov desenvolveram independentemente equaçõesúteis para o cálculo de transições diversas, conhecidas como equações de Chapman-Kolmogorov, descritas por:
P(n+m)i, j = ∑
k∈SP(n)
i,k P(m)i, j (5.40)
ou, na notação matricial, sendo Pn a n-ésima potência da Matriz de Transições:
P(n+m) = P(n)P(m)
Assim, com a distribuição inicial da cadeia (5.38) e com a Matriz de Transição deum passo (5.37), a distribuição parcial do n-ésimo estado da cadeia, isto é, a distribuiçãode probabilidades após n passos da cadeia, pode ser acessada segundo:
P(Xn = j) = ∑i∈S
P(Xn = j|X0 = i)P(X0 = i) = ∑i∈S
P(n)i, j π
(0)i (5.41)
Quando existe n tal que P(n)i, j > 0 diz-se que o estado “j” é acessível a partir de
“i”. Se j é acessível a partir de i e vice-versa, diz-se que os estados “i” e “j” comunicam-se entre si. Deste modo, define-se uma propriedade desejável para as cadeias a seremestudadas, a irredutibilidade: uma cadeia é dita irredutível quando todos os estados dacadeia comunicam-se entre si.
Definindo-se fi, j como sendo a probabilidade de, partindo do estado “i”, visitaro estado “j” pela primeira vez em algum tempo finito, classificamos um estado comorecorrente quando fi, j = 1; ou transiente, quando fi, j < 1.
29
Outra propriedade indispensável para a construção da teoria é a periodicidade deuma cadeia. Um estado é chamado de absorvente quando, partindo dele, só chegamos nopróprio: Pi,i = 1.
Deste modo, uma cadeia é dita aperiódica quando não possui nenhum estado ab-sorvente. Nesses casos, evita-se que, caso a cadeia “pare” em um estado absorvente,“fique” nele para sempre, não visitando outros pontos.
5.6.2 Probabilidades Limites e Distribuições Estacionárias
Em problemas de simulação há interesse no comportamento da cadeia após umnúmero n muito grande de iterações, isto é, deseja-se saber a distribuição assintótica deP(n)
i, j , quando n→∞. Para um número muito grande de iterações, P(n)i, j converge a um valor
π j > 0 dependente somente do próprio estado “j”, restando-nos apenas saber caracterizartal distribuição.
Mas antes da caracterização de Distribuições Limites, será definido o conceitode estacionariedade de cadeias. Seja
π j, j ∈ S
uma distribuição de probabilidades,
satisfazendo, unicamente, o seguinte sistema:Uma cadeia satisfazendo tais condições é chamada de Distribuição Estacionária
ou de Equilíbrio da Cadeia. Pode-se mostrar (ver Gamerman, 1997) que, sendo Xtt∈T
uma cadeia de Markov irredutível e aperiódica e com S finito, ∀ j ∈ S,
limn→∞
Pni, j = π j > 0 , ∀ i ∈ S (5.42)
onde π j é a distribuição estacionária de Xtt∈T . Nessas condições, tal cadeia échamada ergódica e
π j, j ∈ S
é a distribuição Limite da Cadeia.
Nota-se que a existência de distribuição estacionária de uma cadeia não implicanecessariamente na existência de distribuição limite. Mas, se a distribuição limite existee satisfaz (5.42) então, independentemente do estado inicial da cadeia, quando n→ ∞, acadeia tende a π.
Uma cadeia homogênea é dita reversível quando, estando o sistema em equi-líbrio, dado dois estados quaisquer, a taxa de movimentação do sistema de θ′ para θ éa mesma para a movimentação de θ para θ′. Esta propriedade é melhor explicada se-gundo a equação abaixo, denominada equação detalhada de balanceamento ou condiçãode reversibilidade.
π(θ)K(θ,θ′) = π(θ′)K(θ′,θ) , ∀(θ,θ′) (5.43)
30
Agora, o resultado mais importante para a conexão das duas teorias é estabelecido,conhecido como Teorema ergódico (ver [19] )
Considere uma distribuição π(θ) = π(θ1, . . . ,θp) conhecida, a menos de uma con-stante multiplicativa, porém da qual não consegue-se simular valores. Seja agora
θ(t), t ∈ T
uma Cadeia de Markov ergódica com distribuição de equilíbrio π. Então, pelo teoremaergódico:
θ(t) n→∞→ π(θ)
e
g =1n
n
∑t=1
g(θ(t)i ) n→∞→ Eπ(g(θi)) , quase certamente.
isto é, médias de valores da cadeia produzem estimadores fortemente consistentesdos parâmetros da distribuição limite. Embora a cadeia seja por definição dependentetemporalmente, a média aritmética dos valores da cadeia formam um estimador consis-tente da média teórica.
Tal resultado para cadeias de Markov é equivalente à lei dos grandes númerosmencionada anteriormente (Ver seção “Métodos de Integração de Monte Carlo”) para osestimadores de Monte Carlo. Analogamente, há um importante resultado para cadeias deMarkov equivalente ao teorema do limite central.
Considerando, por hipótese, σ2g a variância de g(θ) finita, a variância do estimador
de g é dado por σ2g/n e então o Teorema do Limite Central garante a convergência em
distribuição do erro:
√n
(g−Eπ[g(θ)]
σ2g
)d→ N(0,1) (5.44)
Para a estimativa da variância σ2g, ao considerar a dependência temporal dos val-
ores gerados no processo de amostragem, pode-se mostrar que o erro padrão é dado por:
Var (g) =S2
n
[1+2
n−1
∑k=1
(1− k
n
)ρk
]onde S2 é a variância amostral e ρk a autocorrelação amostral de ordem k dada
anteriormente em (5.4).
31
5.6.3 Espaço de Estados Contínuo
A idéia básica construída até aqui pode ser facilmente estendida para o caso con-tínuo, isto é, quando o espaço de estados S é contínuo. Este é o tipo de cadeias que nopresente trabalho se tem interesse.
O kernel ou função de transição do caso discreto (5.37) é substituído pelo Kernelde probabilidade P(x,y) , satisfazendo:
∫P(x,y)dy = 1 , ∀x ∈ S (5.45)
Já a extensão contínua das equações de Chapman-Kolmogorov é dada por:
πt =∫
y∈Sπt−1(x)P(x,y) , ∀y ∈ S (5.46)
Por fim, a distribuição estacionária da cadeia é dada após a resolução da equação:
π∗ =
∫y∈S
π∗(x)P(x,y) , ∀y ∈ S (5.47)
6 Estimação de Modelos GARCH
Dentre os métodos frequentistas mais utilizados, destaca-se o Método da MáximaVerossimilhança na estimação dos parâmetros do modelo GARCH, inclusive indicadoinicialmente por Engle ([11]). Outros métodos de estimação foram utilizados na lite-ratura para a estimação dos parâmetros de modelos ARCH apenas, como o Método daQuase-Verossimilhança e ainda o método dos momentos generalizado. Para detalhes verBollerslev, Chou e Kroner(1992) ([7])
Será aqui demonstrado o método para os modelos GARCH em geral. Os resulta-dos para os modelos ARCH são facilmente encontrados fazendo-se p = 0 nas fórmulassubsequentes.
Seja ε1,ε2, . . . ,εn uma amostra aleatória (ou realização) de tamanho n da famíliade variáveis aleatórias εtt∈Z com função densidade de probabilidade f (εi|θ) , i =1,2, . . . ,q.
Aqui, θ é o vetor de parâmetros do modelo GARCH(p,q), dado por
θ = (α0,α1, . . . ,αq,β1,β2, . . . ,βp)T ,
tem-se que a Função de Verossimilhança de θ correspondente à amostra aleatóriaobservada é dada por:
L(θ;ε1,ε2, . . . ,εn) = f (ε1,ε2, . . . ,εn|θ) =n
∏t=1
f (εt |θ) (6.1)
Define-se ainda o Estimador de Máxima Verossimilhança de θ como sendo o valorθ que maximiza a Função de Verossimilhança.
Uma função útil para o estudo de tal técnica é a função de Log-Verossimilhança,dada ao se aplicar uma linearização da função de Verossimilhança por meio do logaritmonatural ou neperiano (ln), denotado aqui por log, isto é:
l(θ;ε1,ε2, . . . ,εn) = logL(θ;ε1,ε2, . . . ,εn) (6.2)
32
33
Verifica-se facilmente que o valor de θ que maximiza a função de Verossimilhançatambém maximiza a Função de Log-Verossimilhança. Na maioria das situações é maissimples manipular funções de Log-Verossimilhança, facilitando assim os cálculos para abusca do estimador.
Fato de grande importância quando se utiliza o condicionamento é a decomposiçãoda Verossimilhança. Como mostrado a seguir, qualquer Verossimilhança pode ser decom-posta em função de suas densidades condicionais. Sendo L(θ) = f (ε1, . . . ,εq|θ) a funçãode Verossimilhança exata das q primeiras variáveis aleatórias da amostra, podemos escr-ever a função de verossimilhança da amostra retirada:
L(ε1,ε2, . . . ,εn|θ) = f (ε1, . . . ,εq,εq+1, . . . ,εn|θ)
= f (ε1, . . . ,εq|θ) · f (εn, . . . ,εq+1|ε1, . . . ,εq,θ)
= L(θ) · f (εn, . . . ,εq+2|ε1, . . . ,εq+1,θ) · f (εq+1|ε1, . . . ,εq,θ)...
= L(θ) ·n
∏t=q+1
f (εt |ε1, . . . ,εt−1,θ)
= L(θ) ·Lc(θ) ,
onde Lc(θ) é a função de Verossimilhança condicional conjunta de ε1, . . . ,εn “dadoo passado”, isto é, εt |ψt−1, t = q+1,q+2, . . . ,n
Para um número grande de observações (n) em relação à defasagem (q), pode-se“desprezar” a função de Verossimilhança exata e trabalhar apenas com a Verossimilhançacondicional, considerando a suposição inicial do modelo sobre a distribuição dos erroscondicionais.
Sendo Lt(θ) = f (εt |ε1, . . . ,εt−1,θ) a densidade estacionária marginal de εt , pode-se escrever:
L(θ;ε1,ε2, . . . ,εn) = L(θ) ·Lc(θ)
= L(θ) ·n
∏t=q+1
Lt(θ)
∝
n
∏t=q+1
Lt(θ)
Devido à suposição de normalidade dos erros condicionados às informações pas-sadas (5.14), pode-se escrever:
34
Lt(θ) = f (εt |ψt−1) = (2π)−12 ·(σ
2t)− 1
2 · exp
12· ε2
t
σ2t
(6.3)
e assim,
Lc(θ) =n
∏t=q+1
Lt(θ) =n
∏t=q+1
[(2π)−
12 ·(σ
2t)− 1
2 · exp
12· y2
t
σ2t
](6.4)
Aplicando a função de Log-Verossimilhança, chega-se na seguinte expressão:
lc(θ) = logLc(θ)
= logn
∏t=q+1
Lt(θ)
=n
∑t=q+1
logLt(θ)
=n
∑t=q+1
lt(θ) ,
com,
lt(θ) =−12· log(2π)− 1
2· log
(σ
2t)+
ε2t
2σ2t
, (6.5)
e, portanto,
lc(θ) =n
∑t=q+1
[−1
2· log(2π)− 1
2· log
(σ
2t)+
ε2t
2σ2t
]Como observado, devido ao caráter não-linear da função de Log-verossimilhança
em θ = (α0,α1, . . . ,αq,β1,β2, . . . ,βp)T , a maximização da mesma, dada pela solução de∇lt(θ) = 0, onde ∇lt(θ) é o vetor score ou vetor de derivadas parciais de primeira ordemda Função de Log-Verossimilhança, constitui um problema de otimização não-linear edeve ser resolvido por métodos numéricos.
A maioria dos métodos involve o cálculo tanto do vetor Score, quanto da matrizHessiana (ou matriz de derivadas segundas) da Função de Verossimilhança. Tais Cálculosestão explicitados no apêndice G.
Considerando as suposições sobre zt de nulidade da esperança condicional e davariância condicional igual à 1, além da estacionariedade estrita dos resíduos εt e outrascondições técnicas, o estimador de máxima verrosimilhança é consistente (Engle, 1996).
Outro método bastante útil no contexto dos modelos GARCH é a estimação viaPseudo-Máxima Verossimilhança (ver White, 1982 e Weiss, 1982). Tal método aproxima
35
uma distribuição de probabilidades de variáveis aleatórias, supondo conhecidos e finitosalguns momentos populacionais da distribuição. Em ([22]) o autor demonstra todas aspropriedades de tal método, porém para o modelo infinito.
Em ([18]), Greene (2007) analisa o efeito da não-normalidade utilizando-se apseuda-verossimilhança. Um importante resultado é que o estimador obtido maximizando-se a mesma função de log-verossimilhança convencional é consistente, porém a matriz decovariâncias dos estimadores deve ser corrigida.
6.1 Método Clássico: Otimização de Funções
Não-Lineares
Segundo os modelos GARCH(p,q), considerando o método clássico de estimaçãovia máxima verossimilhança, a função de Log-Verossimilhança, como visto anterior-mente, não possui caráter linear. Assim, as equações de primeira ordem, que são condiçõessuficientes para a maximização da função de Log-Verossimilhança, dadas ao igualar afunção Score a zero: ∇lt(θ) = 0, não podem ser resolvidas analiticamente.
Nessa seção serão comentados os métodos mais utilizados para a resolução doproblema e estimação dos parâmetros.
6.1.1 Método de Newton-Raphson
É um método iterativo usado para maximizar a função de verossimilhança maisfacilmente. Os estimadores obtidos por este processo estarão bem próximos dos esti-madores de máxima verossimilhança.
Dado o vetor atual de estimativas θ(k) = (θ(k)0 ,θ
(k)1 , . . . ,θ
(k)p+q+1)
T dos parâmetrosθ(k), melhora-se essa estimativa através da busca do ponto máximo da expansão de Taylorda função de Log-verossimilhança lt(θ), em torno de θ(k). Para tanto, tal função deve sercontínua, diferenciável e com derivadas de segunda ordem não-nulas.
Como resultado, tem-se o seguinte processo iterativo:
θ(k+1) = θ
(k)−H−1(
θ(k))·∇lt(θ) (6.6)
Assim, dado um vetor de estimativas inicial θ(0), basta executar o processo itera-tivamente até que, dado um ε pequeno, a seguinte condição esteja satisfeita:
36
‖θ(k+1)− θ(k)‖‖θ(k)‖
< ε , (6.7)
onde ‖θ(k)‖ denota o tamanho ou norma do k-ésimo vetor de estimativas θ(k) =(θ(k)
0 ,θ(k)1 , . . . ,θ
(k)p +q+1)T , dado por:
‖θ(k)‖=√
(θ(k)0 )2 +(θ(k)
0 )2 + . . .+(θ(k)p+q+1)2
6.1.2 Método Scoring de Fisher
É o método iterativo mais utilizado para ajustar Modelos Lineares Generaliza-dos (MLG). É, baseado no processo iterativo de Newton-Raphson, porém usando-se umaaproximação da matriz Hessiana ou matriz de derivadas segundas da função de Log-Verossimilhança H(), facilitando assim os cálculos.
A aproximação é feita pela matriz de Informação de Fisher, dada pelo valor nega-tivo do valor esperado da matriz de derivadas de segunda ordem do logaritmo da funçãode Log-verossimilhança, isto é:
I(θ) =−E [H(θ)] (6.8)
Agora, aproximando a matriz de Informação pela Matriz Hessiana, isto é I(θ) ≈H(θ) a equação 6.6 pode ser reescrita da seguinte maneira:
θ(k+1) = θ
(k)− I(θ)−1 ·∇lt(θ) , (6.9)
sendo o critério de parada do algoritmo o mesmo usado no método anterior.Uma importante observação a ser feita é que, em tais métodos, a convergência
não é garantida. Há sempre necessidade de verificação da convergência dos valores quedepende da proximidade que os valores iniciais escolhidos na iteração estejam em relaçãoaos verdadeiros valores dos parâmetros e ainda do tipo de função não-linear que estásendo maximizada.
6.1.3 Métodos Quase-Newton
Os Métodos Quase-Newton possuem como proposta inicial evitar o cálculo damatriz Hessiana requerido em (6.6) a cada iteração, aproximando-a por uma matriz “B”e assim reduzindo o número de cálculos a cada iteração. A intenção é manter ao máximo
37
as propriedades de convergência do Método de Newton-Raphson. Os métodos são entãodescritos como:
θ(k+1) = θ
(k)−B−1 ·∇lt(θ) (6.10)
Em geral, tais métodos conseguem obter a solução em um tempo de execuçãomenor e com menor esforço computacional comparados aos métodos anteriores.
Para cada tipo de escolha da matriz de atualização “B” há um novo método. Dentreos métodos existentes, destacam-se o de Broyden-Fletcher-Goldfarb-Shanno (BFGS), oInverse Column Updating Method (ICUM),entre outros.
6.2 Método Bayesiano - MCMC
6.2.1 Princípios
Muito do avanço e da popularidade dos métodos bayesianos de inferência estatís-tica, nas duas últimas décadas, se deve ao advento de técnicas de simulação estocástica.Até os anos 90, a intratabilidade analítica de distribuições à posteriori em muitas situaçõespráticas limitava a aplicação da metedologia bayesiana.
Dessa forma, inúmeros métodos numéricos foram propostos, como aproximaçõesde Laplace, aproximações via quadratura Gaussiana, entre outros (ver Taner (1996) eGamerman (1997)). Porém, somente após a inserção de métodos de Monte Carlos viaCadeias de Markov (MCMC) com a conseguinte evolução da capacidade tecnológicacomputacional houve de fato a difusão da metodologia.
A idéia básica do método é a construção de uma Cadeia de Markov que pode sersimulada facilmente, sendo que a distribuição de equilíbrio da mesma seja, para o presentecontexto, exatamente a distribuição a posteriori de interesse.
Desta forma, após um número grande de passos ou iterações da cadeia, o estadoda cadeia é utilizado como uma amostra da distribuição de interesse. A qualidade daaproximação está diretamente relacionada ao número de iterações realizadas, isto é, nãohá dificuldade em se construir a cadeia, mas sim na determinação do número de passosque produzem um erro aceitável.
Para tanto, nas próximas seções serão descritos, resumidamente, os principais as-pectos das Cadeias de Markov e da simulação de Monte Carlo. Após tal introdução,serão discutidos e aplicados ao presente problema dois dos principais métodos MCMC: oamostrador de Gibbs e o algoritmo de Metropolis-Hastings.
38
Ehler (2004) atenta à cautela no uso dos métodos de simulação aqui descritos.Sempre que possível sua resolução, deve-se dar preferência à solução analítica das funçõesde interesse em contraste aos métodos de aproximação numérica. Tomar cuidado com o“Erro Tipo ”, isto é, de termos soluções para um problema errado ou ainda uma soluçãonão muito boa para o problema atual.
6.2.2 Amostrador de Gibbs
Em 1982, Geman e Geman publicaram um artigo na área de processamento deimagens onde propuseram um esquema de amostragem para uma certa distribuição à pos-teriori, contextualmente chamada de “distribuição de Gibbs” (daí o nome do método).Apenas em 1990, Gelfand e Smith compararam tal esquema à outros de simulação es-tocástica, desta forma difundindo-o para a comunidade estatística.
O método baseia-se na simulação de uma Cadeia de Markov onde o núcleo de tran-sição é dado pelas distribuições condicionais conjuntas, dadas em (5.27) por π(θi|θ−i) =π(θi,θ−i)
π(θ−i). Supondo-as conhecidas e ainda que amostras podem ser simuladas de tais dis-
tribuições, o algoritmo consiste nos seguintes passos:1 - Inicializar o contador de iterações da cadeia (k=1) e atribuir um valor inicial
para os parâmetros: θ(0) = (θ(0)0 ,θ
(0)1 , . . . ,θ
(0)p+q)T
2 - Atualizar os parâmetros θ(k) de θ(k−1) , segundo a regra de atualização abaixo:
θ(k)0 ∼ π(θ0|θ
(k−1)1 , . . . ,θ
(k−1)p+q )
θ(k)1 ∼ π(θ1|θ
(k)0 ,θ
(k−1)2 , . . . ,θ
(k−1)p+q )
...
θ(k)p+q ∼ π(θp+q|θ
(k)0 ,θ
(k)1 , . . . ,θ
(k)p+q−1)
3 - Incremente o contador de k para k+1 e volte ao passo 2 até que se tenha aconvergência da cadeia.
Deste modo, quando se tem a convergência da cadeia, o valor resultante θ seráuma amostra da distribuição limite π(θ), isto é, da distribuição à posteriori conjunta deinteresse (ver ([9]), George e Casella, 1992).
Assim, para se obter uma amostra de tamanho n da distribuição, basta replicar oalgoritmo n vezes. Pode-se ainda retirar uma amostra de tamanho “n” após o período de“burn-in”, isto é, o período onde o algoritmo ainda não entrou em fase de convergênciae, dessa forma, as amostras são descartadas. Tal propriedade se dá pelas amostras serem
39
retiradas de uma distribuição estacionária após se assegurar a convergência da mesma (verGamerman,1997).
Como visto, o amostrador de Gibbs faz uso das distribuições condicionais com-pletas à posteriori. Porém, na prática nem sempre é possível a geração de amostras detais distribuições, sendo então utilizados outros algoritmos, por exemplo o de Metropolis-Hastings.
6.2.3 Amostrador de Griddy-Gibbs
Para o modelo Garch, tanto supondo normalidade, quanto supondo outras dis-tribuições, não é possível identificar e consequentemente retirar amostras das distribuiçõescondicionais completas, exceto “forçando” prioris (ver detalhes na seção sobre inferênciabayesiana).
Observa-se, porém que há conhecimento do Kernel π(θi|θ−i), podendo este sercalculado para uma grade ou “grid” de pontos de tamanho m, por exemplo θi1,θi2, . . . ,θim.Esta foi a idéia tida por Ritter e Tanner (1992). O algoritmo é então escrito da seguinteforma:
1 - Inicialize o algoritmo selecionando cuidadosamente um intervalo de pontos deθi: θi1,θi2, . . . ,θim
2 - Avalie cada um dos pontos no Kernel π(θi|θ−i) (suposto conhecido), fazendow j = π(θi j|θ−i) j = 1,2, . . . ,m;
3 - Com w1,w2, . . . ,wm obter uma aproximação da função inversa de distribuiçãode probabilidades acumulada Fπ(θi|θ−i) de π(θi|θ−i);
4 - Para cada j = 1,2, . . . ,m, simular um valor u da distribuição uniforme no inter-valo (0,1) e aplicar a inversa F−1
π (u) para chegar em uma amostra da distribuição condi-cional de interesse.
Como comentário, pode-se calcular a função de distribuição correspondente se-gundo uma regra determinística de integração como a regra do trapézio (Bawens e Lu-brano, 1998), ou ainda a transformação integral de probabilidade (ver Gamerman, 1997)
A escolha da grade deve ser cuidadosa, podendo ser avaliada a partir do histogramada amostra gerada. Caso estas estejam muito concentradas nos extremos do intervalo, ointervalo deve ser expandido. Em contrapartida, caso os valores estejam muito distantesdos extremos, o intervalo deve ser diminuído, pois caso contrário o Griddy Gibbs torna-seineficiente visto que a maioria dos w j seriam zero.
40
6.2.4 Algoritmo de Metropolis-Hastings
Apresentado, em 1953 por Metropolis e generalizado por Hastings em 1970, o al-goritmo de Metropolis-Hastings é um algoritmo de rejeição, isto é, utiliza mecanismos de“correção” para que valores que inicialmente possam ter sido gerados de uma distribuiçãoalvo sejam corrigidos para maior representatividade de tal distribuição. Tal mecanismogarante a convergência da cadeia para a distribuição de equilíbrio 1.
Partindo da propriedade natural de cadeias reversíveis vista anteriormente (5.43),a idéia do algoritmo é a construção de um Kernel de transição K, de tal forma que adistribuição de interesse π(θ) seja a distribuição de equilíbrio da cadeia. Para isso, bastaconsiderar cadeias reversíveis onde o kernel K satisfaça:
π(θ)K(θ,θ′) = π(θ′)K(θ′,θ) , ∀(θ,θ′) (6.11)
Supõe-se ainda disponíveis as distribuições condicionais completas, porém nãosabendo gerar amostras de cada uma delas. Pode-se, entretanto, gerar amostras de umnovo valor de θi a partir de uma distribuição proposta condicional ao valor atual de θi, porexemplo q(θi|θ′i).
Hastings propôs definir uma probabilidade de aceitação α(θi,θ′i) de forma que,
quando combinada com um kernel de transição arbitário, defina uma cadeia reversível.Dentre as diversas formas de se definir tal probabilidade, a que caracteriza os algoritmosde Metropolis-Hastings são dadas por:
α(θi,θ′i) = min
1,
π(θ′)q(θ|θ′)π(θ)q(θ′|θ)
= min
1,
p(ε|θ′)p(θ′)q(θ|θ′)p(ε|θ)p(θ)q(θ′|θ)
(6.12)
onde a segunda igualdade é dada pelo teorema de Bayes (5.20). Observe que nãoé necessário o conhecimento de toda a distribuição de interesse, bastando apenas a parteproporcional.
Assim, o algoritmo é construído da seguinte maneira:1 - Inicializar o contador de iterações da cadeia (k=1) e atribuir um valor inicial
para os parâmetros: θ(0) = (θ(0)0 ,θ
(0)1 , . . . ,θ
(0)p+q)T
2 - Atualizar separadamente cada um dos parâmetros de θ(k−1)i para θ(k) , segundo
a regra de atualização, dada em duas etapas: Para cada i = 0,1,2, . . . , p+q:(i) Gerar um valor θ′i da distribuição auxiliar proposta q(|θ(k−1)
i )
1 Detalhes sobre as propriedades de convergência deste amostrador podem ser encontradas em [23] (Ros-ales, (2004))
41
(ii) Aceitar o valor proposto para θ(k)i com probabilidade dada por:
α = min
1,
π(θ′i|θ(a)−i )q(θ(k−1)
i |θ′i)
π(θ(k−1)i |θ(a)
−i )q(θ′i|θ(k−1)i )
onde,
θ(a)−i = (θ(k)
0 , . . . ,θ(k)i−1,θ
(k−1)i+1 , . . . ,θ
(k−1)p+q )
3 - Incremente o contador de k para k+1 e volte ao passo 2 até que se tenha aconvergência da cadeia.
Sob o ponto de vista prático, apesar da generalidade da escolha de distribuição aposteriori e distribuição proposta aceita pelo algoritmo, a escolha da proposta é crucialpara seu desenvolvimento.
Nesse sentido, há algumas escolhas comuns:A distribuição proposta é simétrica, isto é, q(θ(k−1)
i |θ′i) = q(θ′i|θ(k−1)i ). Nesse caso,
a probabilidade de aceitação é dada por:
α = min
1,
π(θ′i|θ(a)−i )
π(θ(k−1)i |θ(a)
−i )
A distribuição proposta não depende do passo anterior, isto é, q(θ′i|θ(k−1)i ) = q(θ′i).
Nesse caso, a probabilidade de aceitação é obtida por:
α = min
1,
π(θ′i|θ(a)−i )
π(θ(k−1)i |θ(a)
−i )
Particularmente, quando há cadeias independentes e ainda quando a distribuição
proposta é a distribuição a priori do parâmetro em questão, a probabilidade de aceitaçãoé função apenas da verossimilhança, isto é:
α = min
1,
p(ε|θ′i,θ(a)−i )
p(ε|θ(k−1)i ,θ
(a)−i )
Por fim, observa-se que o amostrador de Gibbs é um caso particular do Algoritmo
de Metropolis-Hastings, no sentido de “gerar da condicional completa e aceitar sempreem um algoritmo iterativo é a proria definição do amostrador de Gibbs”.
Comentários: Há ainda outros esquemas de transição de cadeia que combinam oamostrador de Gibbs e o Algoritmo de Metropolis-Hastings. Tais algoritmos são chama-dos “híbridos” e podem ser vistos em Gamerman (1997). Recentemente poderosos algo-
42
ritmos MCMC estão sendo melhor divulgados, como o método MCMC com saltos rever-síveis (ver [17] (Green, 1995)). Para uma descrição detalhada, ver ainda [10] (Robert andCasella, 1995))
7 Escolha e adequação do modelo
7.1 Critérios Utilizadas na escolha
Existem muitas variantes que utilizam conceitos da teoria da informação para es-tabelecer critérios que, em geral, são baseados em estatísticas construídas a partir dosvalores dos resíduos correspondentes a um dado modelo ajustado. Tais critérios sugerema escolha de k para minimizar uma função objetivo que gere o trade-off entre o princípioda parcimônia e a redução na soma dos quadrados.
7.1.1 Akaike Information Criterion - AIC
o AIC ou Akaike Information Criterion (Akaike, 1974)é um índice de qualidadedo ajustamento de modelos estatísticos usando, portanto, para a escolha ou seleção entredois ou mais modelos. Quanto menor seu valor, melhor a adequabilidade do modelo emexplicar a realidade. Em sua forma geral, tem-se:
AIC = 2k−2Log(θ) , (7.1)
onde k é dado pelo número de parâmetros do modelo e Log(θ) é a função de Log-Verossimilhança avaliada eu seu valor máximo, isto é, no estimador de máxima verossim-ilhança.
Como os modelos com maior número de parâmetros têm provavelmente a verossim-ilhança mais elevada, o AIC tende a escolher modelos mais parcimoniosos.
7.1.2 Bayesian Information Criterion - BIC
Assim como o AIC, o BIC (Bayesian Information Criterion), também conhecidocomo Critério de Schwarz, é um índice de qualidade do ajuste usado para a seleção de ummodelo ótimo pertencente a uma família paramétrica da modelos com número diferente de
43
44
parâmetros. Desenvolvido por Schwarz (1978), tem nome devido ao argumento bayesianoadotado pelo autor para justificar o uso do modelo.
É um resultado assintótico derivado da suposição de que os dados pertencem àuma família pertencente à distribuição exponencial. A fórmula geral é dada por:
BIC = k ln(n)−2Log(θ) , (7.2)
onde n é o número de dados observados; k é dado pelo número de parâmetros domodelo e Log(θ) é a função de Log-Verossimilhança avaliada eu seu valor máximo, istoé, no estimador de máxima verossimilhança.
Em comparação com o AIC, o BIC “penaliza” mais a adoção de um parâmetroa mais no modelo, portanto há tendência quanto à escolha de valores menores para onúmero de parâmetros (Perron e Vogelsang, 1992).
7.2 Testes de adequação
Escolhido o modelo segundo os critérios da seção anterior e após a estimação dosparâmetros, a adequação dos dados ao modelo de heteroscedasticidade condicional podeser verificada a partir da série de resíduos resultante.
De acordo com a suposição inicial, de normalidade dos erros condicionados aopassado (εt |ψt−1), espera-se que a série de resíduos padronizados, isto é, a série dos resí-duos na qual foi “retirada” a volatilidade ao se dividir pela raiz da volatilidade:
εt =εt√σ2
t, (7.3)
seja independente temporalmente, com cada observação sendo normalmente dis-tribuída com média zero e variância 1, isto é, tenha distribuição Normal Padrão εt
iid∼N(0,1) 1.
Para verificar tal adequação, alguns testes podem ser aplicados à série ε2t , como: o
teste de existência da dependência temporal significativa da série de resíduos padroniza-dos (Testes de Box-Pierc e de Ljung-Box - muito utilizado nesta pesquisa para validaçãoe identificação dos modelos aplicados), e testes de aleatoriedade aos quadrados dos resí-duos; além de testes para a normalidade dos resíduos padronizados: QQ-plot (teste “vi-sual”) e o teste de Jarque-Bera.
1 Ou que possuam adequação com relação à outra densidade proposta considerada, como a t de student.
45
7.2.1 Teste de Ljung-Box
Baseia-se no correlograma para diagnosticar a independência da série não anal-isando cada um dos “lags”, mas sim baseando-se em todos os possíveis. Para os modelosGARCH, espera-se que, caso os dados estejam bem ajustados, os quadrados dos resíduosse comportem de forma aleatória.
Assim, procede-se um teste que tem como estatística seguindo a distribuição dequi-quadrado. Na literatura ([8]) observa-se o uso de um teste equivalente, porém baseadono rank das observações. Observa-se que, quanto mais pesadas as caudas das distribuições,mais é indicado o uso do teste baseado nos Ranks, situação interessante para casos nor-malmente tratados por modelos GARCH.
7.2.2 Teste de Jarque-Bera
O teste de Jarque-Bera é um teste de normalidade porém com um detalhe im-portante: sua estatística é baseada no terceiro e quarto momentos amostrais, assimetria ecurtose respectivamente. Salienta-se novamente que o teste é aplicado à série dos resíduosnormalizados.
7.3 Diagnósticos de Convergência
Definido o número de iterações do modelo, isto é, o número de passos onde espera-se que a cadeia chegue à convergência, resta saber se, tanto para a distribuição marginal àposteriori dos parâmetros π(θ j|θ− j), quanto para distribuição posteriori conjunta π(θ), asequência entrou de fato em convergência.
Tal problema, que trata da convergência de cadeias de Markov irredutíveis, foiamplamente discutido até então (para referências, ver Paulino et all, 2003). Entretanto,devido ao desconhecimento da distribuição estacionária de interesse, não se sabe e nuncase saberá ao certo o número ótimo de iterações a se adotar para alcançar a convergênciada cadeia, mas algumas estratégias podem ser usadas para um possível diagnóstico.
Dois conceitos são importantes para a compreensão dos métodos de análise deconvergência da cadeia e para posterior cálculo de estimativas da mesma. O período deburn-in ou aquecimento da cadeia: é o período a partir do qual acredita-se que tenhaocorrido a convergência da cadeia, geralmente obtido por inspeção visual. Tais valoresdevem ser descartados para evitar viés das estimativas.
E ainda o thinning: é o “lag” o qual acredita-se haver autocorrelação da série.Dada a construção da cadeia obtida, os dados são naturalmente correlacionados e, dessa
46
forma, pode-se haver muita perda de informação no armazenamento de todos esses val-ores. Dessa forma, escolhe-se um valor k, geralmente por inspeção visual do correlo-grama, e usa-se como amostra final, após a retirada do período de burn-in, os valoresobtidos sistematicamente a cada k iterações da cadeia.
Paulino et all (2003) descreve algumas das principais técnicas usuais de diagnós-tico, que podem e devem ser utilizadas conjuntamente. Abaixo estão resumidas váriastécnicas, sendo que apenas algumas serão utilizadas no presente trabalho:
Representações Gráficas das estimativas das densidades:
Para cada parâmetro do modelo, descarta-se o período de “burn-in”, fixa-se umcerto número de parâmetros da cadeia θ(b+1),θ(b+1), . . . ,θ(b+m) e estima-se graficamentesua densidade; Para cada parâmetro amostrado subsequente θ(b+m+1), recalcular a esti-mativa da densidade e sobrepô-la no gráfico anterior. Fazer isso até que não se consigadistinguir novas curvas.
Monitorar médias Ergódicas:
Toma-se valores distintos iniciais para a cadeia e avalia-se médias ergódicas de al-guns momentos amostrais, percentis, etc. Observado aconvergência das mesmas, amostraruma longa cadeia e fazer inferência.
Método dos quantis de Raftery e Lewis:
Baseado em no problema de se obter intervalos de credibilidade para estimativas,é geralmente utilizado quando se quer encontrar o número de iterações necessárias para aconvergência de uma única cadeia. Para detalhes, ver Paulino et all (2003).
Método de Geweke:
Baseado na aplicação de técnicas de séries temporais para o teste de convergência,o método consiste em dividir a amostra em 2 períodos: um inicial com nI valores e outrofinal, com nF valores. Dada a estacionariedade suposta da cadeia, as médias da parteinicial devem ser a mesma, sendo portanto, basicamente um teste de diferenças de médias.As distribuições utilizadas são de caráter assintótico e podem ser vistas com detalhes emGeweke (1992).
47
Método de Gelman e Rubin:
Utiliza componentes da variância de sequências múltiplas da cadeia, simuladasa partir de uma variedade de pontos iniciais dispersos. Requer maior capacidade com-putacional pois são geradas inúmeras cadeias paralelamente e possui melhor desempenhoquando há conhecimento das distribuições condicionais completas.
Diagnóstico de Heibelberger e Welch:
Baseado na teoria do movimento browniano, testa a hipótese de estacionariedadeda cadeia. Caso a cadeia passe nesse teste, é então procedido o teste de “halfwidth”,baseado no coeficiente de variação da cadeia que “passou” no teste anterior. Caso a cadeiatambém passe nesse teste, retira-se uma amostra e então usa-se a média como estimativada média à posteriori do modelo.
7.4 Adequabilidade do modelo
Além da verificação da adequabilidade dos resíduos para os modelos GARCH,no contexto bayesiano outras medidas de ajuste podem ser usadas para a validação domodelo utilizado.
Sendo ε = (ε1,ε2, . . . ,εm) a amostra anterior obtida, utilizada para a estimação dosparâmetros, a partir de uma nova amostra εnovo = (ε1,ε2, . . . ,εm), denominada amostra de
validação, pode-se calcular a distribuição preditiva desses novos valores, isto é:
p(εnovo|ε) =∫
f (ε|θ) f (θ|εnovo)dθ (7.4)
Com ela, pode-se avaliar o modelo é adequado ao esperado após o ajuste. Calculando-se a média e a variância da distribuição preditiva, chega-se aos resíduos bayesianos apósuma padronização dos dados da amostra de validação. Daí, basta comparar os resíduoscom os valores preditos como normalmente se faz no método clássico. Para mais detalhesde métodos de construção dos resíduos bayesianos, ver Paulino et all (2003).
8 Previsão da Volatilidade
Um dos principais interesses em se ajustar modelos é a possibilidade de previsãopara valores futuros. Escolhido o número de parâmetros do modelo, estimados e testadosquanto a hipótese de nulidade e ainda após a verificação de adequabilidade dos dados aomodelos, o próximo passo é fazer predições (inovações).
No caso de modelos da família GARCH, o interesse se dá principalmente na pre-visão da volatilidade da série. Inicialmente, para o modelo GARCH(p,q) calcula-se aprevisão da volatilidade um passo à frente, fixada no tempo t:
σ2t+1|t = α0 +
q
∑i=1
αiε2t−i+1 +
p
∑j=1
βiσ2t− j+1 (8.1)
Se o intuito é fazer previsões da volatilidade “k” passos à frente, a fórmula geral éencontrada de forma recursiva. Sua expressão é dada por:
σ2t+k|t = E(σ2
t+k|ψt) = α0 +p
∑i=1
αi ·E(ε2t+k−i|ψt)+
q
∑j=1
βi ·E(σ2t+k−i|ψt) (8.2)
porém, deve-se estar atento aos índices pois há mudança quanto ao valor da pre-visão, de acordo com o número de defasagens à frente (k). Para k ≤ i, as esperanças sãodadas por:
E(ε2t+k−i|ψt) = ε
2t+k−i
E(σ2t+k−i|ψt) = σ
2t+k−i ,
pois, nessa condição, o conhecimento de todas as informações até o tempo t
garante o conhecimento dos erros e da volatilidade. Agora, para k > i, as esperançassão dadas por:
48
49
E(ε2t+k−i|ψt) = E(σ2
t+k−i · z2t+k−i|ψt)
= E(σ2t+k−i|ψt) ·E(z2
t+k−i|ψt)
= E(σ2t+k−i|ψt) ·1
= E(σ2t+k−i|ψt) ,
que nada mais são que as previsões de volatilidade k− i passos à frente. Destamaneira, a computação de (8.2) é dada de maneira recursiva.
9 Aplicação
9.1 Um estudo de Simulação
Para fins de comparação e principalmente para validar as ferramentas desenvolvi-das durante esse trabalho, um estudo para o ajuste do modelo GARCH a dados simuladosserá feito nesta seção.
A série a ser simulada será proveniente de um modelo do tipo GARCH(1,1), sendoque os parâmetros do modelo serão dados por α0 = 0.02 , α1 = 0.4 e β1 = 0.3 . Observa-se aqui que com tais parâmetros, α0 +β0 +β1 = 0.72 < 1 , validando a restrição de não-negatividade. Abaixo, a programação no software R para a simulação da série, utilizando-se o pacote “garchFit”:
> specifica <- garchSpec(model = list(omega=c(0.02),alpha=c(0.4),beta=c(0.3)))
> simulado <- garchSim(specifica,n=800)[-(1:100)]
> serie_simulada <- as.ts(simulado,start=1,end=length(simulado),frequency=1)
9.1.1 Análise Descritiva
Foram simulados n = 800 valores da referida série, sendo descartados os 100primeiros, a título de “aquecimento” da amostra. Na figura abaixo encontra-se o gráficoda série simulada, além do histograma, qq-plot e correlogramas da mesma.
Como esperado, são observados pequenos “blocos de variação” na série, carac-terística dos modelos GARCH. Antes da tentativa de ajustar o modelo aos dados simula-dos, será feita uma pequena análise descritiva dos para a verificação da possível adequa-bilidade.
Abaixo encontra-se uma tabela com as principais estatísticas descritivas para asérie simulada. Observa-se que média e mediana possuem valores muito próximos, dandocaráter quase simétrico à série, como pode ser observado no histograma acima. Porém,observando o p-valor para o teste de assimetria, apesar de próximo de 0.05 ao nível deconfiança de 5% rejeita-se a hipótese de simetria da série. Rejeita-se ainda a hipótese
51
(a) Série (b) Histograma, QQ-Plot e Correlogramas
Figura 9.1: Série simulada
de excesso de curtose igual a zero, característica (como visto anteriormente) de sériesmodeladas pelo GARCH, isto é, presença de “caudas pesadas” em relação à distribuiçãonormal.
Tabela 9.1: Estatísticas descritivas para a Série SimuladaEstatísticas Série Simulada Estatísticas Série Simulada
Média 0,01032 Curtose 1,78466 (0)Mediana 0,00666 Jarque-Bera 95,93468 (0)Desvio-Padrão 0,2543 Ljung-Box(lag=1) 0,00263 (0,95912)Mínimo -1,29925 Ljung-Box(lag=5) 10,2437 (0,06862)Máximo 1,22388 Ljung-Box(lag=10) 12,57031 (0,2487)Assimetria 0,3171 (0,040021) Ljung-Box(lag=20) 26,09753 (0,16261)Obs.1: O valor apresentado entre parênteses é o p-valor correspondente.
Obs.2: A notação (0) se refere à p-valor encontrado é menor do que 10−6 .
E, como não poderia ser diferente, dado que o teste de normalidade de Jarque-Bera é baseado na assimetria e na curtose da série, foi fortemente rejeitada a hipótese denormalidade da série ao nível de 5%. Para uma inspeção visual, o gráfico dos quantis(QQ-plot) sugere uma distorção da hipótese de normalidade ao observar caudas pesadasfora dos limites de confiança esperados.
Por fim, verificando a autocorrelação da série via FAC e FACP amostrais, observa-se que não parece haver uma estrutura de decaimento exponencial em nenhum dos gráfi-cos ou qualquer tipo de decaimento brusco para o zero, o que caracterizaria uma estruturaauto-regressiva (AR) ou ainda de médias móveis (MA). Poucos lags parecem conter auto-
52
correlação significativa (ultrapassando os limites de confiança) em ambos gráficos, sendoestes considerados “espúrios”.
Analisando-se a o quadrado dos valores simulados, observa-se um padrão auto-regressivo no mesmo, devido ao decaimento exponencial observado na FAC e o decai-mento brusco da FACP a partir do segundo lag. Tal fato, como dito anteriormente, carac-teriza modelos GARCH.
Por fim, para confirmar a evidência do padrão GARCH é usado o teste de Ljung-Box descrito anteriormente. O teste é feito tanto para os valores simulados, quanto parao quadrado dos mesmo para os Lags 1 a 5, e está descrito na tabela 1 acima. Rejeita-sea independência serial dos valores simulados, assim como para o quadrado dos mesmo,justificando a adoção de um modelo heterocedástico para modelar a evidente volatilidadeda série.
Para o ajuste, supõe-se que o a série tenha média nula e que a mesma seja não-autocorrelacionada (estrutura de um retorno). Caso contrário à essa suposição, pode-seespecificar uma estrutura ARMA à série e depois estimá-la conjuntamente com a estruturaGARCH, via máxima verossimilhança.
Nem sempre a estrutura gaussiana dos erros condicionados ao passado pode seraplicada à séries, pois não conseguem explicar todo o excesso de curtose, invalidando asuposição. Nesse caso será usada também a distribuição t de student, com o parâmetro“graus de liberdade” também sendo estimado.
9.1.2 Estimação
Foram procedidos dois tipos de ajuste: clássico e bayesiano. Como possibilidadede ajustamento havia sete rotinas de cálculos disponíveis inicialmente para a escolha dosparâmetros do modelo. Portanto, dada a dimensão do número de modelos GARCH(p,q)disponíveis para o ajuste e ainda as sete metodologias de estimação dos parâmetros alémde todas as séries analisadas, optou-se pelo enfoque na comparação das rotinas.
Isto é, foi feita a escolha inicial das dimensões (p,q) de um modelo e em seguidaforam comparadas as metodologias de cálculo de acordo com a estimação dos parâmetrosdo modelo, além do cálculo de erro padrões e intervalos de confiança/credibilidade paraos mesmos. Por fim, após a seleção de uma metodologia ótima e a subseqüente escolhados parâmetros do modelo para cada uma das séries, avaliou-se de maneira prática talescolha de acordo com a previsão da volatilidade série para o mês de dezembro de 2008,comparando com os valores iniciais.
Para cada um dos parâmetros, segundo o modelo Bayesiano são escolhidos doisconjuntos de distribuições à priori:
53
1. θ j ∼U(0,1) , ∀ j = 0,1,2
2. Sendo Yj∼N(0,h j) e usando a reparametrização de Yj = log θ j1−θ j
para θ j = 11+exp−Y j
,mostra-se facilmente que x ∈ (0,1), sendo portanto uma variável aleatória no inter-valo (0,1), adequada ao modelo.
Para cada uma das prioris foram utilizados os Algoritmos de Metropolis-Hastingse Griddy-Gibbs para a construção de uma amostra da distribuição à posteriori do modelo.Com a amostra obtida, é analisada a convergência das mesmas segundo alguns critériosdescritos anteriormente.
Uma das maiores críticas à modelagem bayesiana é com relação à influência dadistribuição à priori no resultado final da modelagem. Para o presente caso,devido ao usode prioris não informativas, fazendo uso de uma análise de sensitividade, os resultadosmostram a influência quase insignificante das prioris para o modelo em questão.
Os métodos de estimação dos parâmetros utilizados são:
• M1: Método clássico via algoritmo de NR, supondo inovações Normais;
• M2: Método clássico via algoritmo de NR, supondo inovações T de Student;
• M3: Método Bayesiano, via algoritmo de MH supondo inovações Normais comprioris não informativas U(0,1);
• M4: Método Bayesiano, via algoritmo de MH supondo inovações Normais comprioris dadas no item 2 acima;
• M5: Método Bayesiano, via algoritmo de MH supondo inovações T de Student comprioris não informativas U(0,1) e os graus de liberdade com priori cauchy padrão;
• M6: Método Bayesiano, via amostrador de GG supondo inovações Normais comprioris não informativas U(0,1);
• M7: Método Bayesiano, via amostrador de GG supondo inovações T de Studentcom prioris não informativas U(0,1) e os graus de liberdade com priori cauchypadrão;
Para a seleção dos modelos que melhor explicassem o comportamento de cadauma das séries de retornos analisadas foram utilizados os critérios de seleção AIC e BIC,segundo o método M1. Portanto, após a escolha de (p,q), os parâmetros do modelo foramestimados segundo os métodos descritos acima.
54
Para a expressão da Verossimilhança segundo a distribuição t de Student, verapêndice F.
Para o Método frequentista foi utilizado o método de Newton-Raphson para ocálculo das estimativas, tanto supondo a normalidade, quanto o compormento dos erroscondicionados ao passado segundo a distribuição t de student.
A tabela (9.2) abaixo compara os valores para cada um dos critérios AIC e BIC,para os modelos frequentistas, supondo normalidade dos erros condicionados, sendo queos menores valores estão em negrito. Novamente reforça-se aqui que tal escolha, mesmonão sendo a ideal, é viável dada a dimensão de possibilidades de ajustes escolhidas parase fazer nesse trabalho.
Tabela 9.2: Série Simulada - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 -0,0204(30) -0,0415(1) -0,0391(3) -0,038(6) -0,0359(10) -0,0345(13)
-0,0009(6) -0,0155(1) -0,0066(3) 0,001(7) 0,0096(11) 0,0175(16)
2 -0,0401(2) -0,039(4) -0,0363(7,5) -0,0351(11) -0,0331(16) -0,0317(19)-0,0141(2) -0,0065(4) 0,0027(8) 0,0104(12) 0,0189(17) 0,0268(21)
3 -0,0387(5) -0,0363(7,5) -0,0335(15) -0,0323(18) -0,0302(21) -0,0288(24)-0,0061(5) 0,0028(9) 0,012(15) 0,0197(19) 0,0283(22) 0,0362(25)
4 -0,0362(9) -0,0341(14) -0,0312(20) -0,0294(23) -0,0274(25) -0,026(27)0,0029(10) 0,0115(14) 0,0208(20) 0,0291(24) 0,0376(26) 0,0455(28)
5 -0,0348(12) -0,0327(17) -0,0298(22) -0,0271(26) -0,0249(28) -0,0233(29)0,0107(13) 0,0193(18) 0,0287(23) 0,0379(27) 0,0466(29) 0,0548(30)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
Analisando-se a tabela (9.2), os dois crtérios utilizados apontam para a escolha domodelo GARCH(1,1) como o mais adequado para o ajuste da série simulado, como erade se esperar.
Um ponto importante a se destacar é a proximidade do valor dos crtérios paraalguns parâmetros diferentes, como por exmplo o (1,1) e (2,1) para o BIC. Algumas al-ternativas foram propostas na literatura para lidar com casos assim (ver Gamerman eLopes, 1997), como transformações nos índices visando melhor percepção de diferenças,ou ainda no caso bayesiano, como será discutidas na conclusão. Para esse problema,escolhemos portanto o modelo GARCH(1,1) para o ajuste.
A tabela (9.5) mostrada bem abaixo condensa as informações mais importantesresultante dos ajustes para cada um dos modelos M1-M7. Nela encontram-se as médias,erros padrões, intervalos de confiança/credibilidade e ainda o p-valor dos testes de nuli-dade de cada um dos parâmetros de cada um dos modelos analisados.
55
Para os modelos gerados segundo a metodologia frequentista (M1 e M2), escontram-se abaixo os gráficos mais importantes para a análise: Tanto os resíduos, quanto os mesmopadronizados para os dois modelos aparentam aleatoriedade, isto é, nenhum sinal de com-portamento sistemático. Já os correlogramas tanto dos resíduos padronizados, quanto doquadrado dos mesmo, não apresenta quaisquer sinais de autocorrelação das série, comodesejado.
(a) M1 (b) M2
Figura 9.2: Série simulada - Gráficos dos resíduos dos modelos
Por fim, o qq-plot apresenta um diagnóstico visual favorável à suposição de nor-malidade dos resíduos, isto é, para o modelo M1. Já para o qq-plot dos resíduos segundoo modelo M2, isto é, supondo erros distribuídos de acordo com a distribuição t, pare-cem não se ajustar bem à linha central do gráfico, dando indícios de que não houve boamodelagem.
Finalmente, para os dois modelos frequentistas, pode-se dizer que houve bomajuste dos dados, pois em ambos os resíduos padronizados se comportam bem no sentidoda adequação explicada anteriormente.
Como dito, para os modelos Bayesianos M3-M7 também estão explicitados natabela (9.5) mais abaixo os resultados da estimação dos parâmetros. Porém, no casobayesiano é importante a análise dos dignósticos de convergência de cada uma das cadeiasutilizadas para a estimação dos parâmetros à posteriori.
Após inspeções iniciais de convergência dos modelos, para cada um dos algorit-mos foi gerada uma cadeia de tamanho n = 18000 e analisou-se a convergência de taiscadeias segundo os critérios de diagnósticos de Geweke, além dos correlogramas para a
56
escolha do período de “thinning” (no caso 100) e finalmente o teste diagnóstico de Hei-delberg e Welch para a escolha do período de “burn-in” (300 amostras iniciais retiradas).
Feito isso, a análise da convergência da amostra final gerada é procedida a seguir.Abaixo a série para cada um dos parâmetros e ainda a densidade a posteriori aproximadados mesmos.
Na tabela (9.3) abaixo, que contém os testes de Heidelberg e Welch para cadaparâmetro em cada modelo, observa-se que apenas para o modelo M3 o parâmetro α0
encontrou problemas quanto à convergência, falhando já inicialmente no teste de estacio-nariedade da cadeia. Já para o parâmetro “graus de liberdade” denotado por t apesar depassar no teste de estacionariedade, as cadeias falharam no teste de halfwidth, indicandoque as médias à posteiori de tais parâmetros para os modelos M5 e M7 não podem/devemser usadas como estimativa.
Tabela 9.3: Série simulada - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Falhou Passou Passou Passou Passou
NA Passou Passou Passou Passou
α1 Passou Passou Passou Passou PassouPassou Passou Passou Passou Passou
β1 Falhou Passou Passou Passou PassouNA Passou Passou Passou Passou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
Agora, uma tabela equivalente à de cima, porém para o teste de Geweke, é dadaabaixo em (9.4). Nela, como, exceto para o modelo M6, todos os parâmetros apresentamscore menor em módulo que 1,97, isto é, abaixo dos limites de confiança, conclui-se quetodos os parâmetros parecem convergir.
Tabela 9.4: Série simulada - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 0,753460361 -0,941724719 0,388990836 -1,738428241 -1,835987061α1 0,057825527 0,179069885 0,50616884 -1,587164914 -1,385668646β1 -0,821672295 0,730598365 -0,197218795 2,115547835 1,390318888t - - 0,663226317 - -0,265042796
Fonte: Dados da Estação Automática de Brasília - INMET.
Como comentários finais, devido ao custo computacional não foram geradas maiscadeias para então proceder o teste de Gelman e Rubin. Observa-se que é sempre im-portante a realização de mais de um teste de convergência, como no caso. Deste modo,apesar do teste de Heidelberg e Welch indicar convergência não clara para alguns parâmet-ros, pode-se considerar o teste de Geweke como indício de convergência. Por fim, como
57
é garantida a convergência devido à construção natural dos algoritmos, uma saída para aaceitação da convergência é o aumento do número de iterações da cadeia.
Finalmente, na tabela (9.5) abaixo, odos os testes rejeitam a hipótese de nulidade,restando-nos portanto aceitar as estimativas. Observa-se que os valores dos parâmetrossão um pouco diferentes dos valores originais,todos os intervalos contém os verdadeirosvalores dos parâmetros.
58
Tabe
la9.
5:E
stim
ativ
asdo
spa
râm
etro
spa
raa
séri
esi
mul
ada
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0=
0,02
Est
imat
iva
0,02
0141
0,02
1713
0,02
190,
0219
0,03
180,
0211
0,02
15E
rro
Padr
ão0,
0037
570,
0047
740,
0042
0,00
340,
0137
0,00
390,
004
IC(0
,95)
(0,0
1276
4;0
,027
518)
(0,0
1234
;0,0
3108
5)(0
,013
8;0
,029
8)(0
,015
;0,0
285)
(0,0
154
;0,0
548)
(0,0
137
;0,0
291)
(0,0
142
;0,0
302)
P-va
lor
00,
0000
050
00
00
α1
=0,
4E
stim
ativ
a0,
4618
440,
4567
330,
456
0,46
620,
4448
0,45
630,
4531
Err
oPa
drão
0,08
1349
0,09
6042
0,07
430,
0785
0,14
470,
0759
0,08
57IC
(0,9
5)(0
,302
125
;0,6
2156
3)(0
,268
167
;0,6
4529
9)(0
,321
2;0
,599
1)(0
,340
2;0
,619
9)(0
,160
1;0
,725
2)(0
,323
;0,5
88)
(0,2
979
;0,6
294)
P-va
lor
00,
0000
020
00
00
β1
=0,
3E
stim
ativ
a0,
2620
930,
2804
160,
2529
0,24
110,
1529
0,26
040,
2596
Err
oPa
drão
0,08
4367
0,10
0255
0,08
710,
0704
0,13
150,
0828
0,08
5IC
(0,9
5)(0
,096
449
;0,4
2773
7)(0
,083
578
;0,4
7725
4)(0
,107
4;0
,435
2)(0
,119
6;0
,383
)(0
,004
3;0
,498
8)(0
,114
4;0
,455
3)(0
,104
;0,4
411)
P-va
lor
0,00
1893
0,00
5158
00
00
0
tE
stim
ativ
a-
10-
-6,
5076
-14
9,92
59E
rro
Padr
ão-
2,14
4652
--
23,1
873
-96
,549
7IC
(0,9
5)-
(5,7
8924
7;1
4,21
0753
)-
-(0
,036
4;7
1,76
06)
-(3
4,64
3;4
13,2
333)
P-va
lor
-0,
0000
03-
-0,
0003
-0
Font
e:D
ados
daE
staç
ãoA
utom
átic
ade
Bra
sília
-IN
ME
T.
59
Em uma análise mais minuciosa, observa-se que os intervalos de confiança/credibilidadede menor amplitude para cada um dos parâmetros é α0 (M4), α1 (M6), β1 (M4) e t (M2).Exceto para o último parâmetro, todos possuem melhor estimativa de acordo com o mo-delo Bayesiano, atestando a qualidade do método. Observa-se maior robustez dos méto-dos baseados no algoritmo de Metropolis-Hastings quando comparados aos advindos doalgoritmo de Greeddy-Gibbs, porém para o modelo M5, isto é, via MH com inovações t,o algoritmo talvez não convirja muito bem e produza grande variabilidade quando com-parado com outros.
Como comentário final, observa-se a dificuldade dos algoritmos em produzir boasestimativas para o parâmetro t, no sentido da grande variabilidade dos mesmos, pro-duzindo intervalos de credibilidade muito grande. Tal fato pode ser explicado pelo usode uma distribuição à priori que, apesar de não informativa, possibilita a geração de val-ores em um intervalo muito grande (−∞,∞), dificultando a convergência da cadeia para adistribuição limite, como observado nos diagnósticos de convergência anteriores.
9.2 Aplicação a dados reais
Estudadas as técnicas, será procedido o ajuste de modelos da família GARCH paraalgumas séries meteorológicas.
O roteiro da apresentação desta seção será o mesmo utilizado para o estudo desimulação feito na seção anterior, isto é, com uma pequena descrição dos dados, umaanálise descritiva e, por fim, a comparação entre os dois métodos de estimação. Nova-mente, não serão detalhados os aspectos da convergência de todas as amostras simuladaspara cada um dos métodos MCMC utilizados, para cada uma das variáveis, devido aonúmero extenso de cadeias simuladas, mas muitos gráficos importantes para a análise ecompreensão das análises estão disponíveis nos Apêndices deste trabalho.
Foram analisadas séries históricas de 6 variáveis meteorológicas: Temperaturamáxima (C), Temperatura mínima (C), Umidade Relativa do ar(%), Velocidade doVento (média), Insolação (horas de sol), Evaporação(mm) . São dados diários para aestação automática do instituto de meteorologia (INMET), que compreendem o períodode 01/01/2007 a 30/11/2008, isto é, séries com 700 valores utilizados para o ajuste domodelo. O mês de dezembro de 2008 foi utilizado para a validação do modelo, no sentidode se fazer previsões e comparar com os dados reais.
Outras variáveis foram analisadas previamente a fim de serem colocadas nesterelatório, porém alguns problemas ocorreram. Para a série de precipitação, mesmo usandoalgumas técnicas para tentar contornar a situação, o número elevado de zeros complicavao ajuste de quaisquer retornos aos dados, devido à divisão por zero ser impossível. A
60
temperatura média também foi descartada, pois já é usualmente modelada na literatura,sendo que o intuito inicial do trabalho era trabalhar com dados “extremos”.
As variáveis “Radiação de onda Longa” (ROL) e “Radiação de onda curta” (ROC),apesar de serem previamente bem modeladas segundo o GARCH, suas análises não foramcolocadas aqui em razão da indisponibilidade dos dados para períodos atuais, isto é, asestações do INMET não mais coletam tais dados desde o ano 2006. Deste modo, ape-sar de bem modeladas, não há interesse prático para o instituto em se fazer previsões“defasadas”, mas sim previsões atuais, diárias.
Outras variáveis ainda podem/poderiam ser analisadas, como a pressão atmos-férica, porém as estações automáticas não disponibilizam outras informações além dasutilizadas. A princípio fora pensado em utilizar dados das estações convencionais do IN-MET, porém tais dados são sigilosos e requerem autorização para sua utilização. Mesmoautorizado, são disponibilizados apenas dados para, no máximo, 6 meses, o que não émuito bom para a modelagem via GARCH, que requer séries bem longas para devidaatualização das volatilidades.
Por fim, os aspectos computacionais para o cálculo de todas as estimativas serãodiscutidos na próxima seção.
9.2.1 Análise descritiva dos retornos
Abaixo, na figura (9.3) encontram-se as séries temporais analisadas com as re-spectivas séries dos retornos das mesmas. Como dito anteriormente, é comum a análisedos retornos das variáveis em séries financeiras, pois os dados em si não apresentam a es-trutura de “blocos de volatilidade” característica do GARCH. E esta simples visualizaçãográfica exemplifica tal fato, pois em nenhuma das séries originais há aparente comporta-mento de períodos de alta volatilidade seguidos de períodos de alta volatilidade, enquantoperíodos de baixa seguidos de períodos de baixa. Mas as séries dos retornos parecemapresentar tal aspecto de blocos de variação.
Alguns são bem aparentes como a série dos retornos de insolação, onde a vari-abilidade é grande no início, há decaimento brusco, com os dados bem “comportados”,enquanto que ao final do ano de 2007 os dados voltam a possuir grande variabilidade,se estendendo até o meio do ano, quando decai novamente. Para a série dos retornos deevaporação também se percebe um bloco bem aparente entre o ano de 2007 e o ano de2008, porém analisando a série da variável o comportamento errático parece ter grandevariabilidade em todos os aspectos. Para temperatura mínima também é aparente a grandevariabilidade no início e no fim dos anos, o contrário para a temperatura máxima.
61
(a) Evaporação (b) Insolação (c) Temp. Máx.
(d) Temp. Min. (e) Um. Rel. (f) Vel. Vento
Figura 9.3: Séries originais e Retornos
Destaca-se o aparente comportamento cíclico de quase todas as variáveis anal-isadas para a estação automática de Brasília. Uma possível abordagem para a modelagemdessas séries seria considerar uma componente cíclica para a série original, transformaressa série modificada para os retornos e, por fim, analisar segundo a perspectiva de het-erocedasticidade, isto é, aplicar o modelo GARCH. Como já dito, o intuito do trabalhoatual é de analisar apenas os retornos das variáveis sendo, portanto, não observado nen-hum comportamento cíclico nas séries de retornos. Melhor discussão será feita ao finaldo trabalho.
Novamente, abaixo encontra-se a tabela (9.6) com as principais estatísticas des-critivas de cada uma das séries de retornos a ser analisada: média, mediana, desvio padrão,valores máximo e mínimo das séries. Há ainda os coeficientes de assimetria e de curtosecom respectivos p-valores dos testes de adequabilidade e por fim os testes de Jarque-Berapara a normalidade e de Ljung-Box para a independência serial.
Para as séries de temperatura mínima e máxima, evaporação e velocidade dosventos, os valores de média e mediana estão muitos próximos, dando indícios de havercerta simetria em cada uma das séries. Ao nível de 5 % de confiança, confirma-se nãohaver indícios de rejeição de assimetria das séries pela simples inspeção dos p-valores(maiores que 0,05) associados aos testes. Para as outras variáveis, apesar de proximidadeentre média e mediana, os testes indicam a rejeição da hipótese de assimetria das séries.
62
Tabela 9.6: Estatísticas descritivas das séries dos RetornosEstatísticas Variáveis
Temp. Max. Temp. Min. Um. Rel. Vel. Vento Insol. Evap.Média 0,0002 -0,0001 -0,0002 0 0,001 -0,0002Mediana 0,0031 0 -0,0055 0 0 -0,0247Desvio-Padrão 0,0667 0,0797 0,185 0,8283 0,8838 0,5669Mínimo -0,254 -0,3649 -0,5995 -2,6626 -4,6347 -2,8717Máximo 0,2709 0,2603 0,8766 2,6856 4,5951 2,8134Assimetria -0,2272 (0,1067) -0,2026 (0,1492) 0,5336 (0,0003) -0,0228 (0,8699) 0,4083 (0,0047) 0,1817 (0,195)Curtose 5,384 (0) 4,2359 (0) 6,0442 (0) 3,8434 (0,0006) 12,0417 (0) 8,5853 (0)Jarque-Bera 171,5449 (0) 49,2724 (0) 303,0814 (0) 20,7791 (0) 2400,4436 (0) 912,4337 (0)Ljung-Box(lag=1) 32,5086 (0) 52,2485 (0) 13,4058 (0,0003) 77,5709 (0) 125,7671 (0) 89,7195 (0)Ljung-Box(lag=5) 55,9462 (0) 70,5027 (0) 50,2539 (0) 98,8705 (0) 149,017 (0) 99,0024 (0)Ljung-Box(lag=10) 65,0326 (0) 79,1964 (0) 62,347 (0) 111,2179 (0) 153,8936 (0) 99,5436 (0)Ljung-Box(lag=20) 72,367 (0) 109,9602 (0) 80,0921 (0) 131,6383 (0) 182,401 (0) 122,3535 (0)Fonte: Dados da Estação Automática de Brasília - INMET.Obs.1: O valor apresentado entre parênteses é o p-valor correspondente.Obs.2: A notação (0) se refere à p-valor encontrado ser menor do que 10−6 .
Como a simetria é característica de retornos, o comportamento assimétrico das séries deumidade relativa e de insolação pode acarretar problemas futuros quanto aos ajustes.
Quanto à curtose, isto é, o grau de achatamento em relação à distribuição nor-mal, todas as séries apresentam característica leptocúrtica (coeficientes maiores que 3),isto é, caudas mais “pesadas”, ou ainda maior achatamento em relação à normalidade.Como descrito na fundamentação teórica, é de interesse séries com caudas pesadas para amodelagem, pois caracterizam modelos GARCH.
Para melhor inspeção visual das séries dos retornos, na figura (9.4) abaixo se en-contram os histogramas para cada uma delas. É visível o caráter leptocúrtico dos dados, eainda a assimetria. Observa-se ainda a grande concentração dos valores em torno de zero,característica da transformação dos dados para os retornos.
(a) Histograma (b) QQ-Plot
Figura 9.4: Retornos das seis variáveis analisadas
63
Quanto à normalidade, a tabela apresentou apenas o teste de normalidade deJarque-Bera, descrito em seções anteriores. Como se sabe, tal teste é baseado nos co-eficientes de assimetria e de curtose da série sendo que, nesse caso, é natural a rejeiçãode hipótese de ajustamento com relação à distribuição dos dados, em face dos resultadosanteriores.
Na figura acima se encontram também os gráficos “QQ-plot” para a inspeção vi-sual da normalidade dos dados. Neles é possível melhor visualização do grau de achata-mento elevado das séries, caracterizando caudas pesadas e o conseqüente afastamento dacurva em relação ao esperado caso os dados pudessem ser ajustados à distribuição Nor-mal.
Para o final desta análise, é avaliado o grau de correlação serial, isto é, se a sérietemporal dos retornos de cada uma das variáveis apresenta alguma estrutura de correlaçãoao longo dela. Pelos testes de Ljung-Box, rejeita-se a hipótese de independência serialpara os lags 1, 5, 10 e 20. Apesar de não se esperar tais resultados no sentido da mod-elagem, os mesmo são totalmente condizentes no sentido das séries dos retornos seremnaturalmente construídas de acordo com uma estrutura “autoregressiva”, apesar de nãohaver uma estrutura pronta para a mesma.
Para uma análise visual, na figura (9.5) se encontram os correlogramas (FAC ver-sus lags) e os correlogramas parciais (FACP versus lags) de cada uma das séries dosretornos.
Figura 9.5: Correlogramas dos retornos das seis variáveis
64
Observa-se que os correlogramas de todas as variáveis apresentam decaimentobrusco no primeiro lag, sendo este ainda pouco efetivo, apesar de estar fora do intervalo deconfiança. Como dito anteriormente, é esperada certa correlação para uma defasagem emtodos as séries devido à construção natural da transformação via retornos, sendo portantoesperado tal comportamento. Observa-se ainda que, para a maioria das variáveis outroslags parecem ter significância, isto é, estão fora dos limites de confiança, porém podemser considerados espúrios no sentido de não serem significativos a ponto de apontar umaestrutura autoregresseiva ou de médias móveis para as séries.
Quanto aos correlogramas parciais, exceto para a variável umidade relativa, oscorrelogramas parciais parecem se comportar de acordo com um decaimento exponencialpara zero, talvez indicando uma estrutura de médias móveis para tais séries de retornos.Nesse caso poder-se ia tentar analisar as séries segundo outras abordagens como via mod-elos ARMA e, em seguida, tentar acolplar uma estrutura GARCH aos resíduos. Porém,nesse trabalho o intuito é apenas de modelar segundo modelos da família GARCH.
Agora, quando observada a série dos quadrados dos retornos na figura (9.6), não éclara a estrutura de autocorrelação necessária para se considerar um modelo heterocedás-tico, porém algumas ressalvas podem ser feitas.
Figura 9.6: Correlogramas dos quadrados dos retornos das seis variáveis
Inicialmente todos os correlogramas aparentam mesmo decaimento brusco parazero a partir do primeiro lag, assim como fora observado anteriormente, porém agora ovalor da FAC para este lag é mais significativa, principalmente para as séries de temper-atura mínima, insolação e evaporação. Se considerarmos que tal decaimento não é brusco
65
para zero, mas sim de forma exponencial (desconsiderando lags espúrios) e ainda que oscorrelograma parciais aparentem decaimento brusco para zero (lags não significativos),pode-se admitir como razoável a hipótese de autocorrelação da variância, ao menos paraas séries citadas no parágrafo anterior.
Dessa forma, apesar de não haver contundência na análise das séries dos quadra-dos dos retornos da forma como nos dados simulados na seção anterior, isto é, na rápidae clara identificação de uma estrutura autoregressiva para essas séries, serão procedidosos ajustes com justificativa baseada principalmente nos testes de independência de Ljung-Box, pois neles rejeita-se a independência serial das séries, assim como para as séries dosquadrados, permitindo adoção de um modelo heterocedástico para modelar a volatilidade.
9.2.2 Estimação
Novamente supõe-se que as médias de todas as séries sejam nulas, razoável por setratarem de retornos, e ainda que as mesmas não sejam correlacionadas, razoável devidoàs correlações “espúrias” vistas anteriormente.
Foram estimados os parâmetros de cada um dos 7 modelos descritos anteriormenteM1-M7, sendo que análise idêntica será feita para cada variável climatológica. Ao final,foi considerado um modelo ótimo para a análise de cada variável para conseguinte esti-mação da volatilidade da série no período remanescente do ano de 2008 (mês de dezem-bro).
Como no caso da estimação, os parâmetros foram selecionados de acordo comos critérios AIC e BIC para o modelo frequentista considerando erros normais, sendoomitidas aqui as tabelas (ver Apêndice C). Os modelos ótimos baseando-se no princípioda parcimônia quando houvesse discordância entre o valor do BIC e do AIC, para cadavariável, são:
• Evaporação = GARCH(1,0) = ARCH(1)
• Insolação = GARCH(2,5)
• Temperatura Máxima = GARCH(1,3)
• Temperatura Mínima = GARCH(1,1)
• Umidade Relativa = GARCH(1,2)
• Velocidade do Vento = GARCH(1,0)
Novamente serão construídas tabelas que condensam toda a informação a respeitodas estimativas de cada um dos parâmetros, a ser mostrada após a análise da adequação
66
do modelos no caso frequentista (M1 e M2), e da análise de convergência e adequação nocaso bayesiano (M3-M7).
Antes de iniciar as análises um ponto importante quanto à construção das cadeias:para cada modelo bayesiano foram geradas 3 amostras com 4000 iterações cada, iniciandode pontos diferentes: do default (valores bem baixos); usando as estimativas do métodoclássico e finalmente de um ponto aleatório. Foi então procedidos testes de Gelman eRubin para a detecção de dependência de valroes iniciais e, no geral, devido ao númeropequeno de iterações em cada amostra, mostrou-se grande dependência, sendo necessárioportanto um número maior de iterações para a convergência da série.
Depois desse estudo, apesar do diagnóstico de Raftery e Lewis indicar amostrasmuito grandes (variando de 10.000 a 100.000), devido ao custo computacional não foipossível gerar amostras maiores que 10.000, sendo padronizado, após análise dos correl-ogramas e do teste de Raftery, um período de burn-in de 300, com thinning variando de 80a 100. De posse dessa amostra final, foi testada a convergência dos parâmetros segundo odiagnóstico de Geweke e ainda de Heidelberg, com os resultados mostrados aqui.
Evaporação
Para a análise frequentista segundo os modelos M1 e M2, é rejeitada a normali-dade dos resíduos (respectivamente o ajuste quanto à distribuição t), como pode ser bemvisualizado na figura (9.7). Os resíduos não aparentam comportamento autoregressivo,sendo encontradas correlações espúrias por volta do lag 17 para o quadrado dos mesmosem ambos modelos. Pode-se considerar nesse caso um ajuste melhor segundo o modeloM2, dado que o qq-plot parece ser menos violado visualmente.
Para os modelos Bayesianos, como pode ser visto no Apêndice D, apesar doprimeiro parâmetro não aparentar problemas quanto à sua convergência, todos os out-ros tiveram problemas no diagnóstico de Heidelberg. Em todos os casos a adoção de talmodelo parece ser bem complicada para o atual número de iterações. Conclui-se que ex-ceto para o modelo M4, apesar de aceita a estacionariedade, os parâmetros α1 não sãoadequados para a estimativa pontual. O modelo que mais se aproximou da convergênciacompleta foi o M5, porém com problema no parâmetro t.
Segundo o teste de Geweke, o parâmetro α0 do modelo M4 e o α1 dos modelosM4 e M6 não convergem, sendo portanto necessário um número maior para a devidaconvergência das cadeias.
Por fim, a tabela (9.7) resume as estimativas para a série. Segundo ela, para osmodelos Bayesianos, o mais adequado parece ser o modelo M3 que obteve bons diagnós-ticos de convergência da amostra simulada da posteriori e teve intervalos de amplitude
67
(a) M1 (b) M2
Figura 9.7: Evaporação - Gráficos dos resíduos dos modelos
bem menor para seus parâmetros estimados. Observa-se aqui a grande diferença entre osvalores estimados para cada um dos modelos. Enquanto modelos frequentistas atribuiemvalores grandes ao α1, exceto para o modelo M7 os modelos bayesianos estimam comovalores pequenos da ordem de 0.05.
68
Tabe
la9.
7:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Eva
pora
ção
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0E
stim
ativ
a0,
1667
370,
1728
820,
1239
0,09
110,
0725
0,45
610,
4006
Err
oPa
drão
0,01
2348
0,03
9289
0,01
0,00
940,
0161
0,25
620,
2844
IC(0
,95)
(0,1
4249
2;0
,190
981)
(0,0
9574
3;0
,250
021)
(0,1
033
;0,1
431)
(0,0
723
;0,1
004)
(0,0
456
;0,1
024)
(0,0
129
;0,8
906)
(0,0
1;0
,904
4)P-
valo
r0
0,00
0011
00
00
0
α1
Est
imat
iva
0,52
2671
0,95
3502
0,03
710,
1386
0,14
680,
0698
0,26
6E
rro
Padr
ão0,
0845
970,
2777
140,
0222
0,03
330,
0609
0,16
220,
2933
IC(0
,95)
(0,3
5657
5;0
,688
767)
(0,4
0824
6;1
,498
757)
(0,0
037
;0,0
858)
(0,1
;0,2
068)
(0,0
376
;0,2
745)
(3e-
04;0
,717
9)(0
,001
2;0
,888
5)P-
valo
r0
0,00
0596
00
00
0
tE
stim
ativ
a-
2,96
8376
--
82,4
207
-11
,775
5E
rro
Padr
ão-
0,40
5008
--
1047
,265
2-
15,3
393
IC(0
,95)
-(2
,173
194
;3,7
6355
8)-
-(0
,065
3;2
9,11
96)
-(0
,307
6;5
8,81
44)
P-va
lor
-0
--
0,29
65-
0Fo
nte:
Dad
osda
Est
ação
Aut
omát
ica
deB
rasí
lia-I
NM
ET.
69
Portanto, escolhidos os modelos M3 e M1, será procedida a estimação da volatili-dade para os mesmos e depois a previsão para o mês de dezembro.
Insolação
Para a análise frequentista segundo os modelos M1 e M2, é rejeitada a normali-dade dos resíduos (respectivamente o ajuste quanto à distribuição t), como pode ser bemvisualizado na figura (9.7). Os resíduos não aparentam comportamento autoregressivo,sendo encontradas apenas correlações espúrias. Nenhum dos resíduos dos modelos M1 eM2 parecem se adequar bem, sendo então considerada estimativas ruins.
(a) M1 (b) M2
Figura 9.8: Insolação - Gráficos dos resíduos dos modelos
Para os modelos Bayesianos, como pode ser visto no Apêndice D, todos os mod-elos tiveram problemas no diagnóstico de Heidelberg. Em todos os casos a adoção de talmodelo parece ser bem complicada para o atual número de iterações. O ideal seria umnúmero maior de iterações, ou ainda a tentativa da adoção de um modelo mais parcimo-nioso. Mesmo com a adoção de prioris que, apesar de não-informativas, possuem umavariabilidade menor que a U(0,1), o modelo ainda sim não foi capaz de fazer boas estima-tivas. O modelo que talvez tenha melhor se aproximado de uma convergência tenha sidoo M3.
Segundo o teste de Geweke (ver Apêndice E), o modelo M6 parece ser o únicoque aparenta convergência, enquanto que todos os outros apresentam problem em algumacadeia. Nota-se a tamanha divergência dos diagnósticos pois, para o modelo M1 nãohouvera problema nenhum nos parâmetros α0 e α1, ao contrário daqui.
70
Por fim, a tabela (9.7) resume as estimativas para a série. Segundo ela, não parecehaver, ao menos visualmente, um modelo dito ao menos razoável para a explicação dosretornos. Enquanto os modelos clássicos não rejeitam a nulidade de alguns parâmetros,os modelos bayesianos tem inúmeros problemas quanto a convergência das cadeias e,consequentemente, atrapalha a validade da estimação.
Nesse caso, não será adotado modelo algum para as estimativas, porque talvez defato não haja adequação dos dados ao modelo.
71
Tabe
la9.
8:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Inso
laçã
ose
gund
oos
sete
mod
elos
Parâ
met
ros
Est
atís
ticas
Mod
elos
M1
M2
M3
M4
M5
M6
M7
α0
Est
imat
iva
0,00
3234
0,00
193
0,05
740,
022
0,03
730,
3406
0,02
74E
rro
Padr
ão0,
0010
570,
0008
520,
0183
0,02
580,
0095
0,27
580,
0349
IC(0
,95)
(0,0
0116
;0,0
0530
9)(0
,000
256
;0,0
0360
3)(0
,003
5;0
,078
4)(0
,001
4;0
,070
8)(0
,023
8;0
,060
2)(0
,01
;0,8
475)
(0,0
019
;0,1
267)
P-va
lor
0,00
2206
0,02
3573
00
00
0
α1
Est
imat
iva
0,42
318
10,
6041
0,42
270,
6711
0,34
540,
1442
Err
oPa
drão
0,06
9267
0,19
5338
0,12
520,
2134
0,17
740,
2286
0,12
77IC
(0,9
5)(0
,287
18;0
,559
18)
(0,6
1647
3;1
,383
527)
(0,2
726
;0,7
825)
(0,1
881
;0,7
926)
(0,3
147
;0,9
456)
(0,0
293
;0,8
202)
(0,0
095
;0,5
366)
P-va
lor
00
00
00
0
α2
Est
imat
iva
0,14
2434
0,31
3021
0,11
010,
0025
0,23
350,
1448
0,18
61E
rro
Padr
ão0,
0768
150,
3977
950,
0633
0,01
090,
138
0,20
380,
1274
IC(0
,95)
(-0,
0083
85;0
,293
253)
(-0,
4680
11;1
,094
053)
(0,0
034
;0,2
307)
(0;0
,025
7)(0
,017
;0,5
337)
(0,0
01;0
,659
8)(0
,016
4;0
,49)
P-va
lor
0,06
3705
0,43
1346
00,
003
00
0
β1
Est
imat
iva
00,
0744
920,
0243
0,00
10,
0574
0,07
690,
1043
Err
oPa
drão
0,09
7944
0,26
5192
0,02
520,
0045
0,05
470,
1222
0,10
82IC
(0,9
5)(-
0,19
2303
;0,1
9230
3)(-
0,44
6187
;0,5
9517
1)(0
,000
6;0
,096
9)(0
;0,0
075)
(0,0
013
;0,2
05)
(0,0
011
;0,4
324)
(0,0
05;0
,406
5)P-
valo
r1
0,77
8789
00,
0026
00
0
β2
Est
imat
iva
00
0,24
940,
0761
0,05
560,
1102
0,08
28E
rro
Padr
ão0,
0898
550,
1142
460,
101
0,11
830,
0641
0,16
580,
086
IC(0
,95)
(-0,
1764
22;0
,176
422)
(-0,
2243
1;0
,224
31)
(0,0
056
;0,4
113)
(0;0
,379
5)(0
,001
4;0
,204
6)(0
,000
9;0
,541
5)(0
,003
1;0
,303
6)P-
valo
r1
10
00
00
β3
Est
imat
iva
00
0,31
340,
2835
0,12
220,
1331
0,15
81E
rro
Padr
ão0,
0261
140,
0530
240,
088
0,08
440,
0968
0,16
470,
0874
IC(0
,95)
(-0,
0512
72;0
,051
272)
(-0,
1041
07;0
,104
107)
(0,1
051
;0,4
674)
(0,1
609
;0,4
807)
(0,0
077
;0,3
604)
(0,0
006
;0,5
44)
(0,0
27;0
,359
9)P-
valo
r1
10
00
00
β4
Est
imat
iva
0,07
6369
0,02
6933
0,08
70,
3246
0,11
380,
0755
0,10
99E
rro
Padr
ão0,
0334
870,
0381
670,
0651
0,20
250,
0858
0,13
950,
0966
IC(0
,95)
0,01
062
;0,1
4211
8)(-
0,04
8005
;0,1
0187
)(0
,004
3;0
,231
3)(0
,137
2;0
,654
1)(0
,006
3;0
,329
5)(0
,000
6;0
,51)
(0,0
074
;0,3
731)
P-va
lor
0,02
2576
0,48
0404
00
00
0
β5
Est
imat
iva
0,51
1274
0,17
2274
0,07
70,
2123
0,05
910,
0645
0,10
53E
rro
Padr
ão0,
0558
840,
0588
570,
1333
0,22
190,
0583
0,09
80,
0833
IC(0
,95)
(0,4
0155
1;0
,620
997)
(0,0
5671
3;0
,287
835)
(0,0
014
;0,5
065)
(0;0
,531
9)(0
,002
1;0
,218
)(0
,000
5;0
,369
1)(0
,005
4;0
,344
7)P-
valo
r0
0,00
3423
00
00
0
tE
stim
ativ
a-
2,56
4891
--
7,15
41-
7,08
29E
rro
Padr
ão-
0,13
255
--
33,5
115
-9,
0417
IC(0
,95)
-(2
,304
642
;2,8
2513
9)-
-(0
,086
8;4
5,93
11)
-(0
,478
5;2
9,06
44)
P-va
lor
-0
--
0,00
5-
0Fo
nte:
Dad
osda
Est
ação
Aut
omát
ica
deB
rasí
lia-I
NM
ET.
72
Temperatura Máxima
Para a análise frequentista segundo os modelos M1 e M2, é rejeitada a normali-dade dos resíduos (respectivamente o ajuste quanto à distribuição t), como pode ser bemvisualizado na figura (9.7). Os resíduos não aparentam comportamento autoregressivo,sendo encontradas apenas correlações espúrias. Pode-se considerar nesse caso um ajustemelhor segundo o modelo M1, dado que o qq-plot parece ser menos violado visualmente.
(a) M1 (b) M2
Figura 9.9: Temp. Máxima - Gráficos dos resíduos dos modelos
Para os modelos Bayesianos, como pode ser visto no Apêndice D, todos os mod-elos tiveram problemas no diagnóstico de Heidelberg. Em todos os casos a adoção de talmodelo parece ser bem complicada para o atual número de iterações. Para se ter idéia,apenas três parâmetros tiveram êxito no teste, sendo cadeias de modelos diferentes. Destaforma, as cadeias atuais não parecem ser muito satisfatórias para a correta estimação dasérie dos retornos.
Segundo o teste de Geweke, o único modelos em que todos os parâmetros apre-sentam convergência é o M6, isto é, com estimação via amostrador de Griddy-Gibbs cominovações normais. Apesar de não haver unanimidade entre os dois testes, o modelo M6será adotado como o mais adequado.
Por fim, a tabela (9.9) resume as estimativas para a série. Segundo ela, apesarde modelos gerados segundo algoritmo de Metropolis-Hastings possuírem erro-padrãobaixos, os mesmos não são “confiáveis” no sentido de que a cadeia não convergiu satisfa-toriamente. Assim, será comparada a eficiência na predição da volatilidade dos modelosM1 (frequentista) e M6 (bayesiana).
73
Tabe
la9.
9:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Tem
pera
tura
Máx
ima
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0E
stim
ativ
a0,
0012
0,00
0852
0,00
190,
0023
0,00
170,
3376
0,27
Err
oPa
drão
0,00
0474
0,00
0517
0,00
090,
0005
0,00
130,
2252
0,18
64IC
(0,9
5)(0
,000
269
;0,0
0213
1)(-
0,00
0164
;0,0
0186
8)(0
,001
2;0
,002
6)(0
,001
4;0
,003
2)(0
,000
6;0
,005
3)(0
,023
;0,7
995)
(0,0
227
;0,6
678)
P-va
lor
0,01
1372
0,09
9745
00
00
0
α1
Est
imat
iva
0,30
1967
0,34
3308
0,29
320,
3415
0,24
080,
1181
0,11
33E
rro
Padr
ão0,
0877
950,
1222
860,
0699
0,07
470,
1231
0,10
870,
1222
IC(0
,95)
(0,1
2959
2;0
,474
341)
(0,1
0321
2;0
,583
404)
(0,1
847
;0,4
432)
(0,2
157
;0,5
001)
(0,0
695
;0,5
099)
(0,0
049
;0,4
174)
(0,0
043
;0,4
362)
P-va
lor
0,00
0583
0,00
4994
00
00
0
β1
Est
imat
iva
0,16
9969
0,19
2434
0,04
90,
0399
0,07
630,
1107
0,10
64E
rro
Padr
ão0,
0945
610,
1163
030,
0398
0,03
680,
0811
0,10
140,
0981
IC(0
,95)
(-0,
0156
92;0
,355
629)
(-0,
0359
13;0
,420
782)
(0,0
021
;0,1
485)
(0,0
001
;0,1
276)
(0,0
039
;0,3
057)
(0,0
048
;0,3
563)
(0,0
061
;0,3
669)
P-va
lor
0,07
2265
0,09
8006
00
00
0
β2
Est
imat
iva
00
0,03
60,
021
0,09
170,
1104
0,09
13E
rro
Padr
ão0,
1682
450,
1749
310,
027
0,02
70,
0868
0,11
230,
0762
IC(0
,95)
(-0,
3303
3;0
,330
33)
(-0,
3434
58;0
,343
458)
(0,0
02;0
,095
9)(0
,000
1;0
,090
5)(0
,002
6;0
,322
4)(0
,004
4;0
,400
2)(0
,003
1;0
,257
3)P-
valo
r1
10
00
00
β3
Est
imat
iva
0,25
8063
0,31
1259
0,22
520,
1142
0,30
680,
3224
0,18
92E
rro
Padr
ão0,
1208
870,
1650
30,
1474
0,15
750,
1854
0,10
080,
1798
IC(0
,95)
(0,0
2071
6;0
,495
411)
(-0,
0127
59;0
,635
277)
(0,0
17;0
,489
5)(0
;0,4
006)
(0,0
199
;0,6
464)
(0,2
201
;0,4
012)
(0,0
011
;0,6
092)
P-va
lor
0,03
2781
0,05
9285
00
00,
1601
0
tE
stim
ativ
a-
4,57
1261
--
3,57
47-
8,26
73E
rro
Padr
ão-
0,80
6038
--
10,7
861
-10
,667
IC(0
,95)
-(2
,988
695
;6,1
5382
8)-
-(0
,036
7;2
4,64
82)
-(0
,272
7;3
5,16
28)
P-va
lor
-0
--
0-
0Fo
nte:
Dad
osda
Est
ação
Aut
omát
ica
deB
rasí
lia-I
NM
ET.
74
Temperatura Mínima
Para a análise frequentista segundo os modelos M1 e M2, em nenhum dos casosé rejeitada a normalidade dos resíduos (respectivamente o ajuste quanto à distribuiçãot), como pode ser bem visualizado na figura (9.7). Os resíduos não aparentam comporta-mento autoregressivo, sendo encontradas apenas correlações espúrias. Pode-se considerarnesse caso um ajuste melhor segundo o modelo M2, dado que o qq-plot tem um ajusteperfeito visualmente.
(a) M1 (b) M2
Figura 9.10: Temp. Mínima - Gráficos dos resíduos dos modelos
Para os modelos Bayesianos, como pode ser visto no Apêndice D, todos os mod-elos tiveram problemas no diagnóstico de Heidelberg. Apesar de um bom ajuste “visual”da convergência de cada uma das séries, os parâmetros não passaram satisfatoriamentenesse teste. Novamente, apesar de na maioria dos casos a cadeia “passar” no teste deestacionariedade, acaba sendo considerada inadequada para a estimação pontual (médiaà posteriori). Tal problema está intimamente ligado ao número de iterações adotado que,nesse caso, de fato é bem baixo.
Segundo o teste de Geweke, apenas os modelos baseados no algoritmo de Metropolis-Hastings apresentaram problemas quanto à convergência. Tal fato se deve talvez à taxa derejeição ser aparentemente muito baixa para esses modelos (cerca de 20%), o que remetenovamente à questão do número de iterações necessárias para o cálculo.
Por fim, a tabela (9.10) resume as estimativas para a série. Segundo ela, para osmodelos M1 e M2 temos a aceitação da hipótese de nulidade do parâmetro α0 ao nível
75
de 5%. Observamos que o modelo M2, apesar de aparentar ajuste perfeito anteriormente,possui amplitude menor nos intervalos de confiança, além de ser menos parcimonioso.Logo, escolhe-se o modelo M1. Quanto aos modelos Bayesianos, o modelo M6 apresentamelhores condições de convergência e ainda possui menores erros-padrões em compara-ção ao M7 e, portanto, ele será também escolhido.
76
Tabe
la9.
10:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Tem
pera
tura
Mín
ima
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0E
stim
ativ
a0,
0001
960,
0002
190,
0034
0,01
39e
-04
0,36
270,
3511
Err
oPa
drão
0,00
0113
0,00
0165
0,00
230,
0179
0,00
120,
2264
0,22
29IC
(0,9
5)(-
0,00
0026
;0,0
0041
9)(-
0,00
0104
;0,0
0054
2)(0
,000
4;0
,007
5)(0
,000
3;0
,043
4)(0
,000
2;0
,004
)(0
,020
9;0
,836
3)(0
,012
9;0
,842
3)P-
valo
r0,
0833
440,
1834
530
00
00
α1
Est
imat
iva
0,06
9993
0,08
578
0,12
880,
0961
0,13
870,
0973
0,11
36E
rro
Padr
ão0,
0206
080,
0320
190,
0483
0,03
450,
0897
0,09
210,
1129
IC(0
,95)
(0,0
2953
2;0
,110
454)
(0,0
2291
4;0
,148
646)
(0,0
495
;0,2
234)
(0,0
402
;0,1
655)
(0,0
688
;0,2
708)
(0,0
015
;0,2
909)
(0,0
037
;0,3
672)
P-va
lor
0,00
0683
0,00
7384
00
00
0
β1
Est
imat
iva
0,89
9668
0,88
2498
0,36
890,
5226
0,76
310,
231
0,21
Err
oPa
drão
0,03
296
0,05
0467
0,32
20,
3718
0,09
650,
1908
0,20
24IC
(0,9
5)(0
,834
955
;0,9
6438
)(0
,783
413
;0,9
8158
3)(0
,002
8;0
,863
4)(0
,001
7;0
,897
7)(0
,510
4;0
,896
3)(0
,002
5;0
,628
)(0
,001
1;0
,681
5)P-
valo
r0
00
00
00
tE
stim
ativ
a-
8,56
3351
--
10,4
795
-10
,648
4E
rro
Padr
ão-
2,71
2398
--
5,96
88-
17,0
753
IC(0
,95)
-(3
,237
889
;13,
8888
14)
--
(2,7
645
;26,
0562
)-
(0,1
99;6
7,67
15)
P-va
lor
-0,
0015
93-
-0
-0
Font
e:D
ados
daE
staç
ãoA
utom
átic
ade
Bra
sília
-IN
ME
T.
77
Umidade Relativa
Para a análise frequentista segundo os modelos M1 e M2, nos dois casos é rejeitadaa normalidade dos resíduos (respectivamente o ajuste quanto à distribuição t), como podeser bem visualizado na figura (9.7). Os resíduos não aparentam comportamento autore-gressivo, porém para os dois modelos os correlogramas dos resíduos ajustados aparentamcerto comportamento de médias móveis, confirmando a não-aleatoriedade requerida.
Já a adequação quanto ao tipo de distribuição, os qq-plots dos resíduos ajustadosaparentam o mesmo problema: caudas pesadas à direita, indicando que talvez M1 e M2não consigam captar perfeitamente as caudas pesadas.
(a) M1 (b) M2
Figura 9.11: Um. Relativa - Gráficos dos resíduos dos modelos
Para os modelos Bayesianos, como pode ser visto no Apêndice D, todos os mode-los tiveram problemas no diagnóstico de Heidelberg, especificamente os mesmos proble-mas quanto à adequação do uso como média à posteriori dos parâmetros α0 e β2. Já parao teste de Geweke, somente os modelos M5 e M6 não falharam quanto à convergência denenhum parâmetro da cadeia.
Apesar de um bom ajuste “visual” da convergência de cada uma das séries, osparâmetros não passaram satisfatoriamente nesse teste. Novamente, apesar de na maioriados casos a cadeia “passar” no teste de estacionariedade, acaba sendo considerada in-adequada para a estimação pontual (média à posteriori). Tal problema está intimamenteligado ao número de iterações adotado que, nesse caso, de fato é bem baixo.
Por fim, a tabela (9.11) resume as estimativas para a série. Segundo ela, ao nívelde 5% não rejeita-se a hipótese de nulidade para o parâmetro β1 dos modelo M1 e M2.
78
Portanto, segundo tais modelos, parece haver mais persistência da volatilidade diária deacordo com dois dias anterioris, considerando as estimativas altas dos parâmetros β2.
Já para os modelos bayesianos, há clara divisão quanto às estimativas segundo otipo de modelo usado: para os modelos M3, M4 e M5, o padrão de persistência alta noparâmetro β2 segue o padrão frequentista, enquanto para os modelos M6 e M8 parecerhaver menores valores de tais parâmetros, porém com um grande valor de α0. O quepode estar ocorrendo, dado os diagnósticos, é que modelos M3, M4 e M5 podem nãoestar muito certos em seguir o padrão frequentistas, no sentido de não serem muito “con-fiáveis”. Será então escolhido o modelo M6 para representar os Bayesianos.
79
Tabe
la9.
11:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Um
idad
eR
elat
iva
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0E
stim
ativ
a0,
0007
720,
0015
360,
0118
0,01
230,
0073
0,53
230,
2021
Err
oPa
drão
0,00
0387
0,00
10,
0041
0,00
830,
0067
0,23
140,
2535
IC(0
,95)
(0,0
0001
3;0
,001
532)
(-0,
0004
27;0
,003
498)
(0,0
04;0
,023
6)(0
,004
2;0
,026
1)(0
,001
6;0
,024
3)(0
,116
1;0
,897
8)(0
,009
5;0
,749
3)P-
valo
r0,
0458
0,12
4387
00
00
0
α1
Est
imat
iva
0,07
6166
0,13
5216
0,16
070,
1225
0,17
620,
0508
0,19
14E
rro
Padr
ão0,
0152
350,
0488
410,
0486
0,03
450,
1143
0,05
330,
2059
IC(0
,95)
(0,0
4625
4;0
,106
077)
(0,0
3932
3;0
,231
11)
(0,0
821
;0,2
75)
(0,0
708
;0,1
995)
(0,0
095
;0,4
168)
(0,0
032
;0,2
215)
(0,0
052
;0,6
782)
P-va
lor
0,00
0001
0,00
5632
00
00
0
β1
Est
imat
iva
00,
2171
770,
0572
0,00
880,
1352
0,04
930,
1373
Err
oPa
drão
0,02
225
0,25
1041
0,03
510,
0049
0,08
90,
0495
0,15
19IC
(0,9
5)(-
0,04
3686
;0,0
4368
6)(-
0,27
5712
;0,7
1006
6)(0
,007
7;0
,134
6)(0
,002
7;0
,021
3)(0
,011
7;0
,326
1)(0
,005
5;0
,173
4)(0
,001
6;0
,475
1)P-
valo
r1
0,38
6981
00
00
0
β2
Est
imat
iva
0,90
1196
0,60
6682
0,45
120,
4925
0,57
480,
0269
0,13
04E
rro
Padr
ão0,
0277
720,
2568
950,
1453
0,27
810,
1871
0,08
0,16
48IC
(0,9
5)(0
,846
669
;0,9
5572
2)(0
,102
299
;1,1
1106
5)(0
,059
;0,7
696)
(0,1
;0,7
719)
(0,0
762
;0,8
203)
(0,0
001
;0,3
041)
(0,0
008
;0,5
451)
P-va
lor
00,
0181
960
00
0,00
030
tE
stim
ativ
a-
5,13
6593
--
3,30
43-
14,1
877
Err
oPa
drão
-0,
9333
28-
-16
,836
-19
,643
1IC
(0,9
5)-
(3,3
0411
3;6
,969
073)
--
(0,0
409
;14,
1916
)-
(0,1
929
;75,
7154
)P-
valo
r-
0-
-0,
0098
-0
Font
e:D
ados
daE
staç
ãoA
utom
átic
ade
Bra
sília
-IN
ME
T.
80
Velocidade do Vento
Para a análise frequentista segundo os modelos M1 e M2, é rejeitada a normalidadedos resíduos (respectivamente o ajuste quanto à distribuição t), pode ser bem visualizadona figura (9.7). Os resíduos ajustados aparentam comportamento de médias móveis, sendorejeitada a aleatoriedade do mesmo. Pode-se considerar nesse caso um ajuste melhorsegundo o modelo M2, dado que o qq-plot parece ser menos violado visualmente.
(a) M1 (b) M2
Figura 9.12: Vel. Vento - Gráficos dos resíduos dos modelos
Para os modelos Bayesianos, como pode ser visto no Apêndice D, apesar doprimeiro parâmetro não aparentar problemas quanto à sua convergência (exceto para M7),todos os outros tiveram problemas no diagnóstico de Heidelberg. Em todos os casos aadoção de tal modelo parece ser bem complicada para o atual número de iterações.
Segundo o teste de Geweke, o parâmetro α1 do modelo M4 e o α1 do modelo M7não convergem, sendo portanto necessário um número maior para a devida convergênciadas cadeias.
Por fim, a tabela (9.12) resume as estimativas para a série. Segundo ela, para osmodelos Bayesianos, o mais adequado parece ser o modelo M3 que obteve bons diagnós-ticos de convergência da amostra simulada da posteriori e teve intervalos de amplitudebem menor para seus parâmetros estimados. Observa-se aqui a grande diferença entre osvalores estimados para cada um dos modelos.
81
Tabe
la9.
12:E
stim
ativ
asdo
spa
râm
etro
spa
raa
vari
ável
Vel
ocid
ade
doV
ento
segu
ndo
osse
tem
odel
osPa
râm
etro
sE
stat
ístic
asM
odel
osM
1M
2M
3M
4M
5M
6M
7α
0E
stim
ativ
a0,
5304
280,
5136
290,
2563
0,23
570,
2649
0,50
340,
4917
Err
oPa
drão
0,03
7592
0,04
942
0,01
420,
0145
0,03
220,
2793
0,27
83IC
(0,9
5)(0
,456
62;0
,604
235)
(0,4
1659
9;0
,610
658)
(0,2
287
;0,2
836)
(0,2
1;0
,261
4)(0
,205
6;0
,330
7)(0
,036
3;0
,929
8)(0
,010
8;0
,945
2)P-
valo
r0
00
00
00
α1
Est
imat
iva
0,22
6591
0,30
0913
0,00
650,
0667
0,02
230,
2258
0,38
15E
rro
Padr
ão0,
0559
790,
0776
710,
0064
0,01
520,
0217
0,27
090,
2805
IC(0
,95)
(0,1
1668
3;0
,336
499)
(0,1
4841
5;0
,453
412)
(0,0
002
;0,0
267)
(0,0
458
;0,0
942)
(0,0
002
;0,0
75)
(0,0
012
;0,8
821)
(0,0
017
;0,9
019)
P-va
lor
0,00
0052
0,00
0107
00
00
0
tE
stim
ativ
a-
6,22
8321
--
4,04
94-
12,9
666
Err
oPa
drão
-1,
5927
88-
-22
,147
4-
17,0
968
IC(0
,95)
-(3
,101
083
;9,3
5556
)-
-(0
,057
9;2
7,11
16)
-(0
,361
6;6
0,55
41)
P-va
lor
-0,
0000
92-
-0,
016
-0
Font
e:D
ados
daE
staç
ãoA
utom
átic
ade
Bra
sília
-IN
ME
T.
82
Comentário:
Apesar de não explicitadas aqui, foram calculados os Fatores de Bayes (que nadamais são que razões das posterioris avaliadas nos parâmetros estimados por dois modelosdiferentes) dos modelos M3 a M5 dois a dois, formando uma tabela 5x5 para cada sérieclimatológica. Neles os resultados apontam para as escolhas dos modelos citados em cadasituação.
Além disso, os modelos frequentistas selecionados para a etapa da previsão tam-bém foram os “melhores” baseados nos critérios AIC e BIC, omitidas aqui as tabelas quepodem ser vistas também no apêndice.
9.2.3 Previsão de volatilidade
Após a seleção de uma metodologia ótima e a subseqüente escolha dos parâmetrosdo modelo para cada uma das séries, avaliou-se de maneira prática tal escolha de acordocom a previsão da volatilidade para a série para o mês de dezembro de 2008, além dasprevisões de volatilidade da série.
Em cada uma das figuras abaixo encontram-se três gráficos para a previsão devolatilidade dos dois modelos (um clássico e outro bayesiano) para cada uma das sériesde retornos analisada. Tais modelos foram considerados “melhores” que os outros para ocorreto ajuste das séries.
Salienta-se novamente que o intuito é a previsão da volatilidade, não dos valorespontuais das séries dos retornos e, consequentemente, das séries originais. Para uma abor-dagem nesse sentido, apesar de que fora mostrado que o valor esperado para os retornos(estimativas pontuais) é nulo, alguns autores analisaram os modelos via bootstrap parauma boa estimativa dos retornos (ver Pascual, Romo e Ruiz, 2005).
Um comentário importante é com relação aos últimos gráficos analisados. Em to-dos há três linhas que indicam as estimativas das médias de cada uma das volatilidades,sendo construídos intervalos de credibilidade para cada uma das volatilidades (para detal-hes, ver Bawens e Lubrano, 1998) da seguinte maneira: Para cada um dos conjuntos deparâmetros resultantes da amostra final à posteriori são calculadas as volatilidades até otempo desejado. Em seguida, ter-se-á vários valores estimados para uma mesma volatil-idade, sendo então calculadas médias e intervalos de confianças segundo a metodologiabayesiana.
Na verdade, essa metodologia é baseada no uso da distribuição preditiva à pos-teriori, isto é, as volatilidades condicionadas ao conjunto de parâmetros passarão a servariáveis aleatórias (por ser função de parâmetros aleatórios), sendo então passíveis deconstrução de intervalos de credibilidade. Para maiores detalhes ver apêndice.
83
Já no caso frequentista, não é possível a construção de intervalos de confiançada mesma forma, sendo que poder-se-ia recorrer à intervalos bootstrap (ver novamentePascual, Romo e Ruiz, 2005), porém não abordada aqui tal metodologia.
Evaporação
Para a série dos retornos de evaporação, o primeiro gráfico aparenta um bomcomportamento das volatilidades estimadas segundo o modelo frequentista. O modeloBayesiano também conseguem captar bem os períodos de flutuação, porém quando estessão grandes não consegue acompanhar tamanho crescimento, como no caso frequentista.
Analisando as previsões de volatilidade para o mês de dezembro de 2008, o mo-delo frequentista conseguiu prever bem valores para as semanas subsequentes, enquantoo modelo bayesiano pecou muito nas estimativas. Para a evaporação, os valores espera-dos (em verde) sequer estavam entre os limites de credibilidade para as estimativas dasvolatilidaes.
Figura 9.13: Gráficos das previsões de volatilidade para os retornos da Evaporação
Para a série dos retornos de temperatura máxima, já no primeiro gráfico é obser-vada grande disparidade entre os dois modelos de previsão. O modelo frequentista (M1)até “acompanha” muito bem a série dos retornos, porém a variabilidade desta última émuito maior que o previsto no modelo. Pelo segundo gráfico da figura, há maior com-preensão do comportamento: enquanto os valores esperados está muito altos, o modeloprevê pouca variabilidade naquele mês, o que não condiz com a realidade.
84
Já o modelo Bayesiano tem um comportamento diferente, superestimando a varia-bilidade condicional da série, como observado no primeiro gráfico (linha azul), sendo istobem refletido no terceiro gráfico, onde, apesar de não muito nítido na figura (grande vari-abilidade, impossibilitando boa saída gráfica do software), os valores esperados (volatili-dades ajustadas segundo dados originais) estão entre os limites de confiança estabelecidospara a volatilidade, porém a previsão, apesar de estimar bem a grande volatilidade para operíodo no meio do mês, ainda peca na superestimação.
Figura 9.14: Gráficos das previsões de volatilidade para os retornos da Temperatura Máx-ima
O ajuste da série de temperatura mínima foi bem satisfatório para os dois mod-elos considerados. Apesar de visualmente o modelo bayesiano (linha azul do primeirográfico) na figura (9.15) ser um pouco acima dos valores dos retornos, este parece estarapenas superestimando um pouco a variabilidade, porém capta muito bem a “forma” davolatilidade. E isto se reflete no último gráfico onde, apesar da volatilidade “esperada”(observada com os valores originais) para o mês de dezembro estarem entre os limites deconfiança, este tem amplitude muito grande. A própria volatilidade “média” não captacorretamente a variabilidade do período, indicando novamente uma superestimação.
Já no método clássico, há um ajuste muito bom da variabilidade do período se-gundo o modelo M1. observa-se que a linha vermelha (previsão) está um pouco abaixodo esperado, mas no geral indicou ótima previsibilidade da baixa variabilidade do período.
A análise da volatilidade para a série dos retornos da série de Umidade relativaé bem parecida com a encontrada anteriormente para temperatura mínima. Novamente
85
Figura 9.15: Gráficos das previsões de volatilidade para os retornos de Temperatura Mín-ima
o modelo Bayesiano superestima a volatilidade da séries devido ao valor de α0 ser bemalto para o segundo modelo (0.5323), enquanto possui valor bem baixo para o primeiromodelo (0.000772). Porém novamente tanto o modelo M1, quanto o M6 conseguemcaptar as flutuações da série.
O comportamento constante das séries após os primeiros m = max(p,q) valoresé claramente devido ao caráter autoregressivo do modelo, já discutido na parte teórica.Deste modo, as previsões serão “úteis” apenas para esses lags. Os modelos para a previsãoda série dos retornos da umidade relativa podem ser considerados satisfatórios pois até olag especificado (no caso m = max(1,2) = 2) apresentam boa acurácia.
Finalmente, para a os retornos da variável velocidade do vento (figura (9.17)abaixo) muitos problemas ocorreram. O modelo Bayesiano pecou muito na estimaçãoda volatilidade da série para todo o período, tendo muito dificuldade na captação das “nu-ances” da série ou, explicitamente, da heterocedasticidade. Para o modelo clássico, apesarda razoável captação dos picos de variabilidade, notoriamente no final de 2007, em geral,por se tratar de um modelo do tipo ARCH(1), não conseguiu prever bem a variabilidadedo primeiro dia do mês e, como consequência (no caso) superestimou a variabilidade detodo o período subsequente (explicado anteriormente).
O modelo M3 (bayesiano) ainda conseguiu fazer boa previsão da volatilidadequando considerada a magnitude da variabilidade dos intervalos de credibilidade paracada uma das volatilidades estimadas, porém, devido ao caráter autoregressivo de ordem
86
Figura 9.16: Gráficos das previsões de volatilidade para os retornos de Umidade Relativa
1 (ARCH(1)), novamente uma supertimação do primeiro valor já impossibilita a boa pre-visão para o resto do mês que, aliás, não demonstrou grande variabilidade.
Comentário:
Quanto à significância estatística das previsões, como o intuito é prever as volatil-idades, há problema na comparação de valores observados vs valores previstos, dadoque as volatilidades não são observáveis. Caso fossem, outras medidas de diagnósticopoderiam ser adotadas para melhor exemplificação da significância estatística de tal difer-enciação.
9.3 Aspectos Computacionais
Todo a parte computacional do trabalho foi realizada no software livre R, com autilização de vários pacotes tanto para o manejo dos dados, quanto para o cômputo deinúmeras estatísticas.
Na parte frequentista, isto é, para as estimativas dos parâmetros dos modelos M1 eM2 foi utilizada a função “garchFit” do pacote fGarch 1. Já para os modelos Bayesianos(M3 a M7), como não há rotinas prontas exceto para o modelo GARCH(1,1) com ino-
1 disponibilizado em <http://cran.r-project.org/web/packages/fGarch/index.html>
87
Figura 9.17: Gráficos das previsões de volatilidade para os retornos da Velocidade doVento
vações t, seguindo uma metodologia específica de prioris 2 todos os passos do algoritmode Metropolis-Hastings e do Amostrador de Griddy-Gibbs, além de todas as saídas gráfi-cas, foram programados pelo autor.
Essencialmente foram escritas duas funções mestras, uma para cada algoritmo,onde o usuário entra com a série a ser modelada, a dimensão (p,q) do modelo GARCH aser ajustado, o tipo de distribuição suposta (permitido apenas a normal como default e a tde student) e finalmente as prioris para o modelo, obrigatoriamente um vetor de tamanho(p+q+1). 3
Muitos problemas computacionais ocorreram, levando a colocação de estruturascondicionais (if, else, while, etc.) na programação, aumentando bastante o tempo deprogramação. A tabela (9.13) abaixo mostra o tempo médio de processamento do sistema(em minutos) para cada um dos algoritmos. Salienta-se novamente que para cada variável,para cada modelo foram inicialmente geradas cadeias de tamanho 6000, para verificar adependência de valores iniciais e em seguida uma cadeia definitiva de tamanho 10.000(18.000, para a simulação) foi utilizada para a análise.
Os dados foram retirados com a função “system.time” do R, sendo baseados emum computador padrão: processador Dual Core 2 duo, com 3gb de memória RAM.
2 ver pacote bayesGARCH, em http://cran.r-project.org/web/packages/bayesGARCH/index.html3 Devido à extensão, as rotinas foram omitidas dos Apêndices deste trabalho, sendo possível, no entanto, o
acesso às mesmas via contato com autor pelo e-mail: <[email protected]>
88
Tabela 9.13: Tempo de processamento, em minutosVariável M3 M4 M5 M6 M7
Evaporação 22,53 24,40 24,68 756,01 814,30Insolação 26,12 28,24 25,82 823,36 847,84Temp. Máxima 24,46 25,23 28,01 757,67 818,15Temp. Mínima 22,97 24,90 25,21 756,73 815,30Um. Relativa 26,36 27,58 29,09 760,23 815,06Vel. Vento 25,25 25,99 27,10 757,93 814,79
Assim, o tamanho esforço computacional para a geração de todas as amostras decerta forma “atrasou” o trabalho, sendo que o número de variáveis utilizadas (7, ao con-tar com o estudo de simulação), impossibilitar a geração de cadeias maiores devido aotempo e a disponibilidade de computadores com boa capacidade (memória RAM, especi-ficamente).
Há vários softwares mais compactos para análise mais “robusta” de métodos MCMC,como o WinBugs (BUGS - Bayesian Inference Using Gibbs Sampler) 4 e o JAGS (JustAnother Gibbs Sampler) 5. Ambos começaram a ser utilizados para fins de compara-ção com as rotinas escritas em R, porém possuem problema de tempo de processamentoquando as rotinas são muito grandes, inviabilizando comparações.
4 disponível em: <http://www.mrcbsu.cam.ac.uk/bugs/winbugs/contents.shtml> , acessado em 15 denovembro de 2009
5 disponível em: <http://calvin.iarc.fr/ martyn/software/jags/>, acessado em 15 de novembro de 2009
10 Conclusão e Trabalhos Futuros
Foram desenvolvidos modelos de previsão para a volatilidade de várias sériesde retornos de variáveis climatológicas. Tais séries originais, apesar de não possuíremcomportamento análogo ao de séries financeiras comumente modeladas com modelos devolatilidade estocástica em geral, quando utilizada a transformação para os retornos, éaparente a presença de “blocos de variação” na séries dos retornos.
Foi analisada cada uma das séries segundo a ótica frequentista e a bayesiana deestimação, sendo divididos em 7 modelos possíveis, nomeados de M1 a M7 bem especi-ficados no capítulo sobre as previsões. Foram ao total 7 variáveis analisadas, sendo queuma delas era uma série simulada pelo próprio autor utilizando o software R. Para estasérie foram testadas e aprovadas as ferramentas de estimação tanto disponíveis para o soft-ware utilizado (modelos M1 e M2 - frequentistas), quanto as criadas pelo autor (modelosM3 a M7 - bayesianos).
Tais modelos são bem complexos no sentido da estimação e bem diferentes unsdos outros. Enquanto os modelos frequentistas utilizaram o método de Newton-Raphsonpara a estimação dos parâmetros, os métodos Bayesianos utilizaram métodos MCMCpara a geração de amostras da distribuição à posteriori para consequente estimação dosparâmetros.
Apesar de na literatura ser bem evidente a superioridade de método bayesiano deestimação dos modelos GARCH ao produzirem intervalos de credibilidade de menor am-plitude quando comparados aos intervalos de confiança frequentistas, no presente trabalhonão verificou-se nitidamente tal superioridade quando da previsão da volatilidade.
E tal fato pode estar intimamente ligado às amostras da distribuição à posteriorimarginal de cada um dos parâmetros. Como visto, a convergência foi extremamentecomplicada em muito modelos para a maioria das variáveis analisadas. O ideal seriatrabalhar com um número muito maior de iterações (da ordem talvez de 50.000 ou mais),dado que a construção dos algoritmos de Metropolis-Hastings e o amostrador de Griddy-Gibbs gatantem a convergência assintótica para a distribuição limite (posteriori). Destemodo, seriam produzidas amostras confiáveis de acordo com os testes diagnósticos usados
90
e só assim estimativas adequadas poderiam ser utilizadas para a comparação adequada dosmétodos.
Quanto aos modelos Bayesianos, os dois algoritmos tiveram muita dificuldadede produzir boas amostras, devido à taxa de rejeição de todas as cadeias serem muitobaixas para o algoritmo de Metropolis-Hastings, enquanto devido à grande variabilidadedas amostras para o amostrador de Griddy-Gibbs. No geral os dois modelos caracteri-zaram diferentes séries modeladas, sendo que para o segundo sempre resultavam grandesestimativas dos choques iniciais da volatilidade caracterizado pelo parâmetro α0, super-estimando a volatilidade inteira da série.
Foram explorados possíveis combinações de modelos na tentativa de explicar avariabilidade das séries dos retornos climatológicos extremos. Porém, muitas outras al-ternativas podem ser exploradas no presente contexto.
Inicialmente, o uso de outros reparametrizações do GARCH poderiam ser uti-lizados como: o EGARCH (Nelson, 1991) ou GARCH exponencial que possui eficácianotória na literatura, no sentido de analisar retornos positivos e retornos negativos deforma diferente; o IGARCH ou GARCH integrado, que leva em conta diferenciaçõesdentro do modelo original; há ainda inúmeros modelos, como APARCH, FIAPARCH,etc. que podem ser alternativas viáveis (para maiores detalhes, ver Engle (1995) ).
Ainda no contexto da heterocedasticidade, os modelos de volatilidade estocástica(MVE) desenvolvidos por Jacquier et. al. (1994), que levam em consideração um compo-nente aleatório à variância condicional pode ser alternativa, porém deve ser observado ocaráter multiparamétrico do modelo que aumenta a complexidade dos cálculos, tornandotalvez inviável a programação manual dos métodos, como foi feito neste trabalho. Alter-nativa às distribuições supostas pode ser a distribuição generalizada dos erros (GEV) queé bastante utilizada na teoria dos valores extremos (TVE).
Por fim, a adoção de prioris informativas pode também ser alternativa aos métodosbayesianos, apesar da preferência quanto ao uso de prioris como foram utilizadas que nãoparecem interferir abruptamente nos resultados observados, devido ao seu caráter não-informativo.
Um ponto importante foram os critérios de seleção dos modelos, feito apenas se-gundo os critérios AIC e BIC. Para os modelos bayesianos existem outras formas de se-leção ótima de modelos que podem ser utilizadas, pois neles se leva em consideração umainformação à priori dos modelos, tratando portanto os mesmo como aleatórios. Nessecontexto, há destaque para o método MCMC que possibilita a seleção de um modelodentre vários outros, possibilitanto inclusive a mudança de dimensões, como o modelosRJMCMC ou MCMC com saltos reversíveis. Este modelos começou a ser escrito neste
91
trabalho 1, porém, para esse trabalho, não houve tempo hábil para a estimação do extensonúmero de combinações de parâmetros e modelos utilizados.
Finalmente, para as variáveis temperatura mínima, evaporação e Umidade houvebom ajuste considerando as previsões de volatilidade da série. Enquanto para as variáveistemperatura máxima e Velocidade do Vento, os modelos não se ajustam adequadamenteaos dados, como considerado anteriormente.
1 Atualmente as rotinas foram concluídas e estão implementadas no R, auxiliando diariamente a previsãode volatilidade no INMET
Referências Bibliográficas
[1] J. ALBERT. Bayesian Computation with R. Springer, 2007.
[2] L. BAUWENS and M. LUBRANO. Bayesian inference on garch models using thegibbs sampler. Econometrics Journal, 1:23–46, january 1998.
[3] P. J. BICKEL and K.A. DOKSUM. Mathematical Statistical. Basical Ideas and
Selected Topics. Peter Jonh Bickel, 1977.
[4] H. BOLFARINE and M. C. SANDOVAL. Introdução à Inferência Estatística. So-ciedade Brasileira de Matemática, 2003.
[5] T. BOLLERSLEV. Generalized autoregressive conditional heteroskedasticity. Jour-
nal of Econometrics, 31:307–327, 1986.
[6] T.; BOLLERSLEV, R. Y. Chou, and K. F. Kroner. Arch modelling in finance: Areview of the theory and empirical evidence. Journal of Econometrics, 52(1):5–59,1992.
[7] T.; BOLLERSLEV, CHOU R. Y.;, and KRONER K. F. Arch modeling in finance:A review of the theory and empirical evidence. Journal of Econometrics, 52:5–59,1992.
[8] P. BURNS. Robustness of the ljung-box test and its rank equivalent. Correlations
and Volatilities of Asynchronous Data, The Journal of Derivatives, 5:7–18, 2002.
[9] G. CASELLA and E. I. GEORGE. Explaining the gibbs sampler. The American
Statistician, 46(3):167–174, august 1992.
[10] G. CASELLA and C.P ROBERT. Monte Carlo Statistical Methods. Spring tests instatistics. Springer Verlag, Nova York, USA, 2 edition, 2004.
[11] R. F. ENGLE. Autoregressive conditional heteroskedasticity with estimates of thevariance of united kingdom inflation. Econometrica, 50:987–1007, 1982.
93
[12] R. F. ENGLE and BOLLERSLEV T. Modeling the persistence of conditional vari-ances. Econometric Reviews, 5:1–50, 1986.
[13] R.F.; ENGLE, D.; Lilien, and R. Robins. Estimating time-varying risk premia in theterm structure: The arch-m model. Econometrica, 55:251–276, 1987.
[14] W.A. FULLER. Introduction to Statistical Time Series. New York, 1976. Editor:John Wiley and sons.
[15] D. GAMERMAN. Markov Chain Monte Carlo: Stochastic simulation for Bayesian
inference. Chapman & Hall, 1997. London.
[16] A. GELMAN, J. B. CARLIN, H. S. STERN, and D. B. RUBIN. Bayesian Data
Analysis. Chapman & Hall, 2 edition, 1996. CRC Press.
[17] P.J. GREEN. Reversible jump markov chain monte carlo computation and bayesianmodel determination. Biometrika, 82(4):711–732, june 1995.
[18] M. GREENE. Econometric Analysis. Prentice Hall, 6 edition, august 2007. ScottE. Page.
[19] P. G. HOEL, C. S. PORT, and STONE C. J. Introduction to Stochastic Processes.Waveland Press, December 1986.
[20] C. M. MORETTIN, P. A; TOLOI. Análise de séries temporais. Edgard Bluche,2004. São Paulo.
[21] BROCKWELL P.J. and DAVIS R.A. Introduction to Statistical Time Series and
Forecasting. Springer, 2 edition, 2002. New York.
[22] P. M. ROBINSON and P. ZAFFARONI. Pseudo-maximum likelihood estimation ofarch models. The Annals of Statistics, 34(3):1049 – 1074, 2006.
[23] R. A. ROSALES. Mcmc for hidden markov models incorporating aggregation ofstates and liltering. Bulletin of Mathematical Biology, 66(5):1173 – 1199, 2004.
[24] A. A. SALLES. Um modelo garch e a correlação implícita: associação dos retornosdo ibovespa com outros índices americanos. In Anais do XXII Encontro Nacional
de Engenharia de Produção, 2002.
[25] POLITIS D. N. VRONTOS I. D., DELLAPORTAS P. Full bayesian inference forgarch and egarch models. Journal of Business & Economic Statistics, 18(2):187–198, april 2000.
94
APÊNDICE A - Demonstrações das propriedades dosmodelos ARCH e GARCH
Usando a independência assumida entre a variância condicional σt e o ruído brancozt , mostra-se que a esperança condicional é constante e igual a zero:
E[εt |ψt−1] = E[σtzt |ψt−1]
= E[σt |ψt−1] ·E[zt |ψt−1]
= σt ·0
= 0
Usando a propriedade acima tem-se que a esperança incondicional também é con-stante e igual a zero.
E[εt ] = E [E(εt |ψt−1)]
= E[0]
= 0
Como propriedade fundamental, a variância condicional do modelo é igual a σ2t ,
isto é, depende do tempo. Usando novamente a suposição de independência entre σt e zt ,chega-se na expressão da variância condicional:
Var(εt |ψt−1) = E[ε2t |ψt−1]
= E[σ2t z2
t |ψt−1]
= E[σ2t |ψt−1] ·E[z2
t |ψt−1]
= σ2t ·1
= σ2t
Para a prova da variância incondicional é usada a hipótese de estacionariedadeassumida também. Primeiramente observa-se que:
Var(εt) = E[ε2t ] = E
[E(ε2
t |ψt−1)]= E[σ2
t ] , (.1)
agora, usando a equação anterior e a estacionariedade suposta pelo modelo:
95
E[σ2t ] = E[α0 +
q
∑i=1
αiε2t−i]
= E[α0]+q
∑i=1
αiE[ε2t−i]
= α0 +q
∑i=1
αiE[ε2t ]
= α0 +q
∑i=1
αiVar(εt) ,
logo,
Var(εt) =α0
1−∑qi=1 αi
, (.2)
com 0 ≤ ∑qi=1 αi < 1 sendo a restrição de não-negatividade da variância incondi-
cional.Assumindo a estacionariedade do processo, mostra-se que a função de autocovar-
iância é nula, isto é, isto é, a série dos erros é independente entre si. Como será visto, talcaracterística é importante na identificação e verificação do ajuste de modelos ARCH.
γε(h) = Cov(εt ,εt+h)
= E[εt · εt+h]
= E [E(εt · εt+h|ψt+h−1)]
= E [εt ·E(εt+h|ψt+h−1)]
= E [εt ·E(σt+hzt+h|ψt+h−1)]
= E [εt ·E(σt+h|ψt+h−1) ·E(zt+h|ψt+h−1)]
= E[εt ·σt+h ·0]
= E[0]
= 0 ,
Agora, as expressões da variância incondicional para os modelos GARCH. Primeira-mente observa-se que:
Var(εt) = E[ε2t ] = E
[E(ε2
t |ψt−1)]= E[σ2
t ] , (.3)
e, usando a equação anterior e a estacionariedade suposta pelo modelo:
96
E[σ2t ] = E[α0 +
q
∑i=1
αiε2t−i +
p
∑j=1
β jσ2t− j]
= E[α0]+q
∑i=1
αiE[ε2t−i]+
p
∑j=1
β jE[σ2t− j]
= α0 +q
∑i=1
αiE[ε2t ]+
p
∑j=1
β jE[σ2t ]
= α0 +q
∑i=1
αiVar(εt)+p
∑j=1
β jVar(εt) ,
Obtendo, então, a expressão da variância incondicional:
Var(εt) =α0
1−(
∑qi=1 αi +∑
pj=1 β j
) , (.4)
sendo a restrição para a não-negatividade da variância incondicional dada por:k
∑i=1
(αi +βi) < 1 , onde k = maxp,q.
Agora, a representação dos modelos GARCH segundo com o ARMA. Supondo anormalidade, observa-se que:
ηt = ε2t −σ
2t (.5)
E usando o fato de que a soma dos quadrados de n variáveis aleatórias i.i.d comdistribuição normal padrão possui uma distribuição de qui-quadrado com n graus de liber-dade e sendo k = maxp,q, pode-se escrever os erros quadráticos da seguinte forma:
ε2t = σ
2t +σ
2t (z
2t −1)
= α0 +m
∑j=1
(α j +β j)ε2t− j +ηt−
q
∑i=1
βiηt−i ,
Desta maneira, sendo os erros εt pertencentes a um modelo GARCH(p,q), oquadrados dos mesmos, ε2
t , segue um modelo auto-regressivo de médias móveis ARMA(p,q),porém com erros não-gaussianos ηt .
97
APÊNDICE B - Métodos via simulação estocástica
De maneira análoga aos métodos de Monte Carlo vistos anteriormente, os méto-dos a serem mostrados aqui são utilizados para a simulação de valores de distribuiçõesonde tal tarefa é complexa. Inicialmente são geradas observações de uma distribuiçãoauxiliar (como na amostragem via função de importância), nomeada aqui por q(x). Apósesta etapa, utiliza-se um mecanismo de correção para tornar a amostra obtida ao menosrepresentativa da função de interesse p(x).
Método de RejeiçãoÉ considerada a existência de uma constante K, finita, que atenda a p(x) < Kq(x).
Deste modo, após a geração de um valor θ′ da distribuição auxiliar, tal valor é aceitocomo sendo uma amostra da função de interesse com uma probabilidade p(x)/Kq(x).Caso contrário, não é aceito e são geradas novas amostras.
Pode-se mostrar (ver Gamerman, 1997) que o algoritmo funciona mesmo quandonão se conhece completamente a função de interesse e ainda que tal mecanismo efetiva-mente gera valores de p(x).
Para o presente contexto, sendo p(x) = π(θ), isto é, a densidade a posteriori,e ainda q(x) = p(θ) a distribuição à priori, a razão de aceitação será agora dada porπ(θ)/K p(θ). Logo, l(θ)p(θ)/K p(θ) e então a constante K deve satisfazer l(θ) < K. As-sim, sendo K o valor máximo da função de verrosimilhança, isto é, K = l(θ), onde θ é oEstimador de Máxima Verossimilhança (EMV).
Finalmente, chega-se ao seguinte algoritmo:
1. Simular um valor θ′ da distribuição à priori P(θ);
2. Simular um valor u da distribuição uniforme no intervalo (0,1);
3. Aceitar θ′ como valor da posteriori π(θ), caso u < π(θ′)/π(θ) e voltar ao primeiropasso;
Porém com os possíveis problemas do método:
• problemas quanto à maximização da função de verossimilhança;
• Caso as distribuições à priori e à posteriori sejam muito diferentes, o algoritmo correrisco de baixa aceitação sendo, portanto, maior o número de iterações requeridospara obtenção de uma amostra de tamanho razoável;
98
Método de Reamostragem ponderadaConsiste em, após ser gerada uma amostra de q, calcular-se os seguintes pesos:
wi =p(xi)/q(xi)
∑nj=1 p(x j)/q(x j)
i = 1,2, . . . ,n
Então, escolhe-se uma subamostra de tamanho m < n, onde cada uma das obser-vações iniciais xi possui probabilidade wi de ser reamostrada. Para o contexto Bayesiano,seja p(x) = π(θ) e ainda q(x) = p(θ)
wi =π(θi)/p(θi)
∑nj=1 π(θ j)/p(θ j)
=l(θi)
∑nj=1 l(θi)
i = 1,2, . . . ,n
Em relação ao método anterior, não é necessário o conhecimento da constante Kmencionada, não sendo, portanto, necessário a maximização da função de Verossimil-hança. Além disso, não há "perda"de amostra como acontece nos algoritmos de rejeição.Porém, apresentam-se nesse passo dois problemas:
• Não há garantia quanto às observações pertencerem à função alvo, apenas aproxi-madamente.
• Caso as distribuições à priori e à posteriori sejam muito diferentes, apenas poucosvalores gerados da priori terão alta probabilidade de aparecerem na reamostra (Gamer-man, 1997)
99
APÊNDICE C - Critérios de escolha AIC e BIC para asseis variáveis
Tabela .1: Evaporação - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 1,4612(1) 1,464(2) 1,4734(6) 1,474(10) 1,4869(28) 1,4717(5)
1,4742(1) 1,4835(2) 1,4994(5) 1,5066(9) 1,5259(19) 1,5173(15)
2 1,471(3) 1,4739(8,5) 1,4762(12) 1,4769(16) 1,4898(29) 1,4746(11)1,4905(3) 1,4999(6) 1,5088(10) 1,5159(14) 1,5353(24) 1,5267(20)
3 1,4711(4) 1,4739(8,5) 1,4768(15) 1,4796(20) 1,482(21) 1,4775(17)1,4971(4) 1,5065(8) 1,5158(13) 1,5252(18) 1,5341(21,5) 1,5361(25)
4 1,4737(7) 1,4765(14) 1,4794(19) 1,4822(23) 1,4848(25) 1,483(24)1,5062(7) 1,5156(12) 1,5249(17) 1,5343(23) 1,5434(26) 1,5481(28)
5 1,4763(13) 1,4792(18) 1,4821(22) 1,4849(26) 1,4986(30) 1,4858(27)1,5154(11) 1,5248(16) 1,5341(21,5) 1,5435(27) 1,5637(30) 1,5574(29)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
Tabela .2: Insolação - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 2,1458(30) 1,8168(21) 1,7897(16) 1,765(11) 1,7358(6) 1,6544(5)
2,1654(30) 1,8429(19) 1,8223(13) 1,804(9) 1,7814(6) 1,7065(2)
2 2,1097(28) 1,8197(22) 1,7926(17) 1,7678(12) 1,7376(7) 1,6454(1)2,1357(28) 1,8522(21) 1,8316(15) 1,8134(11) 1,7897(7) 1,704(1)
3 2,1098(29) 1,8214(23) 1,7946(18) 1,7707(13) 1,7416(8) 1,6483(2)2,1424(29) 1,8604(22) 1,8401(18) 1,8228(14) 1,8001(8) 1,7134(3)
4 2,0694(27) 1,9506(26) 1,7999(19) 1,7739(14) 1,7433(9) 1,6512(3)2,1085(27) 1,9962(26) 1,852(20) 1,8324(16) 1,8084(10) 1,7227(4)
5 1,9212(24) 1,922(25) 1,8024(20) 1,7743(15) 1,7463(10) 1,654(4)1,9667(24) 1,9741(25) 1,861(23) 1,8394(17) 1,8179(12) 1,7321(5)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
100
Tabela .3: Temperatura Máxima - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 -2,6917(30) -2,6999(18,5) -2,7018(14,5) -2,7096(1) -2,7067(5) -2,7057(7)
-2,6722(2) -2,6739(1) -2,6692(6) -2,6706(3) -2,6612(10) -2,6536(15)
2 -2,6958(26) -2,6984(22,5) -2,6989(20) -2,7072(4) -2,7046(9) -2,7037(11)-2,6697(4) -2,6658(7) -2,6599(13) -2,6617(9) -2,6526(16,5) -2,6451(21)
3 -2,6935(29) -2,6954(27) -2,696(25) -2,7044(10) -2,7018(14,5) -2,7008(17)-2,661(11) -2,6564(14) -2,6504(20) -2,6523(18) -2,6432(23) -2,6357(25)
4 -2,7086(2) -2,706(6) -2,7031(12) -2,7013(16) -2,6999(18,5) -2,6984(22,5)-2,6696(5) -2,6604(12) -2,651(19) -2,6427(24) -2,6348(26) -2,6269(28)
5 -2,7076(3) -2,7047(8) -2,7019(13) -2,6986(21) -2,6971(24) -2,6951(28)-2,662(8) -2,6526(16,5) -2,6433(22) -2,6336(27) -2,6255(29) -2,617(30)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
Tabela .4: Temperatura Mínima - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 -2,2525(30) -2,2984(1) -2,2957(2) -2,2954(3) -2,2949(4,5) -2,2917(10)
-2,2395(17) -2,2789(1) -2,2697(2) -2,2628(4) -2,2558(7) -2,2462(12)
2 -2,2588(29) -2,2949(4,5) -2,2928(8) -2,2934(6) -2,2933(7) -2,2897(14)-2,2393(18,5) -2,2688(3) -2,2603(5) -2,2544(8) -2,2477(11) -2,2376(20)
3 -2,2665(26) -2,2927(9) -2,2907(12) -2,2906(13) -2,2914(11) -2,2889(16)-2,2405(15) -2,2601(6) -2,2516(9) -2,245(13) -2,2393(18,5) -2,2304(24)
4 -2,2644(28) -2,2891(15) -2,2871(19) -2,2877(18) -2,2886(17) -2,2861(20)-2,2318(23) -2,25(10) -2,2415(14) -2,2357(21) -2,23(25) -2,221(28)
5 -2,2657(27) -2,2858(21) -2,2844(23) -2,2843(24) -2,2853(22) -2,2832(25)-2,2266(26) -2,2403(16) -2,2323(22) -2,2257(27) -2,2203(29) -2,2116(30)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
APÊNDICE D - Diagnóstico de Heidelberg para as seisvariáveis
101
Tabela .5: Umidade Relativa - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 -0,5784(30) -0,6406(18) -0,6519(2) -0,6453(12) -0,6531(1) -0,6504(4)
2,1654(30) 1,8429(19) 1,8223(13) 1,804(9) 1,7814(6) 1,7065(2)
2 -0,591(29) -0,6376(22) -0,649(7) -0,6425(16) -0,6503(5) -0,6516(3)2,1357(28) 1,8522(21) 1,8316(15) 1,8134(11) 1,7897(7) 1,704(1)
3 -0,5993(27,5) -0,6345(24) -0,6465(9) -0,6396(19) -0,6474(8) -0,6491(6)2,1424(29) 1,8604(22) 1,8401(18) 1,8228(14) 1,8001(8) 1,7134(3)
4 -0,5993(27,5) -0,6329(25) -0,6455(11) -0,6381(21) -0,6445(13) -0,6463(10)2,1085(27) 1,9962(26) 1,852(20) 1,8324(16) 1,8084(10) 1,7227(4)
5 -0,6323(26) -0,6387(20) -0,6427(15) -0,6351(23) -0,6417(17) -0,6442(14)1,9667(24) 1,9741(25) 1,861(23) 1,8394(17) 1,8179(12) 1,7321(5)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
Tabela .6: Velocidade do Vento - Critérios de escolha AIC e BICp q
0 1 2 3 4 51 2,4212(1) 2,4239(3) 2,4267(5) 2,4363(20) 2,4325(13,5) 2,4299(10)
2,4407(1) 2,4499(3) 2,4593(5,5) 2,4753(10) 2,4781(13,5) 2,482(15)
2 2,4238(2) 2,4267(5) 2,4295(7) 2,4391(25) 2,4354(17,5) 2,4328(15)2,4498(2) 2,4592(4) 2,4686(7,5) 2,4847(16) 2,4874(18) 2,4913(20)
3 2,4267(5) 2,4296(8) 2,4324(11,5) 2,442(29) 2,4382(22) 2,4356(19)2,4593(5,5) 2,4686(7,5) 2,478(12) 2,4941(21) 2,4968(23,5) 2,5007(25)
4 2,4297(9) 2,4325(13,5) 2,4354(17,5) 2,4383(23) 2,4411(27) 2,4385(24)2,4687(9) 2,4781(13,5) 2,4875(19) 2,4968(23,5) 2,5062(27) 2,5101(28)
5 2,4324(11,5) 2,4352(16) 2,4381(21) 2,441(26) 2,4441(30) 2,4413(28)2,4779(11) 2,4873(17) 2,4967(22) 2,5061(26) 2,5157(29) 2,5194(30)
Fonte: Dados da Estação Automática de Brasília - INMET.
Obs.1: O valor apresentado entre parênteses é o rank correspondente.
Obs.2: Para cada parâmetro, o primeiro valor denota o critério AIC e o segundo o critério BIC.
APÊNDICE D - Diagnóstico de Geweke para as seisvariáveis
102
Tabela .7: Evaporação - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Passou Passou Passou Passou Passou
Passou Passou Passou Passou Passou
α1 Passou Falhou Passou Passou PassouFalhou NA Passou Falhou Falhou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
Tabela .8: Insolação - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Passou Falhou Passou Passou Passou
Passou NA Passou Falhou Falhou
α1 Passou Falhou Passou Passou PassouPassou NA Passou Passou Falhou
α2 Passou Passou Passou Passou PassouFalhou Falhou Passou Falhou Falhou
β1 Passou Falhou Passou Passou PassouFalhou NA Falhou Falhou Falhou
β2 Passou Falhou Passou Passou PassouPassou NA Falhou Falhou Falhou
β3 Passou Falhou Passou Passou PassouPassou NA Falhou Falhou Passou
β4 Passou Passou Passou Passou PassouFalhou Falhou Falhou Falhou Falhou
β5 Passou Falhou Passou Passou PassouFalhou NA Falhou Falhou Falhou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
APÊNDICE F - Verossimilhança para a Distribuição T deStudent
Quando a suposição de normalidade é fortemente rejeitada, outras alternativaspodem ser usadas para tentar contornar tal problema. Uma delas, proposta por Bollerslev(1987), é supor outra forma para a distribuição dos erros condicionados ao passado, nocaso a distribuição t de student. Nesse caso, a verossimilhança do modelo GARCH é dadapor:
Lt(θ) = f (εt |ψt−1) =Γ(v+1
2
)Γ( v
2
)[(v−2)σ2
t]1/2
(1+
ε2t
(v−2)σ2t
)−( v+12 )
(.6)
e assim,
103
Tabela .9: Temperatura Máxima - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Falhou Falhou Passou Passou Passou
NA NA Falhou Passou Falhou
α1 Passou Falhou Passou Passou PassouPassou NA Passou Falhou Falhou
β1 Falhou Falhou Passou Passou PassouNA NA Falhou Falhou Falhou
β2 Passou Passou Passou Passou PassouFalhou Falhou Falhou Falhou Falhou
β3 Falhou Falhou Passou Passou PassouNA NA Falhou Falhou Falhou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
Tabela .10: Temperatura Mínima - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Falhou Passou Passou Passou Passou
NA Falhou Falhou Passou Falhou
α1 Falhou Falhou Passou Passou PassouNA NA Falhou Falhou Falhou
β1 Falhou Passou Passou Passou PassouNA Falhou Passou Falhou Falhou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
Lc(θ) =n
∏t=q+1
Lt(θ) =n
∏t=q+1
Γ(v+1
2
)Γ( v
2
)[(v−2)σ2
t]1/2
(1+
ε2t
(v−2)σ2t
)−( v+12 ) (.7)
104
Tabela .11: Umidade Relativa - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Passou Passou Passou Passou Passou
Falhou Falhou Falhou Falhou Falhou
α1 Passou Passou Passou Passou PassouPassou Passou Passou Passou Passou
β1 Passou Passou Passou Passou PassouPassou Passou Passou Passou Passou
β2 Passou Passou Passou Passou PassouFalhou Falhou Falhou Falhou Falhou
t NA NA NA NA NANA NA NA NA NA
Fonte: Dados da Estação Automática de Brasília - INMET.
Tabela .12: Velocidade do Vento - Diagnóstico de Heidelberg para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 Passou Passou Passou Passou Passou
Passou Passou Passou Passou Falhou
α1 Passou Falhou Passou Passou PassouFalhou NA Falhou Falhou Falhou
t NA NA Passou NA PassouNA NA Falhou NA Falhou
Fonte: Dados da Estação Automática de Brasília - INMET.
APÊNDICE G - Cálculo do Vetor Score e da MatrizHessiana
A Função Score da Função de Log-Verossimilhança condicional lt(θ) de εt , t =q+1,q+2, . . . ,n é dada por:
∇lt(θ) =∂lt(θ)
∂θ=(
∂lt(θ)∂α0
,∂lt(θ)∂α1
, . . . ,∂lt(θ)∂αq
,∂lt(θ)∂β1
,∂lt(θ)∂β2
, . . . ,∂lt(θ)∂βp
)(.8)
Usando a equação (6.5), a função ou vetor score pode ser reescrita como:
Tabela .13: Evaporação - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 -0,499712532 4,037866113 -0,810546054 0,72800625 -0,625039761α1 1,108566592 -17,42682392 0,486890924 -2,484867164 1,257398563t - - -1,003856014 - 0,254881755
Fonte: Dados da Estação Automática de Brasília - INMET.
105
Tabela .14: Insolação - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 1,905369417 NA -0,99365044 0,884142419 0,837722321α1 1,977725326 NA -0,06054314 -0,23108795 1,811069906α2 0,151501197 NA -4,21308811 -0,57141367 -3,31742013β1 0,346112933 NA -0,96965717 -1,3638849 2,964873397β2 0,924148338 NA 0,079654407 0,570325181 -0,47346685β3 0,256641384 NA 1,971142791 -1,64263443 -1,90809948β4 -1,43651319 NA 1,221946973 1,063260636 -1,97132403β5 -1,54374874 NA 0,135836135 -2,5662906 0,537994955t - - -0,41150866 - 0,691088423
Fonte: Dados da Estação Automática de Brasília - INMET.
Tabela .15: Temperatura Máxima - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 2,287131433 1,581754672 2,246287117 0,515252021 -0,778201923α1 -0,212944365 -0,800563625 -0,338879968 0,134572606 4,015450528β1 0,492684344 -0,23575007 -1,517553283 -0,158381683 -0,074487724β2 -0,731939494 -0,110893511 -2,425187228 1,095059089 -2,219090072β3 -4,741452335 -2,843331071 -2,444864777 0,731971509 0,526078212t - - -1,858360893 - -1,039180949
Fonte: Dados da Estação Automática de Brasília - INMET.
∂lt(θ)∂θ
=∂
∂θ
(−1
2· log(2π)− 1
2· log
(σ
2t)+
ε2t
2σ2t
)(.9)
= −12·− 1
σ2t
[∂σ2
t∂θ
]− ε2
t2· ∂
∂θ
(1
σ2t
)(.10)
mas,
∂
∂θ
(1
σ2t
)=−
(1
σ2t
)· ∂
∂θ
(σ
2t)
(.11)
e então, pode-se reescrever (.10) como:
∂lt(θ)∂θ
=− 12σ2
t·[
∂σ2t
∂θ
](1− ε2
t
σ2t
)(.12)
Tabela .16: Temperatura Mínima - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 3,44462621 NA 0,08689021 0,839009838 -2,16345033α1 0,45885879 NA 0,569161 -0,78872379 -1,58171771β1 -5,19909554 NA -0,30619309 0,593729829 1,837395746t - - -1,84476241 - -2,79774588
Fonte: Dados da Estação Automática de Brasília - INMET.
106
Tabela .17: Umidade Relativa - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 3,125324875 33,82907193 1,715676443 0,31461933 -3,12399725α1 0,738578822 1,005564995 1,672374693 1,44121981 0,175419801β1 0,178840822 1,332583665 -0,41086345 1,07449999 -0,12624451β2 -2,61634416 -46,3829294 -1,4701308 0,20056103 0,277572938t - - 0,683320935 - -0,20312409
Fonte: Dados da Estação Automática de Brasília - INMET.
Tabela .18: Velocidade do Vento - Diagnóstico de Geweke para a amostra da posterioriParâmetros Modelos
M3 M4 M5 M6 M7α0 0,829739308 -0,943293045 0,418860263 0,308902001 2,417992362α1 0,387376352 6,219706767 -0,512973087 -0,483826033 -0,650653996t - - 0,912690093 - 0,047311687
Fonte: Dados da Estação Automática de Brasília - INMET.
Agora, lembrando da definição de volatilidade dada por (5.18), a derivada parcialresultante é resolvida como:
∂σ2t
∂θ= (1,ε2
t−1,ε2t−2, . . . ,ε
2t−q,σ
2t−1,σ
2t−2, . . . ,σ
2t−p)
T +p
∑j=1
∂
∂θ
(σ
2t− j)
(.13)
E, finalmente, basta substituir a equação acima em (.12) para obter a expressãoanalítica do vetor score. Ao final, o vetor score é dado por:
∇lc(θ) =∂lc(θ)
∂θ=
n
∑t=q+1
∂lt(θ)∂θ
(.14)
A matriz Hessiana da função de Log-Verossimilhança condicional é a matriz dassegundas derivadas de lc(θ) em relação ao vetor de parâmetros θ =(α0,α1, . . . ,αq,β1,β2, . . . ,βp)T ,isto é:
H(θ) =∂2lc(θ)∂θ∂θT =
∂
∂θ
(∂lc(θ)∂θT
)(.15)
Para resolvê-la, usando o resultado anterior do vetor score, aparece:
Ht(θ) =∂
∂θ
(− 1
2σ2t·[
∂σ2t
∂θ
](1− ε2
t
σ2t
))=
∂
∂θ
(1
2σ2t
(ε2
t
σ2t−1)[
∂σ2t
∂θ
])mas,
107
∂
∂θ
(1
2σ2t
(ε2
t
σ2t−1))
=12·
(−(
1σ2
t
)2)·(
ε2t
σ2t−1)
+1
2σ2t· ε2
t ·
(−(
1σ2
t
)2)
=12·(
1σ2
t
)2
·(
1− ε2t
σ2t
)·[
∂σ2t
∂θ
]− 1
2·(
ε2t
σ2t
)·[
∂σ2t
∂θ
]=
12·(
1σ2
t
)2
·(
1−2ε2
t
σ2t
)· ∂σ2
t∂θ
e também,
∂
∂θ
(∂σ2
t∂θ
)=
∂2σ2t
∂θ∂θT (.16)
e, portanto, pode-se reescrever a matriz Hessiana em cada tempo t como:
Ht(θ) =1
2σ2t
(ε2
t
σ2t−1)· ∂2σ2
t∂θ∂θT +
12·(
1σ2
t
)2
·(
1−2ε2
t
σ2t
)· ∂σ2
t∂θ· ∂σ2
t∂θT (.17)
Por fim, a expressão da matriz hessiana é encontrada segundo:
H(θ) =n
∑t=q+1
Ht(θ)
