Aplica˘c~ao de Modelos de Volatilidade Estoc astica em ... · Social da Faculdade de Medicina de Ribeir~ao Preto da Universidade de Sao Paulo para a obtenc~ao do t tulo de Mestre

UNIVERSIDADE DE SAO PAULO

FACULDADE DE MEDICINA DE RIBEIRAO PRETO

DEPARTAMENTO DE MEDICINA SOCIAL

Aplicacao de Modelos de Volatilidade Estocastica emDados de Poluicao do Ar de Duas Grandes Cidades:

Cidade do Mexico e Sao Paulo

HENRIQUE CERETTA ZOZOLOTTO

Ribeirao Preto

2010

HENRIQUE CERETTA ZOZOLOTTO

Aplicacao de Modelos de Volatilidade Estocastica emDados de Poluicao do Ar de Duas Grandes Cidades:

Cidade do Mexico e Sao Paulo

Dissertacao apresentada ao Departamento de Medi-cina Social da Faculdade de Medicina de RibeiraoPreto da Universidade de Sao Paulo para a obtencaodo tıtulo de Mestre em Saude na Comunidade.Area de Concentracao: Saude na Comunidade

Orientador: Prof. Dr. Jorge Alberto Achcar

Ribeirao Preto

2010

Autorizo a reproducao e divulgacao total ou parcial deste trabalho, por qualquer meioconvencional ou eletronico, para fins de estudo e pesquisa, desde que citada a fonte.

Ficha Catalografica

Zozolotto, H. C.

Aplicacao de modelos de volatilidade estocastica em dados depoluicao do ar de duas grandes cidades: Cidade do Mexico e Sao Paulo.Ribeirao Preto, 2010.

105 p.:il.; 30cm

Dissertacao de Mestrado apresentada a Faculdade de Medicina deRibeirao Preto - USP. Area de concentracao: Saude na Comunidade.

Orientador: Jorge Alberto Achcar

1. Poluicao do ar. 2. Modelos de volatilidade estocastica. 3. ModelagemBayesiana.

Folha de Aprovacao

Henrique Ceretta Zozolotto

Aplicacao de modelos de volatilidade estocastica em dados de poluicao do ar de duasgrandes cidades: Cidade do Mexico e Sao Paulo

Dissertacao apresentada ao Departamento de MedicinaSocial da Faculdade de Medicina de Ribeirao Preto daUniversidade de Sao Paulo para a obtencao do tıtulo deMestre em Saude na Comunidade.

Area de Concentracao: Saude na Comunidade

Aprovado em: / /

Banca Examinadora

Prof.(a) Dr.(a):

Instituicao: Assinatura:

Prof.(a) Dr.(a):


Prof.(a) Dr.(a):


Dedicatoria

Ao meu pai, Paulo, e as minhas maes,

Maria Amalia, Cristina e Sonia, com amor e

gratidao, pelo empenho e amor dedicado a

mim em todos esses anos de vida.

Agradecimentos

A Deus, primeiramente, que sempre esteve ao meu lado.

A minha namorada Jaina, por ser tao especial e hoje, ser a pessoa mais importante em

minha vida.

Ao meu orientador, Prof. Dr. Jorge Alberto Achcar, pela otima e competente orientacao,

pela paciencia e por ter acreditado em mim quando poucos o teriam feito.

Ao Prof. Dr. Edson Zangiacomi Martinez, por todos os momentos de ajuda e compa-

nheirismo prestados, pela otima orientacao que recebi em meu estagio e principalmente,

pela oportunidade oferecida no Centro de Metodos Quantitativos, em um dos momentos

mais difıceis de minha vida.

Ao meu irmao Tulio, que por mais diferente que seja de mim, e igual a mim. Muito

obrigado por sempre estar ao meu lado.

Aos amigos Davi, Emılio e Roberto e pela ajuda e compreensao tanto nos bons quanto

nos maus momentos.

Ao eterno amigo Gush, pela amizade e pelas aulas de bateria/secoes de terapia.

A todos os meus tios, tias e primos, por serem como pais, maes e irmaos para mim.

E especialmente a minha mae, que mesmo nao estando mais entre nos, consegue ainda

ser o meu maior exemplo de vida.

Resumo

ZOZOLOTTO, H. C. Aplicacao de Modelos de Volatilidade Estocastica em Da-dos de Poluicao do Ar de Duas Grandes Cidades: Cidade do Mexico e SaoPaulo. 2010. 105f. Dissertacao (mestrado) - Faculdade de Medicina de Ribeirao Preto,Universidade de Sao Paulo, Ribeirao Preto, 2010.

Estudos recentes relacionados ao meio ambiente vem ganhando grande destaque em todoo mundo devido ao fato dos nıveis de poluicao e a destruicao das reservas naturais teremaumentado de maneira alarmante nos ultimos anos. As grandes cidades sao as que maissofrem com a poluicao e aqui serao estudados os nıveis de poluicao do ar em duas cidadesem particular, a Cidade do Mexico e Sao Paulo. A Cidade do Mexico apresenta serios pro-blemas com os nıveis de ozonio e Sao Paulo e a cidade brasileira com os maiores problemasrelacionados a poluicao. Entre os diferentes modelos considerados para analisar dados depoluicao do ar, pode-se considerar o uso de modelos de series temporais para modelar asmedias diarias ou semanais de poluicao. Nessa direcao pode-se usar modelos de volatili-dade estocastica. Essa famılia de modelos estatısticos tem sido extensivamente usada paraanalisar series temporais financeiras, porem nao se observa muitas aplicacoes em dadosambientais e de saude. Modelos de volatilidade estocastica bivariados e multivariados,sob a aproximacao Bayesiana, foram considerados para analisar os dados, especialmenteusando metodos MCMC (Monte Carlo em Cadeias de Markov) para obter os sumarios aposteriori de interesse, pois pode-se ter muitas dificuldades usando metodos classicos deinferencia estatıstica.

Palavras - chave: Poluicao do ar, Modelos de volatilidade estocastica, ModelagemBayesiana.

Abstract

ZOZOLOTTO, H. C. Application of Stochastic Volatility Models to Air PollutionData of Two Big Cities: Mexico City and Sao Paulo.. 2010. 105p. Dissertation(master degree) - Faculty of Medicine of Ribeirao Preto, University of Sao Paulo, RibeiraoPreto, 2010.

Recent studies related to environmental has been considered in all world due to increasinglevels of pollution and of natural resources destruction especially, in the last years. Thelargest cities in the world are the ones been mostly affected by pollution and in this workwe consider the analysis of air pollution data of two important cities: Mexico City andSao Paulo. The Mexico City presents serious problems of ozone levels and Sao Paulo isthe Brazilian city with the largest problems related to air pollution. Among the differentmodels which could be used to analyze air pollution data, we consider the use of timeseries modeling to the weekly or daily levels of pollution. In this way, we consider theuse of volatility stochastic models. This family of models has been well explored withfinancial data but not well explored to analyze environmental and health data. Bivariateand multivariate stochastic models under the Bayesian approach were considered to ana-lyze the data, especially using MCMC (Markov Chain Monte Carlo) methods to obtainthe posterior summary of interest, since we usually have big difficulties using standardclassical inference methods.

Keywords: Air pollution, Stochastic volatility models, Bayesian modeling.

Lista de Figuras1 Medias semanais de ozonio para a cidade do Mexico . . . . . . . . . . . . . 18

2 Medias semanais de poluentes para a cidade de Sao Paulo . . . . . . . . . 19

3 Log-retornos centrados em suas medias para as regioes NE, CE e SW de01 de Janeiro de 1990 a 31de Dezembro de 2005. . . . . . . . . . . . . . . . 59

4 Raızes quadradas das volatilidades para as regioes NE, CE e SW quando omodelo 1 e utilizado e quando as regioes sao pareadas em regioes NE e CE(linha superior) e regioes CE e SW (linha inferior). . . . . . . . . . . . . . 61






10 Log-retornos centrados em suas medias para os poluentes de Maio de 1996a Dezembro de 2006. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

11 Raızes quadradas das volatilidades para os poluentes quando o modelo 1 eutilizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

12 Raızes quadradas das volatilidades para os poluentes quando o modelo 2 eutilizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

13 Raızes quadradas das volatilidades para os poluentes quando o modelo 2e utilizado e numero de internacoes por doencas do aparelho circulatorio erespiratorio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

Lista de Tabelas1 Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de cre-

dibilidade (95%) para as quantidades de interesse quando o modelo 1 eutilizado e os pares de regioes NE-CE e CE-SW sao considerados. . . . . . 60

2 Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de cre-dibilidade (95%) para as quantidades de interesse quando o modelo 2 eutilizado e os pares de regioes NE-CE e CE-SW sao considerados. . . . . . 62





7 Deviance Information Criterion para os Modelos 1, 2, 3, 4, 5 e 6 quando ospares de regioes NE-CE e CE-SW sao considerados . . . . . . . . . . . . . 72

8 Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de cre-dibilidade (95%) para as quantidades de interesse quando o modelo 1 eutilizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

9 Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de cre-dibilidade (95%) para as quantidades de interesse quando o modelo 2 eutilizado. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

10 Deviance Information Criterion para os Modelos 1 e 2. . . . . . . . . . . . 80

Sumario

1 Introducao 12

2 Introducao a Metodologia Bayesiana 21

2.1 Formula de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2 Inferencia Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2.1 Estimacao por Intervalo . . . . . . . . . . . . . . . . . . . . . . . . 22

2.3 Distribuicoes a Priori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.3.1 Distribuicoes a Priori Nao-informativas . . . . . . . . . . . . . . . . 24

2.4 Metodos de Simulacao para Obtencao dos Resumos a Posteriori . . . . . . 25

2.4.1 Metodo de Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . 25

2.4.2 Metodo de Monte Carlo em Cadeias de Markov . . . . . . . . . . . 27

2.4.3 Amostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4.4 Algoritmo Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . 30

3 Modelos de Volatilidade Estocastica 32

3.1 Modelos ARCH e GARCH . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2 Modelos de Volatilidade Estocastica Univariados . . . . . . . . . . . . . . . 35

3.2.1 Analise Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3 Modelos de Volatilidade Estocastica Bivariados . . . . . . . . . . . . . . . 39


3.4 Modelos de Volatilidade Estocastica Multivariados . . . . . . . . . . . . . . 50


3.5 Criterios para Selecao de Modelos . . . . . . . . . . . . . . . . . . . . . . . 54

4 Aplicacoes 56

4.1 Analise dos Dados de Poluicao da Cidade do Mexico . . . . . . . . . . . . 56

4.1.1 Descricao dos Dados . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.1.2 Uso dos Modelos e Resultados . . . . . . . . . . . . . . . . . . . . . 59

4.1.3 Selecao do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.1.4 Discussao dos Resultados Obtidos . . . . . . . . . . . . . . . . . . . 72

4.2 Analise dos Dados de Poluicao da Cidade de Sao Paulo . . . . . . . . . . . 73

4.2.1 Descricao dos Dados . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4.2.2 Uso dos Modelos e Resultados . . . . . . . . . . . . . . . . . . . . . 75

4.2.3 Selecao do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.2.4 Discussao dos Resultados Obtidos . . . . . . . . . . . . . . . . . . . 80

5 Conclusoes e Perspectivas Futuras 84

Referencias Bibliograficas 86

A Programas 94

A.1 Cidade do Mexico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.1.1 Modelo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.1.2 Modelo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

A.1.3 Modelo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

A.1.4 Modelo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

A.1.5 Modelo 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

A.1.6 Modelo 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

A.2 Sao Paulo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

A.2.1 Modelo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

A.2.2 Modelo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

12

1 Introducao

Estudos recentes relacionados ao meio ambiente vem ganhando grande destaque em

todo o mundo devido ao fato dos nıveis de poluicao e a destruicao das reservas naturais

terem aumentado de maneira alarmante nos ultimos anos.

As alteracoes climaticas, que estao ocorrendo em nosso planeta, e o aumento da

poluicao podem causar serios problemas a saude da populacao; em especial, as doencas

relacionadas ao aparelho respiratorio e cardiovascular tem sido objetivo de muitas pesquisas

realizadas por profissionais da area medica (GOUVEIA, 2006; BRAGA, 2002).

Os poluentes representam qualquer substancia presente no ar e que, pela sua concen-

tracao, possam torna-lo improprio, nocivo ou ofensivo a saude, causando inconveniente ao

bem estar publico, danos aos materiais, a fauna e a flora ou prejudicial a seguranca, ao uso

e gozo da propriedade e as atividades normais da comunidade (http://www.cetesb.sp.gov.b

r/Ar/ar saude.asp).

Os poluentes sao divididos em duas categorias: primarios e secundarios. Os poluentes

primarios sao aqueles emitidos diretamente de uma fonte, como por exemplo, atraves

de automoveis. Ja os poluentes secundarios sao aqueles que se formam quando alguns

poluentes (primarios) entram em reacao com componentes naturais da atmosfera. Entre

os diversos poluentes primarios existentes pode-se citar os oxidos de nitrogenio (NO e

NO2), dioxido de enxofre (SO2), monoxido de carbono (CO) e os materiais particulados

(MP). E como maior exemplo de poluente secundario, tem-se o ozonio (O3).

Alem desses cinco poluentes citados anteriormente utiliza-se os hidrocarbonetos (HC)

como indicadores da qualidade do ar. Isso se deve ao fato destes serem os mais frequentes

no ar e tambem, os que podem ocasionar maiores efeitos adversos. O nıvel de poluicao

atmosferica e medido atraves da quantidade de substancias desses poluentes presentes no

ar.

Os oxidos de nitrogenio sao gases incolores que resultam de processos em alta tempe-

13

ratura, como a combustao da gasolina, oleo, etc. O NO2 e formado atraves da acao da

luz solar sobre o NO.

O dioxido de enxofre, assim como os oxidos de nitrogenio, e um gas incolor, e e

produzido, principalmente, atraves da combustao de combustıveis que contem enxofre.

O SO2 e liberado para a atmosfera por meio de usinas de energia, metais fundidos e em

maior escala de veıculos.

O monoxido de carbono e resultado da combustao incompleta de gasolina, carvao,

oleo, etc. E um gas muito toxico e que traz grande preocupacao para as cidades, pois os

veıculos automotores sao a principal fonte de emissao deste.

Os materiais particulados sao poluentes (partıculas) que ficam suspensos no ar devido

ao seu pequeno tamanho, como por exemplo poeira e fumaca. Veıculos automotores,

processos industriais e ressuspensao de poeira do solo sao alguns exemplos das fontes de

emissao deste poluente na atmosfera. O MP tambem pode ser resultado de reacoes quımi-

cas de outros poluentes como o dioxido de enxofre e os oxidos de nitrogenio. Os materiais

particulados sao classificados em categorias de acordo com o tamanho das partıculas,

como por exemplo o MP2,5 e o MP10, que representam respectivamente as partıculas cujo

diametro aerodinamico e menor que 2,5µm e 10µm. E no caso deste poluente, quanto

menor o tamanho das partıculas, maior e o dano que pode ser causado a saude.

O ozonio e resultado da mistura de poluentes secundarios formados atraves de reacoes

oxidos de nitrogenio e compostos organicos volateis, na presenca da luz solar. Alem do

ozonio, outras substancias sao formadas atraves dessa mistura. Essas substancias sao

chamadas de oxidantes fotoquımicos. O ozonio e a principal substancia obtida destas

reacoes, por esse fato este e utilizado como parametro indicador da presenca de oxidantes

fotoquımicos na atmosfera.

O aumento da poluicao pode trazer diversos problemas a natureza e a saude da popu-

lacao. Um nıvel elevado dos poluentes pode ocasionar desde irritacoes dos olhos, nariz e

da garganta, bronquite e pneumonia ate doencas respiratorias cronicas, cancer de pulmao,

problemas cardıacos, etc. A exposicao ao dioxido de nitrogenio pode levar ao decrescimo

14

das funcoes respiratorias e aumentar o risco de sintomas respiratorios. As consequencias

da exposicao ao dioxido de enxofre vao desde tosse e dor de garganta ate danos perma-

nentes do sistema respiratorio. Ja a exposicao ao monoxido de carbono pode levar a serios

problemas de coracao e do sistema circulatorio. Caso mulheres gravidas sejam expostas a

um elevado nıvel deste poluente, os bebes poderao nascer com serios problemas de saude.

A exposicao aos materiais particulados pode ocasionar em cancer cardio-respiratorio e

de pulmao, reducao da funcao pulmonar, doencas pulmonares cronicas obstrutivas, etc.

Um alto nıvel de ozonio na regiao perto do solo pode ocasionar em diversos problemas

respiratorios como asma e cancer de pulmao.

Em todo o mundo existem diversos paıses que enfrentam serios problemas relacionados

a poluicao. Recentemente, pode-se observar os elevados nıveis de poluicao que existem

na China. Em 2008, devido aos jogos olımpicos, os orgaos de imprensa deram grande

destaque ao estado alarmante em que se encontravam as cidades chinesas. Outro caso

onde a poluicao se encontrava em um nıvel preocupante, mas com relacao ao ozonio era

a Cidade do Mexico (ACHCAR et al., 2008). No perıodo de estudo que foi de Janeiro

de 1998 a Dezembro de 2004, em 2063 dias o nıvel de ozonio ultrapassou o nıvel padrao

mexicano de ozonio de 0,11ppm.

Quando direcionamos ao Brasil, a cidade de Sao Paulo e a que mais se destaca quanto

aos problemas relacionados a poluicao. Sao Paulo e uma das maiores cidades do mundo

possuindo uma area de 1522,986km2 com mais de 10 milhoes de habitantes (IBGE, 2008).

O seu clima e o subtropical e possui temperatura media anual de 19oC com nıveis de

umidade considerados aceitaveis. Sao Paulo tambem e conhecida como um grande centro

financeiro que hospeda um grande numero de empresas e industrias. E pode-se dizer que

a principal causa dos elevados nıveis de poluicao em Sao Paulo e o excesso de veıculos

automotivos que circulam pelas ruas da cidade.

Diversos estudos foram feitos para avaliar a poluicao e os seus efeitos na saude da

populacao da cidade de Sao Paulo. Farhat et al. (2005) avalariam o efeito da poluicao

no numero de criancas com visita a sala de emergencia e que foram admitidas no hospital

com problemas respiratorios. O perıodo de estudo foi de Agosto de 1996 a Agosto de

15

1997 e os dados mostraram que neste perıodo os padroes diarios foram ultrapassados tres

vezes para o MP10 e o NO2. O padrao de 1 hora para o O3 foi ultrapassado 14 vezes e o

padrao de 8 horas do CO, 4 vezes. Somente os nıveis de SO2 se mantiveram abaixo do

padrao estabelecido, durante o perıodo de estudo. Os resultados apresentados indicam

a existencia de correlacoes moderadas e altas entre os poluentes, exceto para o ozonio.

Alem disso, indıcios de existencia de associacao entre a poluicao do ar e a morbidade por

doenca respiratoria em criancas foram obtidos atraves da analise realizada.

Martins et al. (2002) investigaram os efeitos causados pela poluicao atmosferica na

morbidade por pneumonia e por gripe em idosos entre 1996 e 1998. Assim como para

o estudo feito em criancas, os resultados indicaram que a poluicao atmosferica pode re-

sultar em problemas a saude de idosos. Como outros exemplos de artigos onde foram

estudadas as relacoes entre os nıveis dos poluentes e os ındices de morbi-mortalidade em

cidades brasileiras temos, Cendon et al. (2006) estudaram a relacao entre a poluicao e

a morbidade por infarte do miocardio na cidade de Sao Paulo. Freitas et al. (2004),

Conceicao et al. (2001), Gouveia et al. (2006) e Martins et al. (2006) avaliaram a relacao

entre alguns poluentes atmosfericos com internacoes e obitos na cidade de Sao Paulo. Ja

Marcilio e Gouveia (2007) avaliaram a relacao entre alguns poluentes atmosfericos com

internacoes e obitos em sete cidades brasileiras. Moura et al. (2008) pesquisaram sobre

a relacao entre a qualidade do ar e transtornos respiratorios agudos em criancas no Rio

de Janeiro, Nascimento et al. (2006) investigaram a relacao entre a poluicao e o numero

de internacoes por pneumonia em criancas no interior de Sao Paulo e Medeiros e Gouveia

(2005) estudaram a relacao entre o baixo peso ao nascer e a poluicao do ar na cidade

de Sao Paulo. Em todos esses casos foram utilizadas metodologias classicas de analise

estatıstica, principalmente modelos lineares generalizados e coeficientes de correlacao.

Apesar da grande quantidade de estudos que buscam relacionar os nıveis de poluicao

com o numero de internacoes e obitos nas cidades brasileiras e ate em outros paıses,

pouco foi feito ate entao com relacao ao estudo da volatilidade para essas series de dados

relacionadas a poluicao do ar (ver ACHCAR et al., 2008a; ACHCAR et al., 2008b).

A volatilidade e um termo muito empregado na area financeira. Mas antes de definir

16

volatilidade e preciso introduzir, neste contexto de meio ambiente e saude, um conceito

bastante utilizado na area de series temporais financeiras.

O conceito mais comum quando se trabalha com series temporais financeiras e o con-

ceito de retorno. Mas primeiro e necessario entender o conceito de ativo. Os ativos sao

todos os bens, direitos e valores que uma empresa possui. Os ativos podem ser imoveis,

dinheiro, mercadorias, entre outros. E o retorno e definido como a variacao relativa

de precos de um determinado ativo entre dois instantes de tempo (MORETTIN, 2008).

Isto e, seja P (t) o preco de um ativo em um instante de tempo t, o retorno e definido

como P (t)−P (t−1)P (t−1) . Porem, na pratica, e mais comum se trabalhar, por questoes estatısti-

cas, com o logaritmo do retorno. Nesta dissertacao serao analisados dados de poluicao

do ar, e neste contexto podemos definir o retorno como sendo a variacao relativa entre as

concentracoes de um determinado poluente, entre dois instantes de tempo.

Com isso, pode-se definir a volatilidade como a variancia condicional de um retorno,

isto e a variacao de um retorno em um dado instante de tempo, dependente dos retornos e

outras informacoes passadas. Trazendo este conceito para o tipo de dados que sera usado

neste trabalho, a volatilidade e definida como a variacao de um retorno de um poluente

em questao em um dado instante de tempo, onde essa variacao depende dos valores dos

retornos e de outras informacoes passadas. Os ganhos que se pode ter ao estudar a

volatilidade neste tipo de estudo sao varios. Atraves dos valores estimados da volatilidade

pode-se saber se estao ocorrendo grandes variacoes dos poluentes em determinadas epocas

do ano. Tambem podemos estimar a probabilidade de que a concentracao de um poluente

venha a cair ou subir em um perıodo de tempo especifico. A aplicacao de metodologias

Bayesianas para a analise de dados de volatilidade vem ganhando grande destaque, tanto

na area financeira como na saude, principalmente com o uso dos modelos de volatilidade

estocastica.

O uso de metodologias Bayesianas para analisar dados relacionados a poluicao do

ar e seus efeitos na saude da populacao tem se tornado uma pratica mais frequente em

estudos realizados em diferentes cidades do mundo. Berhane e Molitor (2008) estudaram

os efeitos da poluicao do ar nas funcoes pulmonares em criancas do sul da California. Os

17

autores destacam neste artigo a dificuldade de se trabalhar com as taxas de crescimento

das funcoes pulmonares, devido ao fato destas variarem de indivıduo para indivıduo, por

exemplo quando as criancas passam da pre-adolescencia a puberdade. Alem disso os

autores comentam que ainda nao existem metodologias bem estabelecidas para relacionar

as curvas de crescimento das funcoes pulmonares com os nıveis de poluicao. Os autores

fazem uso de tecnicas Bayesianas para obter estimativas a posteriori das distribuicoes de

caracterısticas importantes das curvas de crescimento. Lin et al. (2008) estudaram a

relacao entre a concentracao de ozonio e as internacoes respiratorias de criancas no estado

de Nova York em um perıodo de onze anos. Era de interesse dos pesquisadores avaliar os

efeitos da exposicao ao ozonio na saude das criancas em 11 regioes do estado e de maneira

geral. Para isso, estes utilizaram modelos Bayesianos hierarquicos de dois estagios. Os

resultados indicaram a existencia de associacao em 5 das 11 regioes estudadas. Welty et

al. (2009) utilizaram modelos Bayesianos de defasagem distribuıda para estudar os efeitos

da concentracao de material particulados na mortalidade diaria na cidade de Chicago

entre 1987 e 2000.

Alem destes trabalhos acima mencionados, diversos outros artigos envolvendo o uso

de metodologias Bayesianas com aplicacao em estudos epidemiologicos envolvendo o meio

ambiente vem sendo publicados (DOMINIC et al., 2008; BURSTYN et al., 2008).

Como motivacao e exemplo inicial vamos considerar os dados da concentracao semanal

de ozonio (medias diarias) na Cidade do Mexico no perıodo compreendido entre Janeiro

de 1990 a Dezembro de 2005 medidos pelo Instituto de Ecologia da Cidade do Mexico

(www.sma.df.gob.mx/simat/). Para esses dados foram consideradas medidas de 3 regioes

da Cidade do Mexico: Nordeste (NE), Central (CE) e Sudoeste (SW).

18

(a) NE (b) CE

(c) SW

Figura 1: Medias semanais de ozonio para a cidade do Mexico

Observar (ver Figura 1) que para as regioes da Cidade do Mexico, ha uma forma

decrescente para o ozonio.

Outro exemplo a ser utilizado sao os dados de poluicao do ar da cidade de Sao

Paulo, onde as medidas diarias de poluentes (SO2, NO2, PM10, O3 e CO) foram cole-

tadas pelo Instituto de Astronomia, Geofısica e Ciencias Atmosfericas da USP (http://ww

w.iag.usp.br/) no perıodo compreendido entre Maio de 1996 a Dezembro de 2006 (ver

Figura 2).

19

(a) CO (b) NO2

(c) PM10 (d) SO2

(e) O3

Figura 2: Medias semanais de poluentes para a cidade de Sao Paulo

Observamos (ver Figura 2) que alguns poluentes (CO, SO2) tem um comportamento

decrescente mais evidente, mas os outros poluentes aparentemente nao tem esse compor-

tamento. Isso pode ser uma informacao relevante que requer uma analise mais detalhada

dos dados.

Para analisar os dados de poluicao das figuras 1 e 2 podemos usar varias metodologias

estatısticas, em especial sob o enfoque Bayesiano. Em especial vamos considerar modelos

20

de volatilidade estocastica para as series temporais dos poluentes.

Em resumo, os objetivos deste trabalho sao:

Aplicar e comparar os diferentes modelos de volatilidade estocastica para modelar

os dados de poluicao, sob o enfoque Bayesiano, a serem apresentados;

Analisar dados de poluicao para duas grandes cidades: Cidade do Mexico e Sao

Paulo;

Com base nos resultados obtidos dos modelos aplicados estudar como esta o com-

portamento da volatilidade para os poluentes estudados em cada cidade;

Por meio destes modelos verificar se as medidas que estao sendo tomadas pelas

autoridades ambientais das cidades em estudo estao surtindo um efeito positivo

para reduzir os nıveis de poluicao;

Estudar a possıvel relacao existente entre maior volatilidade de alguns poluentes

com o numero de internacoes/obitos associados a diferentes doencas.

21

2 Introducao a Metodologia Bayesi-

ana

Os metodos Bayesianos tem sido amplamente utilizados como alternativa aos metodos

classicos de inferencia estatıstica. Um ponto importante da inferencia Bayesiana e que

esta permite a incorporacao de uma informacao a priori de um especialista. Outra carac-

terıstica que merece destaque dos metodos Bayesianos e que os parametros de interesse

sao considerados como sendo quantidades aleatorias.

A fundamentacao da teoria de inferencia Bayesiana e baseada na formula de Bayes.

Esta combina a informacao dos dados com a informacao a priori, obtida atraves de um

especialista, e entao se obtem a distribuicao a posteriori, onde sera realizado todo o

processo inferencial.

2.1 Formula de Bayes

Seja A1, A2, ..., Ak uma sequencia de eventos mutuamente exclusivos e exaustivos for-

mando uma particao do espaco amostral Ω, isto e,k⋃j=1

Aj = Ω e Ai∩ Aj = φ (conjunto

vazio) para i 6= j tal que P

(k⋃j=1

Aj

)=

k∑j=1

P (Aj) = 1.

Entao para qualquer outro evento B (B ⊂ Ω), temos,

P (Ai | B) =P (B | Ai)P (Ai)k∑j=1

P (B | Aj)P (Aj)

(2.1)

para todo i variando de 1 ate k.

Podemos interpretar a formula de Bayes (2.1) da seguinte forma: antes do conheci-

mento de qualquer informacao sobre o evento Ai, atribuımos uma probabilidade a priori

22

para Ai, dada por P (Ai); essa probabilidade e atualizada a partir da ocorrencia do evento

B. Essa probabilidade atualizada, ou probabilidade condicional do evento Ai dado a ocor-

rencia do evento B, ou seja, P (Ai | B) e dada pela formula de Bayes (2.1).

Assumindo agora que temos um vetor de dados y = (y1, . . . , yn)′

e quantidades desco-

nhecidas θ representando os parametros de uma distribuicao de probabilidade associada

com a variavel aleatoria Yi com valores observados yi, i = 1, . . . , n.

Considerando uma amostra aleatoria y = (y1, . . . , yn)′, isto e, os dados sao indepen-

dentes e identicamente distribuıdos com uma distribuicao conjunta dada por f (y | θ),

tambem definida como funcao de verossimilhanca para θ quando os dados foram obser-

vados e uma distribuicao a priori para θ, dada por π (θ), assumindo os valores discretos

θ1, . . . , θk, temos de (2.1), a distribuicao a posteriori para θi dado y,

π (θi | y) =f (y | θi)π (θi)k∑j=1

f (y | θj) π (θj)

(2.2)

Observar que o parametro θ tambem e considerado como uma quantidade aleatoria

sob o enfoque Bayesiano.

Supondo agora que o parametro θ assume valores contınuos num dado intervalo, pode-

mos escrever (2.2) por

π (θ | y) =f (y | θ) π (θ)∫f (y | θ) π (θ) dθ

(2.3)

em que a integral no denominador de (1.4) e definida no intervalo de variacao de θ.

2.2 Inferencia Bayesiana

2.2.1 Estimacao por Intervalo

Na inferencia classica, assim como na inferencia Bayesiana, e de interesse do pes-

quisador estimar um intervalo de valores no qual se espera que esteja o parametro de

23

interesse, com um certo grau de confianca, ao inves de utilizar apenas o valor estimado

do parametro. Na inferencia classica esse intervalo e chamado de intervalo de confianca.

Considerando um grau de confianca, por exemplo igual a 95%, a interpretacao de um

intervalo de confianca para um parametro θ e dada por: se retirassemos da populacao um

numero grande de amostras de tamanho n, 95% destas amostras iriam gerar intervalos de

confianca que contem o verdadeiro valor de θ.

Diferente da inferencia classica, a construcao e interpretacao de um intervalo para

um parametro de interesse, em inferencia Bayesiana, e feita de maneira muito mais facil.

Seja θ um parametro unidimensional e assumir que a distribuicao a posteriori para θ seja

unimodal. Um estimador por intervalo Bayesiano com probabilidade (1− α) sera dado

por (θ∗, θ∗) para ∫ θ∗

−∞π (θ | y) dθ =

α

2

e ∫ ∞θ∗

π (θ | y) dθ =α

2

O intervalo (θ∗, θ∗) e chamado um intervalo de credibilidade para θ com probabilidade

(1− α). E sua interpretacao e: a probabilidade do verdadeiro valor de θ estar no intervalo

(θ∗, θ∗) e igual a (1− α) .

2.3 Distribuicoes a Priori

Na inferencia Bayesiana, os parametros sao considerados como quantidades aleatorias

e a estes sao atribuıdas distribuicoes a priori. O uso de informacoes a priori sao a grande

vantagem dos metodos Bayesianos quando comparados aos metodos classicos, pois estas

permitem a incorporacao de informacoes de especialistas em adicao a informacao contida

nos dados. A elicitacao da distribuicao a priori para os parametros de interesse pode ser

feita de varias formas:

(i) Podemos assumir distribuicoes a priori definidas no domınio de variacao do pa-

rametro de interesse. Como caso particular, poderıamos considerar uma distribuicao a

24

priori Beta que e definida no intervalo (0, 1) para proporcoes que tambem sao definidas no

intervalo (0, 1) ou considerar uma priori normal para parametros definidos em toda reta;

(ii) Podemos construir uma priori baseada em informacoes de um ou mais especialistas;

(iii) Podemos considerar metodos estruturais de elicitacao de distribuicoes a priori

(PAULINO et al., 2002);

(iv) Podemos considerar distribuicoes a priori nao-informativas quando temos total

ignorancia sobre parametros de interesse;

(v) Podemos usar metodos Bayesianos empıricos em dados ou experimentos previos

para construir a priori de interesse.

Dentre as formas apresentadas acima, uma que merece um destaque e o item (iv), o

uso de priori nao-informativas.

2.3.1 Distribuicoes a Priori Nao-informativas

Muitas vezes nao temos informacoes que podem ser obtidas de especialistas, ou entao

as informacoes que obtemos nao sao dignas de muita credibilidade. Nesses casos, uma

alternativa e fazer uso de districoes a priori nao-informativas. Os objetivos de se usar

distribuicoes a priori nao-informativas podem ser:

(i) Deduzir crencas a posteriori para quem parte de um conhecimento escasso, isto e,

quando os dados fornecem grande parte da informacao sobre o parametro (ignorancia a

priori).

(ii) Permitir a comparacao com os resultados obtidos da inferencia classica que so usa

a informacao amostral.

(iii) Averiguar a influencia de uma priori subjetiva quando comparada com os resul-

tados obtidos usando uma priori nao-informativa.

25

2.4 Metodos de Simulacao para Obtencao dos Re-

sumos a Posteriori

Para obtencao dos sumarios a posteriori de interesse, geralmente precisamos resolver

integrais Bayesianas que nao apresentam solucao analıtica. Para isso varias alternativas

foram introduzidas na literatura para resolver essas integrais Bayesianas, como por exem-

plo a aproximacao de Laplace. Porem, muitas vezes o uso da aproximacao de Laplace

para resolucao das integrais Bayesianas pode nao ser o mais indicado, pois dependendo

da complexidade das integrais, o tempo gasto nessas resolucoes pode ser muito grande.

Com o desenvolvimento computacional, os metodos de simulacao puderam ser utilizados

para a obtencao dos resumos de interesse. Um metodo bastante utilizado e o de Monte

Carlo.

2.4.1 Metodo de Monte Carlo

Supor que estamos interessados em aproximar uma integral na forma,

E [g (θ) | y] =

∫g (θ) π (θ | y) dθ (2.4)

em que y e θ podem ser vetores.

Pelo metodo de Monte Carlo ordinario, simular uma amostra θ1, . . . , θk da distribuicao

a posteriori π (θ | y). Assim, (2.4) e aproximado por,

E [g (θ) | y] =1

n

n∑i=1

g (θi) (2.5)

Observar que pela lei forte dos grandes numeros, E [g (θ) | y] converge quase certa-

mente para E [g (θ) | y].

Intervalos de credibilidade para θ podem ser obtidos usando o metodo de Monte Carlo

26

ordinario. Ordenar a amostra simulada de π (θ | y) : θ(1) < θ(2) < . . . < θ(n). Um intervalo

de credibilidade 100(1− α)% para θ e dado por,

Rc (α) =[θ (1−α)

2

; θ (1+α)2

](2.6)

cujos extremos definem quantis de probabilidade a posteriori (1−α)2

e (1+α)2

de θ.

Isto e, P[θ ≤ θ (1−α)

2

| y]

= (1−α)2

e P[θ ≤ θ (1+α)

2

| y]

= 1− (1−α)2

= (1+α)2

.

Quando nao conseguimos gerar valores diretamente da posteriori π (θ | y), utilizamos

uma extensao do metodo de Monte Carlo, que e o metodo de Monte Carlo via funcao de

importancia.

Seja p (θ) uma densidade da qual seja facil simular amostras e que aproxime a dis-

tribuicao π (θ | y).

Assim, podemos escrever (2.4) na forma,

∫g (θ)π (θ | y) dθ =

∫g (θ) f (y | θ) π (θ) dθ∫f (y | θ) π (θ) dθ

(2.7)

=

∫g (θ) f(y|θ)π(θ)

p(θ)p (θ) dθ∫ f(y|θ)π(θ)

p(θ)p (θ) dθ

=

∫g (θ)w (θ) p (θ) dθ∫w (θ) p (θ) dθ

em que w (θ) = f(y|θ)π(θ)p(θ)

, f (y | θ) e a funcao de verossimilhanca para θ e π (θ) e a

distribuicao a priori para θ.

Obtendo uma amostra θ1, . . . , θn de p (θ), encontramos a aproximacao de Monte Carlo

para E [g (θ) | y] dada por,

E [g (θ) | y] =1∑ni=1wi

n∑i=1

wig (θi) (2.8)

em que

wi =f (y | θi) π (θi)

p (θi)

27

2.4.2 Metodo de Monte Carlo em Cadeias de Markov

Supor que temos interesse em gerar uma amostra de uma distribuicao a posteriori

π (θ | y), θ ∈ Rk, mas nao podemos fazer isso diretamente. Entretanto, supor que pode-

mos construir uma cadeia de Markov com espaco de estados no espaco parametrico Θ

(conjunto de todos valores possıveis de θ) que e simples para simular e cuja distribuicao

de equilıbrio seja dada por π (θ | y). Se temos muitas simulacoes dessa cadeia, os valores

simulados da cadeia podem ser usados como uma base para sumarizar caracterısticas da

posteriori π (θ | y).

Resultado: (BESAG, 1994). Se a distribuicao conjunta a posteriori π (θ | y) for po-

sitiva em Θ1× Θ2×· · ·× Θk, com Θi sendo suporte para a distribuicao de θi, i = 1, . . . , k,

entao a distribuicao a posteriori π (θ | y) e unicamente determinada pelas distribuicoes

condicionais completas π(θi| y,θ(i)

)para i = 1, . . . , k em que θ = (θ1, . . . , θk) e θ(i) e o

vetor de todos os componentes de θ exceto θi, isto e, θ(i) = (θ1, . . . , θi−1, θi+1, . . . , θk).

Sob algumas condicoes de regularidade e facil ver que os resultados simulados da cadeia

com distribuicao de equilıbrio π (θ | y) podem ser supostos com uma amostra aleatoria

de π (θ | y).

Se θ(1),θ(2), . . . ,θ(t), . . . e uma realizacao de uma cadeia, temos,

θ(t)D→ θ ∼π (θ | y) (2.9)

em que o sımbolo “D” significa convergencia em distribuicao.

Da mesma forma para estimar o valor esperado de g (θ) com respeito a π (θ | y), isto

e,

E [g (θ | y)] =

∫g (θ) π (θ | y) dθ (2.10)

observamos que

1

t

t∑i=1

g(θ(i))q.c.→ E [g (θ | y)] (2.11)

(“q.c.”: convergencia quase certa).

28

Na pratica, θ(i) pode estar correlacionado, mas poderıamos considerar espacos ade-

quados entre os θ(i) gerados para garantir uma amostra aleatoria de π (θ | y).

E necessario gerarmos cadeias de Markov cuja distribuicao de equilıbrio seja exata-

mente π (θ | y) , para isso podemos utilizar os algoritmos Metropolis-Hastings e o amos-

trador de Gibbs.

2.4.3 Amostrador de Gibbs

Supor que estamos interessados em obter inferencias da distribuicao a posteriori con-

junta, π (θ | y), θ = (θ1, . . . , θk). Para isso simulamos quantidades aleatorias de dis-

tribuicoes condicionais completas π(θi| y,θ(i)

)que produzem uma cadeia de Markov.

Observar que π(θi| y,θ(i)

)sao facilmente identificadas como funcoes de θi por ins-

pecao da forma de π (θ | y) a distribuicao a posteriori para θ dado y. (GAMERMAM,

1997).

Supor que atribuımos um conjunto arbitrario de valores iniciais θ(0)1 , θ

(0)2 , . . . , θ

(0)k para

o vetor de parametros θ.

Daı, escrevemos o algoritmo:

(i) Gerar θ(1)1 de π

(θ1 | y, θ(0)2 , . . . , θ

(0)k

); (2.12)

(ii) Gerar θ(1)2 de π

(θ2 | y, θ(1)1 , θ

(0)3 , . . . , θ

(0)k

);

(iii) Gerar θ(1)3 de π

(θ3 | y, θ(1)1 , θ

(1)2 , θ

(0)4 , . . . , θ

(0)k

);

.

.

.

(k) Gerar θ(1)k de π

(θk | y, θ(1)1 , θ

(1)2 , . . . , θ

(1)k−1

)

Entao, substituir os valores iniciais com uma nova realizacao θ(1) =(θ(1)1 , θ

(1)2 , . . . ,

θ(1)k

)de θ e repetir o processo acima.

29

Para t suficientemente grande, observar que o valor θ(t)1 , θ

(t)2 , . . . , θ

(t)k converge para

um valor da quantidade aleatoria com distribuicao π (θ | y) (GEMAN e GEMAN, 1984).

Alem disso, θ(t)j pode ser considerado como uma observacao simulada da distribuicao a

posteriori marginal π (θj| y), j = 1, 2, . . . , k.

Replicando o processo acimaB vezes obtemosB vetores (θ(g)1 , θ

(g)2 , . . . , θ

(g)k ); g = 1, 2, . . . ,

B.

Da convergencia do amostrador de Gibbs, qualquer caracterıstica da densidade a pos-

teriori marginal π (θj| y) pode ser obtida.

Em particular, se π(θj| θ(j),y

)e dada em forma fechada, entao

π (θj| y) =1

B

B∑g=1

π(θj| θg(j),y

)(2.13)

em que j = 1, . . . , k.

Para verificar a convergencia do algoritmo, podemos considerar varias tecnicas. Gelfand

e Smith (1990) sugerem o uso de tecnicas graficas; dessa forma considerar varias cadeias

paralelas geradas a partir de valores iniciais diferentes. Apos um grande numero de

interacoes em cada cadeia, comparar os histogramas para cada componente θj de θ. His-

togramas similares, indicam convergencia da cadeia.

Geweke (1992) sugere metodos graficos baseados em series temporais das amostras

selecionadas.

Uma tecnica para monitorar a convergencia do algoritmo e proposta por Gelman e

Rubin (1992) baseada na analise de variancia.

Nota: Na geracao de amostras de Gibbs devemos considerar as l primeiras itera-

coes como perıodo de aquecimento (“burn-in-samples”), que devem ser descartadas para

eliminar o efeito de valores iniciais.

30

2.4.4 Algoritmo Metropolis-Hastings

O amostrador de Gibbs nos permite gerar amostras da distribuicao a posteriori con-

junta desde que as distribuicoes condicionais completas possuam formas fechadas. O

algoritmo Metropolis-Hastings nos permite gerar amostras da distribuicao a posteriori

conjunta com distribuicoes condicionais completas possuindo ou nao forma fechada.

Quando usamos uma priori conjugada, as distribuicoes condicionais em geral sao de

forma conhecida como, por exemplo, as distribuicoes normal, gama, Poisson, beta, etc, e

a simulacao de amostras dessas distribuicoes sao disponıveis em qualquer software.

Quando as distribuicoes condicionais nao sao facilmente identificadas, devemos usar o

algoritmo de Metropolis-Hastings ou metodos de amostragem por importancia.

Supor que desejamos simular amostras de uma densidade nao-regular π(θi | θ(i),y

),

ou simplesmente π(θi | θ(i)

), em que θ(i) = (θ1, . . . , θi−1, θi+1, . . . , θk).

Definir o nucleo de transicao q (θ, β) da distribuicao p (θ) que representa π(θi | θ(i)

)e

que transforma θ em β. Se θ e uma variavel real com amplitude em toda reta R, podemos

construir q tal que β = θ+σz, com Z ∼ N (0, σ2), em que σ2 reflete a variancia condicional

de θ em p (θ).

Se θ e limitado com amplitude (a, b) usar uma transformacao que leva (a, b) em

(−∞,∞) e daı usar o nucleo de transicao q e aplicar o algoritmo de Metropolis para

a densidade da variavel transformada. O algoritmo de Metropolis e dado por:

(i) Iniciar com um valor θ(0) e indicador de estagio, j = 0;

(ii) Gerar um ponto β do nucleo de transicao q(θ(j), β

);

(iii) Atualizar θ(j) por θ(j+1) = β com probabilidade,

p = min

1,p (β) q

[θ(j), β

]p[θ(j)]q[β, θ(j)

] (2.14)

Ficar com θ(j) com probabilidade 1− p;

31

(iv) Repetir os estagios (ii) e (iii) ate conseguir uma distribuicao estacionaria. Ob-

servar que:

(a) O algoritmo de Metropolis Hastings e especificado pela densidade candidata para

geracao q (x, y) ;

(b) Se um valor candidato e rejeitado, o valor atual e considerado na proxima etapa;

(c) O calculo de p em (2.14) nao depende da constante normalizadora;

(d) Se a densidade candidata para geracao das amostras e simetrica, isto e, q (x, y) =

q (y, x), a probabilidade de movimento se reduz a p(β)

p[θ(j)]. Assim, se p (β) > p

[θ(j)], a

cadeia se move para β; em caso contrario, ela se move para β com probabilidade p(β)

p[θ(j)].

Em outras palavras: um salto na direcao“ascendente” e sempre aceito; um salto na direcao

“descendente” e aceito com uma dada probabilidade.

32

3 Modelos de Volatilidade Estocas-

tica

Uma serie temporal e definida como uma sequencia de observacoes medidas ao longo

do tempo. A analise de uma serie temporal tem como objetivos entender a estrutura que

gerou a serie e prever valores futuros desta. Os modelos de series temporais sao, principal-

mente, utilizados nas areas financeiras e de economia. Porem, como citado anteriormente,

ja temos na literatura diversos casos onde foram analisadas series de dados de poluentes

buscando relacionar estes com algumas doencas. Atualmente, devido a crise ambiental,

muitos estudos tem sido feitos para melhorar a condicao do meio ambiente, e o uso de

modelos de series temporais tem se mostrado de grande importancia nessa area.

Entre os diferentes modelos considerados para analisar dados de poluicao do ar, pode-

mos considerar o uso de modelos de series temporais para modelar as medias diarias

ou semanais de poluicao (LOOMIS et al., 1996). Nessa direcao podemos usar modelos

de volatilidade estocastica (MVE) (GHYSELS et al., 1996; KIM et al., 1998; MEYER

e YU, 2000). Essa famılia de modelos estatısticos tem sido extensivamente usada para

analisar series temporais financeiras (DANIELSSON,1994; YU, 2002), como uma melhor

alternativa quando comparada com modelos usuais de series temporais financeiras de

tipo ARCH (autoregressive conditional heteroscedastic) introduzidos por Engle (1982) e

os modelos GARCH (generalized autoregressive conditional heteroscedastic) introduzidos

por Bollerslev (1986). O uso de modelos de volatilidade estocastica tem muitas vantagens

para analisar series temporais, pois esses modelos assumem dois processos para mode-

lar as series:um processo para modelar as observacoes e um processo para modelar as

volatilidades latentes. Inferencia Bayesiana usando metodos de Monte Carlo em Cadeias

de Markov (MCMC) (GELFAND e SMITH, 1990; SMITH e ROBERTS, 1993) tem sido

considerada para analisar modelos de volatilidade estocastica, pois podemos ter muitas di-

ficuldades usando metodos classicos de inferencia estatıstica, como alta dimensionalidade,

33

funcao de verossimilhanca com formas nao fechadas e grande custo computacional.

3.1 Modelos ARCH e GARCH

Muitas vezes, pode ser de interesse do pesquisador estudar a volatilidade da serie de

dados ao inves de simplesmente investigar o comportamento da serie ao longo do tempo.

A volatilidade pode ser definida como a variancia condicional de uma observacao, isto e a

variacao de uma observacao em um dado instante de tempo, dependente das observacoes

passadas. Uma caracterıstica importante da volatilidade e que esta nao e diretamente

observavel. Entre os diversos modelos utilizados para estudar a volatilidade de uma serie

de dados, os mais usuais sao os modelos ARCH e GARCH.

Usualmente, devido a algumas propriedades estatısticas como por exemplo a esta-

cionariedade, trabalha-se com o log-retorno da serie de dados ao inves desta em sua escala

original. Seja z(t), t = 1, ..., N uma serie temporal. O seu log-retorno pode ser definido

como:

y(t) = logz(t)

z(t− 1)

Os modelos de volatilidade, diferentemente dos modelos tradicionais de series tem-

porais, nao consideram uma variancia constante para a serie. Os modelos ARCH intro-

duzidos por Engle (1982) permitem que a variancia condicional (volatilidade) varie com

tempo como uma funcao dos erros passados. Nos modelos ARCH a volatilidade depende

dos retornos passados por meio de uma funcao quadratica (MORETTIN, 2008). Um

modelo ARCH de ordem p pode ser definido como:

y(t) = ε(t)√σ(t)

onde

σ(t) = α0 + α1y2(t−1) + α2y

2(t−2) + ...+ αpy

2(t−p) = α0 +

p∑l=1

αly2(t−l),

34

e onde y(t) sao os log-retornos de z(t) centrados em suas medias, e ε(t)iid∼ N (0, σ2

ε) denota

uma distribuicao Normal com media 0 e variancia σ2ε , porem na pratica considera-se σ2

ε =

1, ou ε(t) ∼ tv, onde tv denota uma distribuicao t de Student com v graus de liberdade.

Uma desvantagem desses modelos e que estes consideram os retornos ao quadrado, isto e os

retornos negativos e positivos sao tratados da mesma forma. Isso pode ser um problema,

pois a volatilidade nao e a mesma para retornos positivos e negativos. Outra desvantagem

desses modelos e que pelo fato de trabalhar com os retornos ao quadrado, nas situacoes

onde for observado algum retorno grande, corre-se o risco de fazer super-previsoes.

Outros modelos bastante utilizados para modelar volatilidades sao os modelos GARCH

propostos por Bollerslev (1986). Estes modelos sao uma generalizacao dos modelos ARCH

e a sua diferenca com relacao a estes e que os modelos GARCH incorporam os valores

passados da volatilidade, alem dos valores passados do processo ao quadrado. Estes

modelos representam uma tentativa de expressar, com mais parcimonia, a dependencia

temporal da variancia condicional. Um modelo GARCH (p, q) pode ser definido como:

y(t) = ε(t)√σ(t)

onde

σ(t) = α0 +

p∑i=1

αiy2(t− 1) +

q∑j=1

βjσ(t− j)

e onde y(t) sao variaveis aleatorias i.i.d. com medias iguais a zero, p > 0, q ≥ 0, α0 > 0,

αi ≥ 0, i = 1, ..., p, βj ≥ 0, j = 1, ..., q,p∑i=1

αi +q∑j=1

βj < 1.

Alem dos modelos GARCH existem diversas outras generalizacoes dos modelos ARCH,

que eliminam alguns dos problemas citados anteriormente. Porem, o objetivo desta disser-

tacao sera a aplicacao dos modelos de volatilidade estocastica com o auxılio de inferencia

Bayesiana para o estudo da volatilidade.

35

3.2 Modelos de Volatilidade Estocastica Univariados

Diferentes classes de modelos de volatilidade estocastica sao introduzidos na litera-

tura (YU e MEYER, 2005). Os modelos de volatilidade estocastica (MVE) constituem

boas alternativas a analise de series de tempo quando comparados aos modelos ARCH e

GARCH. Diferente dos modelos ARCH, que supoem que a variancia condicional depende

dos retornos passados, os modelos de volatilidade estocastica supoem que a volatilidade

depende de seus valores passados, porem e independente dos retornos passados.

Considerando uma serie de dados z(t) define-se y(t) como sendo os log-retornos de

z(t) centrados em suas medias. Na presenca de heterocedasticidade vamos assumir que o

modelo para a serie temporal y(t) seja dado por,

y(t) = σ(t)ε(t) (3.1)

em que ε(t) e um ruıdo suposto independente e identicamente distribuıdo com distribuicao

Normal N (0, σ2ε). Denotar por P (t) = (y(t− 1), y(t− 2), ...) os valores passados de y(t).

A volatilidade da serie y(t) e a variancia condicional σ2(t) = E (y2(t) | P (t)). Tambem

assumir que o desvio-padrao e dado pelo modelo:

σ(t) = exp

(h(t)

2

)(3.2)

sendo h(t) uma variavel latente definida por um modelo auto regressivo AR(1) dado por,

h(t) = µ+ φ (h(t− 1)− µ) + η(t) (3.3)

para t = 2, 3, ..., N . Onde φ e o parametro auto regressivo do modelo. Vamos assumir

que h(1) seja uma variavel aleatoria com uma distribuicao conhecida e η(t) e um ruıdo

associado as variaveis latentes com uma distribuicao normal N(0, σ2

η

). Se |φ| < 1, a media

e a variancia nao-condicional de h(t) sao respectivamente E(h(t)) = µ e var(h(t)) =σ2η

1+φρ1,

sendo ρ1 o coeficiente de autocorrelacao entre h(t) e h(t− 1).

36

Com as suposicoes (3.1), (3.2) e (3.3) temos:

y(t) ∼ N[0, σ2

ε exp(h(t))]

(3.4)

h(1) ∼ N[µ, σ2

η

]

h(t) | h(t− 1) ∼ Nµ+ φ (h(t− 1)− µ) ;σ2

η

para t = 2, 3, ..., N .

Uma possıvel generalizacao para o modelo de volatilidade estocastica definido pelas

expressoes (3.1) e (3.2) pode ser obtida definindo-se um modelo para as variaveis latentes

dado por:

h(t) = µ+

p∑j=1

φj (h(t− j)− µ) + η(t) (3.5)

para t = p + 1, ..., N, com as raızes do polinomio φ(B) = 1 −∑p

j=1 φjBj fora do cırculo

de raio unitario (B e o operador retardo definido por Bkh(t) = h(t − k)). Denota-

se o modelo definido por (3.1), (3.2) e (3.5) como modelo generalizado de volatilidade

estocastica. Neste caso temos:

h(t) | h(t− 1), ..., h(t− p) ∼ N

µ+

p∑j=1

φj (h(t− j)− µ) ;σ2η

(3.6)

para t = 2, 3, ..., N .

A funcao de verossimilhanca do modelo generalizado de volatilidade estocastica con-

siderando as equacoes (3.1) e (3.2) e dada por:

L =T∏i=1

p (y(t) | h(t))

37

Da equacao (3.4) temos que:

L =T∏i=1

1√2πσ2

ε exp(h(t))exp

[− y2(t)

2σ2ε exp(h(t))

]

3.2.1 Analise Bayesiana

Para uma analise Bayesiana do modelo de volatilidade estocastica definido em (3.1) e

(3.2) com as variaveis latentes definidas por (3.5) vamos assumir as seguintes distribuicoes

a priori independentes para os parametros Φ =(φ1, ..., φp

)′, σ2

ε , σ2η e µ :

φj ∼ U (−1, 1) , j = 1, ..., p

σ2η ∼ IG (c1, d1)

σ2ε ∼ IG (c2, d2)

µ ∼ N(e; f 2

)sendo que U (−1, 1) denota uma distribuicao Uniforme em (−1, 1), N(e, f 2) denota uma

distribuicao normal com media e e variancia f 2 e IG (c, d) denota uma distribuicao gama

inversa com media d/ (c− 1) e variancia d2/[(c− 1)2 (c− 2)

], c > 2. Vamos considerar

os hiperparametros (ci, di, e), i = 1, 2; j = 1, ..., p conhecidos. Fazendo hl = 0 para

l = 0,−1, ...,−p+ 1 na equacao (3.6), a densidade condicional de h(t) dado h(p)(t− 1) =

(h(t− 1), ..., h(t− p)), para t = 1, ..., N , e dada por:

p[h(t) | h(p)(t− 1)

]=

1√2πσ2

η

exp

− 1

2σ2η

[h(t)− µ−

p∑j=1

φj (h(t− j)− µ)

]2Denotando por Φ =

(φ1, ..., φp

)′, θ = (µ,Φ,σ2

ε , σ2η)′, h = (h(1), ..., h(n))′ e assumindo

independencia a priori entre os parametros µ,Φ,σ2ε , σ

2η a densidade a posteriori conjunta

38

para ϕ = (θ,h) e dada por:

π (ϕ | y) ∝N∏t=1

1√2πσ2

ε exp(h(t))exp

[− y2(t)

2σ2ε exp(h(t))

]×

× 1√2πσ2

η

exp

− 1

2σ2η

[h(t)− µ−

p∑j=1


]2 π (θ)

Podemos escrever a densidade a posteriori como:

π (ϕ | y) ∝ (σ2ε)−N

2 exp

[−1

2

N∑t=1

h(t)− 1

2σ2ε

N∑t=1

y2(t) exp(−h(t))

]×

×(σ2η)−N

2 exp− 1

2σ2η

N∑t=1

[h(t)− µ−

p∑j=1


]2 π (θ)

Considerando a notacao vetorial Φ =(φ1, ..., φp

)′, Y = (y2(1), y2(2),..., y2(N)), h =

(h(1), ..., h(n))′, G(h) = [exp(−h(1)), exp(−h(2)), ..., exp(−h(N))]′ e a matriz X dada

por:

X =

h(0) h(−1) · · · h(−p+ 1)

h(1) h(0) · · · h(−p+ 2)

· · · ·

· · · ·

· · · ·

h(p) h(p− 1) · · · h(1)

h(p+ 1) h(p) · · · h(2)

· · · ·

· · · ·

· · · ·

h(N − 1) h(N − 2) · · · h(N − p)

N×p

39

Podemos escrever a densidade a posteriori na forma,

π (ϕ | y) ∝ (σ2ε)−N

2 exp

[−1

21′h− 1

2σ2ε

Y′G(h)

]×

×(σ2η)−N

2 exp

[− 1

2σ2η

(h− µ−XΦ)′ (h− µ−XΦ)

]π (θ)

em que 1 = (1, 1, ..., 1)′ e o vetor N × 1 de 1′s e µ e um vetor N × 1 de µ′s.

Para gerar amostras da distribuicao a posteriori conjunta para ϕ = (θ,h) usamos

metodos de MCMC com algoritmo Gibbs Sampling (GELFAND e SMITH, 1990) ou o

algoritmo Metropolis-Hastings (SMITH e ROBERTS, 1993). Essas amostras sao geradas

a partir das distribuicoes condicionais π[θj | θ(j),y

]em que θ(j) denota o vetor de todas

as componentes de θ exceto o j - esimo componente.

Uma grande simplificacao e obtida usando o software WinBugs (SPIEGELHALTER

et al., 1999), no qual apenas necessitamos definir a distribuicao para os dados e as dis-

tribuicoes a priori para os parametros do modelo.

3.3 Modelos de Volatilidade Estocastica Bivariados

Sejam y(t) = (y1(t), y2(t))′ , t = 1, 2, ..., N , duas series temporais onde y1(t) e y2(t),

usualmente sao os logaritmos dos retornos centralizados em suas medias, com modelo,

y(t) = H(t)ε(t) (3.7)

onde ε(t) = (ε1(t), ε2(t))′ e o vetor de componentes de erros para o vetor de series e

H(t) e uma matriz 2× 2 diagonal dada por,

H(t) = diag(eh1(t)/2, eh2(t)/2

)(3.8)

Isto e,

40

y1(t) = eh1(t)/2ε1(t) (3.9)

y2(t) = eh2(t)/2ε2(t)

para t = 1, 2, ..., N.

Assumir que ε(t) = (ε1(t), ε2(t))′ tem uma distribuicao normal bivariada, dada por,

ε(t) ∼ N (0,Σε) (3.10)

onde a matriz de variancia-covariancia e dada por,

Σε =

1 ρε

ρε 1

(3.11)

Tambem assumir que h(t) = (h1(t), h2(t)) , t = 1, 2, ..., N e um vetor de variaveis

latentes definidas por um modelo AR(1),

h1(1) = µ1 + η1(1) (3.12)

h2(1) = µ2 + η2(1)

h1(t) = µ1 + φ1,1 (h1(t− 1)− µ1) + η1(t)

h2(t) = µ2 + φ2,2 (h2(t− 1)− µ2) + η2(t)

para t = 2, 3, ..., N ; 0 < φ1,1 < 1; 0 < φ2,2 < 1.

η(t) = (η1(t), η2(t)) e suposta com uma distribuicao normal dada por,

η(t) ∼ N

0; diag(σ2η1, σ2

η2

)(3.13)

41

onde 0 = (0, 0)′ ; e diag(σ2η1, σ2

η2

)e uma matriz 2× 2 diagonal.

De (3.7), observar que E (y(t)) = 0 a matriz de variancia-covariancia para y(t) e dada

por,

var (y(t)) = H′(t)ΣεH(t) =

eh1(t) ρεeh1(t)/2eh2(t)/2

ρεeh1(t)/2eh2(t)/2 eh2(t)

(3.14)

Assim, y(t) = (y1(t), y2(t))′ , tem uma distribuicao normal bivariada com densidade,

f (y1(t), y2(t) | h1(t), h2(t)) =1

2π√

(1− ρ2ε) eh1(t)+h2(t)× (3.15)

× exp

− 1

2 (1− ρ2ε)

[y21(t)

eh1(t)+y22(t)

eh2(t)− 2ρεy1(t)y2(t)

eh1(t)/2eh2(t)/2

]

Observar que as funcoes de h(t) = (h1(t), h2(t))′ fornecem a volatilidade da serie

temporal.

Tambem observar de (3.12) e (3.13), que as variaveis latentes h(t) = (h1(t), h2(t))′ tem

distribuicoes normais,

hi(1) ∼ N(µi;σ

2ηi

)(3.16)

hi(t) | hi(t− 1) ∼ Nµi + φi,i (hi(t− 1)− µi) ;σ2

ηi

para i = 1, 2; t = 2, ..., N.

Serao utilizados 6 modelos diferentes considerando a estrutura bivariada apresentada

acima, para modelar a volatilidade das aplicacoes a serem estudadas posteriormente.

Modelo 1 Um caso especial do modelo de volatilidade estocastica bivariada apresentado

acima e quando assumimos que os componentes de erros ε(t) sao independentes, isto e

42

ρε = 0. Portanto y(t) seguira uma distribuicao Normal com um vetor de medias 0 e matriz

de variancia-covariancia diagonal dada por Σε = diag(eh1(t), eh2(t)

).

Modelo 2 Nesta versao do modelo de volatilidade estocastica bivariada, a covariancia

entre os componentes do erro, ρε, e considerada como uma quantidade desconhecida que

precisa ser estimada. Portanto y(t) seguira uma distribuicao Normal com um vetor de

medias 0 e matriz de variancia-covariancia diagonal dada por (3.14).

Modelo 3 Em uma terceira versao do modelo sera mantido o que foi assumido no

modelo 2, exceto pela maneira como definimos a variavel latente h2(t), t = 1, 2, ..., N .

Nesta versao, e assumida presenca da causalidade de Granger quando modelamos h2(t),

isto e, a variavel latente h2(t) e agora dada por:

h2(t) = µ2 + φ2,1 (h1(t− 1)− µ1) + φ2,2 (h2(t− 1)− µ2) + η2(t) (3.17)

onde t = 2, 3, ..., N e −1 < φ2,1, φ2,2 < 1.

Observacao: Notar que de (3.17), se φ2,1 6= 0, entao a volatilidade do segundo retorno e

considerada com a causalidade de Granger induzida pela volatilidade do primeiro retorno.

Modelo 4 Nesta versao do modelo sera tomado o que foi assumido no modelo 2 exceto

que sera considerada uma hipotese adicional na covariancia entre ε1(t) e ε2(t). Nesta

direcao, sera assumido que ε(t) = (ε1(t), ε2(t))′ siga uma distribuicao normal bivariada

com um vetor de medias 0 e a matriz de variancia-covariancia Σε(t) e uma matriz 2 × 2

dada, para t = 1, 2, ..., N, por:

Σε(t) =

1 ρε(t)

ρε(t) 1

onde ρε(t) =

(eq(t) − 1

)/(eq(t) + 1

), q(1) = ψ0 + σρυ(t) e para t = 1, 2, ..., N, q(t) =

ψ0 + ψ1 [q(t− 1)− ψ0] + σρυ(t), com υ(t), t = 1, 2, ..., N, quantidades independentes e

43

identicamente distribuıdas com uma mesma distribuicao Normal N (0, 1) (YU e MEYER,

2005). Tambem sera assumido que q(1) siga uma distribuicao Normal N(ψ0, σ

2ρ

)e para

t = 2, 3, ..., N pode-se tambem, assumir que dado q(t − 1), a quantidade q(t) tera uma

distribuicao Normal N(ψ0 + ψ1 [q(t− 1)− ψ0] , σ

2ρ

).

Observacao: Neste modelo de volatilidade estocastica com correlacao dinamica entre

os componentes do erro pode-se observar que a covariancia entre estes muda ao longo do

tempo. Tambem observar que e necessario ter −1 < ρε(t) < 1 para se ter uma matriz de

variancia-covariancia Σε(t) bem definida.

Modelo 5 Serao consideradas as mesmas hipoteses de modelo 4 com a excecao da

hipotese feita para h2(t). No caso desta variavel latente sera considerado o que foi assumido

no modelo 3, a presenca da causalidade de Granger, isto e h2(t) e dado por (3.17).

Modelo 6 Neste modelo sera considerado um ajuste similar ao do modelo 1. A diferenca

e que agora os componentes do erro ε(t) seguem uma distribuicao t-Student multivari-

ada com v graus de liberdade, t = 1, 2, ..., N. Por isso, sera assumido que ε(t) segue

uma distribuicao t-Student multivariada com vetor de medias 0 e a matriz de variancia-

covariancia [v/ (v − 2)] Σε, para v > 2. Portanto, o vetor ε(t) = (ε1(t), ε2(t))′ tem como

funcao densidade uma t (0,Σε, v), onde 0 representa o vetor de medias, Σε a matriz de

variancia-covariancia e v os graus de liberdade, isto e, sua funcao de densidade e dada

por:

f (ε(t)) =Γ [(v + 2) /2]

vπΓ [v/2]|Σε|−1/2

(1 +

1

vε′(t)Σ−1ε ε(t)

)−(v+2)/2

(3.18)

onde Γ (x) denota uma funcao Gama.

Observacao: Usando esta distribuicao t-Student de caudas pesadas, pode-se ter a

presenca de extra Curtose para as distribuicoes dos retornos. Tambem e interessante ob-

servar que muitas outras modificacoes para os modelos de volatilidade estocastica podem

ser consideradas.

44


A analise Bayesiana para os modelos de volatilidade estocastica bivariados sera divi-

dida em tres categorias. Aqueles na assim chamada Classe 1 sao os modelos que possuem

o vetor de erros ε(t) = (ε1(t), ε2(t))′ normalmente distribuıdos com vetor de medias 0 e

uma correlacao constante para os componentes do erro do modelo, isto e a correlacao e

zero ou uma constante nao dependendo de t. Os modelos que se encaixam nessa classe

sao os modelos 1, 2 e 3. Na Classe 2, nos incluımos os modelos onde ε(t) e Normalmente

distribuıdo com vetor de medias 0, mas a correlacao entre ε1(t) e ε2(t) tambem depende

do tempo, isto e os modelos 4 e 5. Finalmente, na Classe 3 nos temos o modelo 6.

A inferencia Bayesiana sera aplicada atraves de uma amostra retirada da distribuicao

a posteriori dos parametros do modelo. Entretanto, nos precisamos especificar em cada

caso quais as distribuicoes a priori nos estamos considerando e quais as funcoes de verossi-

milhanca utilizadas. Em todos os casos nos assumimos independencia entre os parametros

dos modelos.

Classe 1 Se o modelo 1 for considerado entao o vetor de parametros a ser estimado e

θI = (φ1,1, φ2,2, σ2η1, σ2

η2, µ1,µ2). Nos assumimos que φi,i, σ

2ηi

e µi tem como distribuicoes

a priori uma Uniforme, uma Gama Invertida e uma Normal, respectivamente, i = 1, 2,

isto e, φi,i, σ2ηi

e µi tem, respectivamente, distribuicoes a priori U(aii, bii), IG(ci, di) e

N(ei, f2i ), onde os hiperparametros aii, bii, ci, di, ei e f 2

i sao conhecidos, i = 1, 2.

Quando o modelo 2 e considerado, o vetor de parametros e θII = (θI, ρε). Nos assumi-

mos que ρε tem como distribuicao a priori uma distribuicao Uniforme U(−1, 1) e usamos

as mesmas distribuicoes a priori do modelo 1 para θI, com possıveis diferentes valores

para os hiperparametros.

Nos modelos 1 e 2, as funcoes densidades conjuntas das variaveis latentes h(t) =

(h1(t), h2(t)) dado o vetor de parametros sao dadas por

45

g (h(1) | θ) ∝2∏i=1

(σ2ηi

)−1/2exp

[− 1

2σ2ηi

(hi(1)− µi)2

], t = 1 (3.19)

g (h(t) | h(t− 1),θ) ∝2∏i=1

N∏t=2

(σ2ηi

)−1/2exp

− 1

2σ2ηi

[hi(t)− µi−

−φi,i (hi(t− 1)− µi)]2

, t = 2, 3, ..., N

Seja θ = θII e tomando ϕ = (θ,h) , com h(t) = (h(1), h(2), ..., h(N)) . Assim, nos

temos que a funcao de verossimilhanca conjunta de θ e h para o modelo 2 e dada, para

y(t) = (y(1), y(2), ..., y(N)) , por

L(ϕ | y(t)) ∝N∏t=1

p (y(t) | h(t),θ) ∝ (3.20)

∝(1− ρ2ε

)−N/2exp

−1

2

[N∑t=1

h1(t) +N∑t=1

h2(t)

]×

×

exp− 1

2(1− ρ2ε)

[N∑t=1

y21(t)e−h1(t) +N∑t=1

y22(t)e−h2(t)−

−2ρε

N∑t=1

y1(t)y2(t)e−h1(t)/2e−h2(t)/2

]

No caso de θ = θI nos simplesmente ajustamos ρε = 0 em (3.20).

Entretanto, para ϕ = (θI,h) ou ϕ = (θII,h) nos temos que a distribuicao a posteriori

conjunta do vetor de parametros e h e dada por

π (ϕ | y(t)) ∝ π (θ) g (h(1) | θ)N∏t=2

g (h(t) | h(t− 1),θ)L(ϕ | y(t)) (3.21)

onde π (θ) e a distribuicao a priori do vetor de parametros com θ = (θI,θII) , L(ϕ | y(t)) e

a funcao de verossimilhanca do modelo dada por (3.20), e g (h(1) | θ) , g (h(t) | h(t− 1),θ),

t = 2, 3, ..., N sao dadas pelo conjunto de funcoes recursivas (3.19).

Quando assumimos a presenca da causalidade de Granger para h2(t) e uma correlacao

46

constante para os componentes do erro, isto e, o modelo 3, o vetor de parametros e θIII =

(θII, φ2,1). Nos tomamos as mesmas distribuicoes a priori para θII com possıveis diferentes

hiperparametros. Adicionalmente para a variavel φ2,1 nos tomamos como sua distribuicao

a priori uma distribuicao U (a21, b21) com hiperparametros a21 e b21 conhecidos.

A funcao de verossimilhanca do modelo 3 e dada por (3.20). A distribuicao a posteriori

conjunta de ϕ = (θ,h) tambem tem a mesma expressao como nos modelos 1 e 2, isto

e, a expressao e dada por (3.21), mas agora tomando θ = (θIII) e substituindo a funcao

densidade g (h(t) | h(t− 1),θ) , t = 2, 3, ..., N por

g (h(t) | h(t− 1),θ) ∝

(N∏t=2

(σ2η1

)−1/2× (3.22)

× exp

− 1

2σ2η1

[h1(t)− µ1 − φ1,1 (h1(t− 1)− µ1)

]2)××

(N∏t=2

(σ2η2

)−1/2exp

− 1

2σ2η2

[h2(t)− µ2−

−φ2,1 (h1(t− 1)− µ1)− φ2,2 (h2(t− 1)− µ2)]2)

Classe 2 Os modelos nesta classe tem um modelo de volatilidade estocastica com cor-

relacao dinamica para os componentes dos erros. Quando o modelo 4 e considerado nos

tempos que o vetor de parametros e θIV = (θI, ψ0, ψ1, σ2ρ). A distribuicao a priori para

os componentes de θI sao as mesmas tomadas na classe 1. Os parametros ψ0, ψ1 e

σ2ρ tem como distribuicoes a priori uma distribuicao Normal N(g, h2), uma distribuicao

Beta(m,n) e uma distribuicao Gama Invertida IG(c3, d3), onde os hiperparametros c3, d3

, g, h, m e n sao considerados conhecidos.

Assim, a funcao de verossimilhanca conjunta desta classe e dada, para θ = θIV, por,

47

L(ϕ | y(t)) ∝N∏t=1

p (y(t) | h(t),θ) ∝

N∏t=1

(1− ρ2ε(t)

)−1/2× (3.23)

× exp

−1

2

[N∑t=1

h1(t) +N∑t=1

h2(t)

]×

×

exp−1

2

[N∑t=1

y21(t)e−h1(t)

1− ρ2ε(t)+

N∑t=1

y22(t)e−h2(t)

1− ρ2ε(t)−

−2N∑t=1

ρε(t)y1(t)y2(t)e−h1(t)/2e−h2(t)/2

1− ρ2ε(t)

]

A distribuicao a posteriori conjunta de ϕ = (θ,h,q) , onde θ = θIV, q = (q(1), q(2), ...,

q(N)), e dada por

π (ϕ | y(t)) ∝ π (θ) g (h(1) | θ)N∏t=2

g (h(t) | h(t− 1),θ)× (3.24)

×f (q(1) | θ)

(N∏t=2

f (q(t) | q(t− 1),θIV)

)L(ϕ | y(t))

onde π (θ) e a distribuicao a priori conjunta para θ; g (h(1) | θ) e g (h(t) | h(t− 1),θ)

sao definidas por (3.19), f (q(1) | θ) e f (q(t) | q(t− 1),θIV) sao as funcoes densidades

Normal de q(1) e a funcao densidade condicional Normal de q(t) dado q(t− 1), e L(ϕ |

y(t)) e a funcao de verossimilhanca definida em (3.23).

Quando consideramos o modelo de volatilidade estocastica com a correlacao dinamica

para os componentes do erro e a presenca da causalidade de Granger para h2(t), t =

2, 3, ..., N, isto e, modelo 5, o vetor de parametros aqui e θV = (θIV, φ2,1). Nos tomamos

a mesma distribuicao a priori para θIV exceto para os parametros ψ0 e ψ1 onde agora

terao como distribuicoes a priori, distribuicoes Gama com hiperparametros apropriados.

Nos continuaremos tomando uma distribuicao a priori U (a21, b21) para φ2,1.

A verossimilhanca do modelo e tambem dada por (3.23). A distribuicao a posteriori

e similar a (3.24), a diferenca e que g (h(t) | h(t− 1),θ) , t = 2, 3, ..., N sao dadas por

48

(3.22) com θ substituido por θV.

Classe 3 Nesta classe de modelos nos assumimos que os componentes do erro tem

uma distribuicao t-Student bivariada com grau de liberdade > 2. As variaveis latentes

h(t) = (h1(t), h2(t)) , t = 1, 2, ..., N sao definidas como no modelo 1.

A funcao de verossimilhanca pode ser obtida diretamente de (3.18) ou representando a

distribuicao t-Student bivariada como uma mistura de uma distribuicao Normal com uma

distribuicao Gama (BERNARDO e SMITH, 1995). Seguindo esta ultima representacao,

primeiro nos temos que a distribuicao para ε(t) e dada por,

fε (ε(t)) =

∫fNormal (ε(t) | µ,Σyz(t)) fGama (z(t) | α, β) dz(t) (3.25)

onde fNormal (ε(t) | µ,Σ) denota uma funcao densidade Normal bivariada com vetor de

medias µ e matriz de variancia-covariancia Σ e fGama (z | α, β) e uma densidade Gama

com media α/β e variancia α/β2.

Entretanto, tomar z(t) = 1/√w(t), onde w(t) segue uma distribuicao Gama(λw, λw)

com E(w(t)) = 1 e var(w(t)) = λ−1w . Assim, a distribuicao condicional de y(t) dado Σy e

w(t) e uma distribuicao Normal bivariada com vetor de medias 0 e matriz de variancia-

covariancia 1√w(t)

Σy, onde Σy e dada por (3.14).

O vetor de parametros para o modelo 6 e θVI = (θII,h,w) e, para w =(w(1), w(2), ...,

w(N)), a funcao de verossimilhanca do modelo e dada, para θVI = (θII,h,w), por

49

L(ϕ | y(t)) ∝N∏t=1

p (y(t) | h(t),w(t)) (3.26)

∝(1− ρ2ε

)−N/2 N∏t=1

w1/2(t)

exp

−1

2

[N∑t=1

h1(t) +N∑t=1

h2(t)

]×

× exp

− 1

2 (1− ρ2ε)

[N∑t=1

y21(t)w1/2(t)e−h1(t)+

+N∑t=1

y22(t)w1/2(t)e−h2(t) − 2ρε

N∑t=1

y1(t)y2(t)w1/2(t)e−h1(t)/2e−h2(t)/2

]

As distribuicoes a priori de θII sao as mesmas usadas nos modelos da Classe 1 (com

possıveis hiperparametros diferentes) e nos tomamos uma distribuicao a priori Gama(o, p)

para w. Nos tambem temos que o e p sao hiperparametros conhecidos.

Assim, a distribuicao a posteriori conjunta de θVI = (θII,h,w), e dada por

π (ϕ | y(t)) = π (θVI) g (h(1) | θVI)

(N∏t=2

g (h(t) | h(t− 1),θVI)

)× (3.27)

×

(N∏t=1

gw (w(t) | θVI)

)L(ϕ | y(t))

onde π (θVI) e a distribuicao a priori conjunta de θVI, g (h(1) | θVI) sao dadas por (3.19),

eN∏t=1

gw (w(t) | θVI) e o produto de densidades Gama(λw, λw), isto e,

N∏t=1

g (w(t) | θVI) ∝N∏t=1

wλw−1(t)e−λww(t)

e L(ϕ | y(t)) e a funcao de verossimilhanca (3.26).

50

3.4 Modelos de Volatilidade Estocastica Multivaria-

dos

Alem dos modelos bivariados, tambem serao utilizados modelos de volatilidade es-

tocastica multivariados para analisar os dados de poluicao a serem apresentados.

Seja y(t) = (y1(t), y2(t), ..., yk(t))′ , t = 1, 2, ..., N um vetor de series temporais, onde

y1(t), y2(t), ..., yk(t), usualmente sao os logaritmos dos retornos centralizados em suas me-

dias, com modelo,

y(t) = H(t)ε(t) (3.28)

onde ε(t) = (ε1(t), ε2(t), ..., εk(t))′ e o vetor de componentes de erros para o vetor de

series e H(t) e uma matriz k × k diagonal dada por,

H(t) = diag(eh1(t)/2, eh2(t)/2, ..., ehk(t)/2

)(3.29)

Isto e,

yi(t) = ehi(t)/2εi(t) (3.30)

para i = 1, 2, ..., k; t = 1, 2, ..., N.

Assumir que ε(t) = (ε1(t), ε2(t), ..., εk(t))′ tem uma distribuicao normal multivariada,

dada por,

ε(t) ∼ N (0,Σε) (3.31)

onde a matriz de variancia-covariancia e dada por,

51

Σε =

1 ρ1,2 ρ1,3 · · · ρ1,k

1 ρ2,3 · · · ρ2,k

· ·

· ·

· ·

1 ρk−1,k

1

(3.32)

onde ρij, i, j = 1, 2, ..., k; i 6= j e a covariancia entre εi(t) e εj(t). Tambem observar

que estamos assumindo que as variancias dos erros εi(t), i = 1, 2, ..., k sao iguais a 1.

Tambem assumir que h(t) = (h1(t), h2(t), ..., hk(t)) , t = 1, 2, ..., N e um vetor de vari-

aveis latentes definidas por um modelo AR(1),

hi(1) = µi + ηi(1) (3.33)

hi(t) = µi + φi,i (hi(t− 1)− µi) + ηi(t)

para i = 1, 2, ..., k; t = 2, 3, ..., N ;∣∣φi,i∣∣ < 1.

η(t) = (η1(t), η2(t), ..., ηk(t)) e suposta com uma distribuicao normal dada por,

η(t) ∼ N

0; diag(σ2η1, σ2

η2, ..., σ2

ηk

)(3.34)

onde 0 = (0, 0, ..., 0)′ ; e diag(σ2η1, σ2

η2, ..., σ2

ηk

)e uma matriz k × k diagonal.

De (3.28), observar que E (y(t)) = 0 a matriz de variancia-covariancia para y(t) e

dada por,

52

var (y(t)) =

eh1(t) ρ1,2eh1(t)/2eh2(t)/2 · · · ρ1,ke

h1(t)/2ehk(t)/2

eh2(t) · · · ρ2,keh2(t)/2ehk(t)/2

· ·

· ·

· ·

ehk(t)

(3.35)

Assim, y(t) = (y1(t), y2(t), ..., yk(t))′ , tera distribuicao normal multivariada como em

(3.15).

Observar que as funcoes de h(t) = (h1(t), h2(t), ..., hk(t)) fornecem a volatilidade da

serie temporal.

Tambem observar de (3.33) e (3.34), que as variaveis latentes h(t) = (h1(t), h2(t), ...,

hk(t)) tem distribuicoes normais,

hi(1) ∼ N(µi;σ

2ηi

)(3.36)

hi(t) | hi(t− 1) ∼ Nµi + φi,i (hi(t− 1)− µi) ;σ2

ηi

para i = 1, 2, ..., k; t = 2, ..., N.

Devido ao grande custo computacional, nas aplicacoes a serem apresentadas nos capı-

tulos a seguir, serao utilizados apenas os modelos 1 e 2 apresentados na secao anterior.


Se o modelo 1 for considerado entao o vetor de parametros a ser estimado e θI =

(φ,σ2η,µ), onde φ = (φ11, φ22, ..., φkk) , σ

2η =

(σ2η1, σ2

η2, ..., σ2

ηk

)e µ = (µ1, µ2, ..., µk). As-

sumimos que φii, σ2ηi

e µi tem como distribuicoes a priori uma Beta, uma Gama Invertida

53

e uma Normal, respectivamente, i = 1, 2, ..., k, isto e, φii, σ2ηi

e µi tem, respectivamente,

distribuicoes a priori U(aii, bii), IG(ci, di) e N(ei, f2i ), onde os hiperparametros aii, bii, ci,

di, ei e fi sao conhecidos, i = 1, 2, ..., k.

Quando o modelo 2 e considerado o vetor de parametros e θII = (θI,ρ), onde ρ = (ρ12,

ρ13, ..., ρk−1,k).Assumimos que ρij tem como distribuicao a priori uma distribuicao Uni-

forme U(gi, hi) e usamos as mesmas distribuicoes a priori para θI do modelo 1, com

possıveis diferentes valores para os hiperparametros.

Nos modelos 1 e 2, as funcoes densidades conjuntas das variaveis latentes h(t) =

(h1(t), h2(t), ..., hk(t)) dado o vetor de parametros sao dadas por

g (h(1) | θ) ∝k∏i=1

(σ2ηi

)−1/2exp

[− 1

2σ2ηi

(hi(1)− µi)2

], t = 1 (3.37)

g (h(t) | h(t− 1),θ) ∝k∏i=1

N∏t=2

(σ2ηi

)−1/2exp

− 1

2σ2ηi

[hi(t)− µi−

−φi,i (hi(t− 1)− µi)]2

Seja θ = θII e tomando ϕ = (θ,h) , com h(t) = (h1(t), h2(t), ..., hk(t)) , onde t =

1, 2, ..., N. Assim, temos que a funcao de verossimilhanca conjunta de θ e h para o modelo

2 sera definida, para y(t) = (y1(t), y2(t), ..., yk(t)) , onde t = 1, 2, ..., N , por

L(ϕ | y(t)) ∝N∏t=1

p (y(t) | h(t),θ) (3.38)

onde p (y(t) | h(t),θ) e dado em (3.15)

No caso de θ = θI nos simplesmente ajustamos ρij = 0 em (3.38).

Entretanto, para ϕ = (θI,h) ou ϕ = (θII,h) nos temos que a distribuicao a posteriori

conjunta do vetor de parametros e h e dada por

54

π (ϕ | y(t)) ∝ π (θ) g (h(1) | θ)N∏t=2

g (h(t) | h(t− 1),θ)L(ϕ | y(t)) (3.39)

onde π (θ) e a distribuicao a priori do vetor de parametros com θ = (θI,θII) , L(ϕ | y(t)) e

a funcao de verossimilhanca do modelo dada por (3.38), e g (h(1) | θ) , g(h(t) | h(t−1),θ),

t = 2, 3, ..., N sao dadas pelo conjunto de funcoes recursivas (3.37).

3.5 Criterios para Selecao de Modelos

Varios criterios para escolha do modelo de volatilidade estocastica que melhor se ajusta

a uma serie de dados podem ser adotados no contexto Bayesiano (BERG et al. (2004)).

Dentre estes vamos usar aqui o seguinte criterio:

Criterio Desvio-Informacao (DIC) O DIC (Deviance Information Criterion) e uma

generalizacao do BIC (Criterio de Informacao Bayesiana). Este criterio e particularmente

usual nos problemas Bayesianos de selecao de modelos para os quais amostras da dis-

tribuicao a posteriori dos parametros dos modelos foram obtidas por simulacao de Monte

Carlo em Cadeia de Markov (MCMC). Semelhante ao BIC este criterio e uma aproxi-

macao assintotica para amostras grandes e e valido quando a distribuicao a posteriori e

aproximadamente uma distribuicao normal multivariada.

Define-se o desvio como:

D (θ) = −2 lnL (θ) + C (3.40)

em que θ e o vetor de parametros desconhecidos do modelo e L (θ) a funcao de verossi-

milhanca. O criterio DIC introduzido por Spiegelhalter et al. (2002) e dado por:

DIC = D (θ) + 2pD (3.41)

sendo D (θ) o desvio avaliado na media a posteriori e pD e o numero efetivo de parametros

55

no modelo, que e dado por pD = D − D (θ), em que D = E[D(θ)] e o desvio medio a

posteriori que mede a qualidade do ajuste do modelo aos dados. Valores menores para

DIC indicam melhores modelos, podendo estes serem negativos. O DIC e usado para a

discriminacao de diferentes modelos propostos para a analise de series financeiras.

56

4 Aplicacoes

4.1 Analise dos Dados de Poluicao da Cidade do Me-

xico

4.1.1 Descricao dos Dados

A Cidade do Mexico possui serios problemas com a poluicao do ar, especialmente com

o ozonio.As autoridades ambientais tem implementado medidas buscando reduzir o nıvel

dos poluentes em geral. Tais medidas sao importantes porque quando a concentracao de

ozonio fica acima de um dado limite por um certo perıodo de tempo, os indivıduos expostos

ao poluente podem sofrer serios problemas de saude (WILSON et al., 1980; LOOMIS et

al., 1996; BELL et al., 2004; BELL et al., 2005; BELL et al., 2007). Em particular, sabe-se

que para nıveis de ozonio acima de 0,11 parte por milhao (0,11ppm) uma porcao muito

sensıvel da populacao (como os recem-nascidos e os mais velhos) sofre uma deterioracao

de sua saude. Portanto, ser capaz de entender o comportamento a longo termo deste

poluente e de grande importancia para o planejamento de polıticas publicas de saude. Se

existe uma tendencia nas medidas em diminuir com o tempo com variacoes nao muito

grandes, ou por outro lado se estas estao aumentando ou apresentam grandes variacoes,

entao as autoridades ambientais podem ter uma ideia se, no longo perıodo de observacao,

as medidas preventivas tomadas por elas estao produzindo o resultado desejado, ou nao.

Diversos metodos tem sido utilizados para predizer a violacao de um padrao de quali-

dade do ar. Entre eles podemos nos referir a Roberts (1979a, 1979b), Horowitz (1980) e

Smith (1989) quando o interesse consiste em aplicar a teoria de valores extremos para re-

alizar previsoes. Entretanto, outras tecnicas tambem podem ser utilizadas para estes tipos

de problemas. Como exemplos nos podemos citar, analise multivariada (GUARDANI et

57

al., 2003), redes neurais (GUARDANI et al., 1999), modelos de Poisson (JAVITS, 1980;

RAFTERY, 1989; LEADBETTER, 1991; ACHCAR et al., 2007) e modelos de cadeia de

Markov (LARSEN et al., 1990; AUSTIN e TRAN, 1999; ALVAREZ et al., 2005).

Outra possibilidade e considerar a modelagem das medias diarias ou semanais da serie

temporal (LOOMIS et al., 1996; LANFREDI e MACCHIATO, 1997; ZOLGHADRI e

HENRY, 2004; PAN e CHEN, 2007).

Aqui nos utilizaremos os modelos de volatilidade estocastica (GHYSELS et al., 1996;

KIM et al., 1998; MEYER e YU, 2000) para estudar o comportamento das medias mensais

das medidas de ozonio obtidas pela rede de monitoramento da Cidade do Mexico.

Os dados foram obtidos atraves da rede de monitoramento da Area Metropolitana da

Cidade do Mexico (www.sma.df.gob.mx/simat/). A Area Metropolitana e dividida em

cinco regioes, ou setores correspondendo a Nordeste (NE), Noroeste (NW), Centro (CE),

Sudeste (SE) e Sudoeste (SW) e as estacoes de monitoramento de ozonio estao localizadas

por toda a cidade (ver ALVAREZ et al., 2005; ACHCAR et al., 2007).

Os dados utilizados na analise correspondem a dezesseis anos de medicoes do maximo

diario de ozonio, indo de 01 de Janeiro de 1990 ate 31 de Dezembro de 2005, nas regioes

NE, CE e SW. As medidas sao obtidas minuto a minuto e o resultado da media horaria

e reportada em cada estacao. A medida do maximo diario para uma dada regiao e o

maximo entre todos os maximos dos valores medios registrados de hora em hora durante

um perıodo de 24 horas por cada estacao localizada na regiao. A medidas medias nos

dezesseis anos nas regioes NE, CE e SW sao 0,1, 0,141 e 0,132, respectivamente, com

desvios-padrao 0,042, 0,055 e 0,048. As regioes NE, CE e SW foram escolhidas devido

a direcao do vento na Area Metropolitana da Cidade do Mexico ser principalmente da

regiao NE para SW. Muitos dos precursores do ozonio sao produzidos nas regioes NE e

CE. Em geral, a regiao SW e a regiao com maior incidencia de nıveis elevados de ozonio.

Nos realizamos duas analises em separado. Na primeira analise nos consideramos os

dados das regioes NE e CE conjuntamente. Na segunda nos consideramos a informacao

conjunta provida pelas regioes CE e SW. A opcao por estes dois pares de regioes e devido

58

as caracterısticas geograficas dessas regioes. Essas regioes sao montanhosas e a direcao

dos ventos e um fator importante para definir esses pares. Os dados utilizados na analise

foram gerados pela media semanal das medidas de ozonio para as regioes NE, CE e SW.

Por isso, por exemplo, se x1(t), t = 1, 2, ..., N, indica as medidas dos maximos diarios de

ozonio na regiao NE, entao z1(t) e a medida media na semana t na regiao NE, isto e,

z1(1) = (1/7)7∑i=1

x1(i), z1(2) = (1/7)14∑i=8

x1(i), e assim por diante.

A inferencia foi realizada utilizando uma amostra gerada das distribuicoes a posteriori

conjuntas apos um perıodo de burn-in de tamanho 5000. Apos o perıodo de burn-in, para

a obtencao das estimativas dos parametros a posteriori, nos selecionamos 1000 amostras

de Gibbs tomando cada decima amostra com o objetivo de ter aproximadamente valores

nao correlacionados.

Na Figura 1, nos temos os graficos das medias semanais de ozonio versus as regioes

NE, CE e SW.

Observando a Figura 1 e possıvel notar uma tendencia decrescente para as medias

semanais de ozonio para as tres regioes da Cidade do Mexico, especialmente apos a

quadragesima semana (proxima ao ano de 1998).

Na Figura 3, nos temos os graficos dos log-retornos yi(t), t = 1, 2, ..., N ; i = 1, 2, 3,

centrados em suas medias para as regioes NE, CE e SW. (Aqui nos temos N = 384).

59

(a) NE (b) CE

(c) SW

Figura 3: Log-retornos centrados em suas medias para as regioes NE, CE e SW de 01 deJaneiro de 1990 a 31de Dezembro de 2005.

4.1.2 Uso dos Modelos e Resultados

Nos vamos analisar cada modelo e par de dados separadamente.

Modelo 1 Para se obter a convergencia do algoritmo foi assumido uma priori Beta(aii,bii)

para o parametro φi,i. Essa distribuicao tambem foi assumida nos demais modelos desta

aplicacao. Portanto, os hiperparametros das distribuicoes a priori sao dados como se

segue. Quando o par NE e CE ou o par CE e SW e considerado, nos temos que os hiper-

parametros das distribuicoes a priori de φi,i e µi sao aii = bii = 1, ci = di = 1, i = 1, 2.

Quando o par NE e CE e considerado, os hiperparametros da distribuicao a priori de µi

sao ei = 0, fi = 10, i = 1, 2. Entretanto, quando o par CE e SW e considerado, nos temos

ei = 0, fi = 1, i = 1, 2. Na Tabela 1 nos temos o sumario das quantidades estimadas

60

quando o par NE-CE e considerado assim como quando as regioes CE e SW sao levadas

em consideracao. (Note que na Tabela 1, e sempre que adequado, nos reportamos as

estimativas de τ ηi = 1/σ2ηi

e a estimativa de τ ρ = 1/σ2ρ ao inves de reportar as estimativas

de σ2ηi

e σ2ρ, i = 1, 2).

Tabela 1: Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de credibili-dade (95%) para as quantidades de interesse quando o modelo 1 e utilizado e os pares deregioes NE-CE e CE-SW sao considerados.

Regiao Parametro Media D.P. Intervalo de Credibilidade (95%)NE-CE τ η1 4,83 1,39 (2,86; 8,13)

τ η2 3,96 1,09 (2,39; 6,48)µ1 -3,07 0,06 (-3,19; -2,95)µ2 -2,99 0,08 (-3,13; -2,84)φ1,1 0,49 0,15 (0,18; 0,75)φ2,2 0,63 0,09 (0,46; 0,78)

CE-SW τ η1 3,92 0,90 (1,93; 5,21)τ η2 4,06 1,19 (2,26; 6,95)µ1 -2,94 0,08 (-3,09; -2,78)µ2 -2,97 0,08 (-3,12; -2,82)φ1,1 0,63 0,09 (0,44; 0,77)φ2,2 0,63 0,10 (0,41; 0,80)

Para o primeiro caso, onde estamos estudando o comportamento da volatilidade para

as regioes NE e CE, φ1,1 e um parametro auto-regressivo que indica que o valor da variavel

latente h para a regiao NE no instante de tempo t depende do valor de h no tempo t− 1.

Ja φ2,2 e o parametro auto-regressivo para a analise feita na regiao CE. Para o segundo

caso, φ1,1 representa o parametro auto-regressivo para a regiao CE e φ2,2 e o parametro

auto-regressivo para a regiao SW.

No primeiro caso vemos que a regiao CE apresenta valores diferentes para os parame-

tros φ e τ η em relacao a regiao NE, enquanto que no segundo caso, as regioes CE e SW

apresentam estimativas muito parecidas para os parametros analisados.

Na Figura 4, nos temos os graficos das raızes quadradas da volatilidade estimada

quando as regioes NE e CE (os dois graficos superiores) e quando as regioes CE e SW (os

dois graficos infeiores) sao levadas em conta quando o modelo 1 e utilizado.

61

(a) NE (b) CE

(c) CE (d) SW

Figura 4: Raızes quadradas das volatilidades para as regioes NE, CE e SW quando omodelo 1 e utilizado e quando as regioes sao pareadas em regioes NE e CE (linha superior)e regioes CE e SW (linha inferior).

Nos podemos observar atraves da Figura 4, que no primeiro par de regioes, a regiao

NE apresenta menor volatilidade ao longo do perıodo estudado, quando comparada a

regiao CE. E possıvel identificar a maior variacao e os picos que ocorrem na regiao CE.

Pode-se verificar tambem, que apos a quadragesima semana (proximo ao ano de 1998),

para ambas as regioes, ha uma queda e uma estabilizacao da volatilidade. O mesmo pode

ser observado para as regioes CE e SW. Pode-se observar tambem que as regioes CE e

SW apresentam nıveis bem similares de volatilidade.

Modelo 2 No caso do modelo 2, os hiperparametros das distribuicoes a priori sao dados

como se segue. Em ambos os casos, isto e, quando as regioes NE e CE sao consideradas

e quando o par CE e SW e levado em conta, nos temos que os hiperparametros das

distribuicoes a priori de φi,i e µi sao aii = bii = 1 e ei = −3, fi = 10, i = 1, 2. No caso dos

parametros σ2ηi, i = 1, 2, nos temos que se considerarmos as regioes NE e CE entao c1 = 5,

62

c2 = 4 e d1 = d2 = 1, i = 1, 2. Quando os dados das regioes CE e SW sao utilizados, nos

temos que c1 = 3, c2 = 4 e d1 = d2 = 1. Nos gostarıamos de chamar a atencao para o

fato que a mudanca nos hiperparametros foi feita, pois utilizamos a informacao provida

pelos resultados do modelo 1. Os novos valores dos hiperparametros foram calculados

igualando os valores obtidos pelo modelo 1 as respectivas media e variancia da distribuicao

a priori de interesse. Assim, nos estamos utilizando uma metodologia de Bayes empırica

para analisar o problema (ver CARLIN e LOUIS, 2000). Esta metodologia tambem foi

utilizada nos demais modelos desta aplicacao. Na Tabela 2 nos apresentamos o sumario

das quantidades estimadas.



τ η2 7,44 1,86 (4,66; 11,68)µ1 -3,02 0,07 (-3,16; -2,89)µ2 -2,92 0,07 (-3,07; -2,77)φ1,1 0,73 0,08 (0,54; 0,86)φ2,2 0,70 0,08 (0,53; 0,83)ρε 0,75 0,02 (0,72; 0,78)

CE-SW τ η1 7,96 1,49 (5,25; 10,99)τ η2 9,63 2,12 (6,43; 14,38)µ1 -2,85 0,06 (-2,98; -2,73)µ2 -2,89 0,06 (-3,02; -2,77)φ1,1 0,69 0,07 (0,54; 0,80)φ2,2 0,72 0,07 (0,56; 0,83)ρε 0,88 0,01 (0,87; 0,90)

Observando os resultados da Tabela 2 (modelo com a presenca do parametro ρε, que

captura a correlacao entre os componentes dos erros), podemos notar que as estimativas

para as variancias dos erros associados a variavel latente h (τ η) diferem bastante daquelas

encontradas no modelo 1, inclusive apresentando intervalos de credibilidade menos precisos

que no primeiro modelo estudado. Porem os intervalos de credibilidade para os parametros

auto-regressivos sao mais precisos utilizando este modelo. Neste modelo, as estimativas

para os parametros auto-regressivos das regioes NE e CE nao diferem tanto como no

modelo 1.

63

A Figura 5 mostra graficos similares aos da Figura 4, mas agora nos utilizamos o

modelo 2, ao inves do modelo 1.

(a) NE (b) CE

(c) CE (d) SW


Com base na Figura 5 e possıvel observar que o modelo 2 apresentou graficos mais

suaves que os obtidos pelo modelo 1. Os resultados obtidos por este modelo nao nos

mostram a grande diferenca dos valores de volatilidade entre as regioes NE e CE que o

modelo 1 registrou. Aqui, nota-se que os valores de volatilidade para este par de regioes

sao bem similares, com a regiao CE apresentando apenas algumas variacoes um pouco

mais elevadas em alguns perıodos. Neste caso esta mais evidente a queda e a estabilizacao

da volatilidade apos a quadragesima semana para os pares de regioes considerados.

Modelo 3 Considerando agora o modelo 3. Os hiperparametros das distribuicoes a

priori de φi,i, µi, i = 1, 2 e ρε tem distribuicoes a priori com os mesmos hiperparametros

64

para ambos os pares de regioes NE-CE e CE-SW. Os hiperparametros das distribuicoes

a priori de σ2η1

e σ2η2

sao c1 = 9, c2 = 7 e d1 = d2 = 1 quando as regioes NE e CE

sao consideradas e sao c1 = 8, c2 = 10 e d1 = d2 = 1 quando as regioes NE e SW sao

levadas em conta. A distribuicao a priori Beta da quantidade φ2,1 tem hiperparametros

a21 = b21 = 1 quando as regioes NE e CE ou as regioes CE e SW sao consideradas. A

Tabela 3 mostra o sumario das quantidades estimadas.



τ η2 9,23 2,31 (5,65; 14,32)µ1 -3,04 0,08 (-3,21; -2,88)µ2 -2,92 0,07 (-3,06; -2,78)φ1,1 0,83 0,04 (0,74; 0,90)φ2,2 0,40 0,20 (0,04; 0,75)φ2,1 0,49 0,23 (0,14; 0,9634)ρε 0,75 0,02 (0,71; 0,781)

CE-SW τ η1 13,83 2,76 (9,06; 19,34)τ η2 13,39 2,69 (8,94; 19,41)µ1 -2,84 0,10 (-3,04; -2,64)µ2 -2,90 0,09 (-3,08; -2,73)φ1,1 0,88 0,03 (0,82; 0,93)φ2,2 0,33 0,13 (0,06; 0,57)φ2,1 0,56 0,13 (0,31; 0,83)ρε 0,88 0,01 (0,86; 0,79)

Considerando o sumario das quantidades estimadas para o modelo 3 que incorpora o

parametro auto-regressivo φ2,1, que indica que o valor da variavel latente h2 no instante de

tempo t depende do valor de h2 no tempo t− 1 e do valor de h1 no tempo t− 1, notamos

que este modelo oferece estimativas mais precisas para os parametros auto-regressivos.

Na Figura 6, nos temos os graficos das raızes quadradas da volatilidade estimada

quando o modelo 3 e considerado.

65

(a) NE (b) CE

(c) CE (d) SW


Pode-se observar atraves da Figura 6 que, assim como no modelo 2, nao parece existir

muita diferenca entre os nıveis de volatilidade em cada um dos pares de regioes. A

diferenca mais evidente nos resultados obtidos por este modelo, quando comparado ao

modelo 2, e que aqui existe uma variacao muito maior, com a presenca de varios saltos

na volatilidade.

Modelo 4 Quando o modelo 4 e utilizado, os hiperparametros das distribuicoes a priori

de σ2ηi, i = 1, 2 sao c1 = 13, c2 = 9 e d1 = d2 = 1, quando as regioes NE e CE sao

consideradas, sao c1 = 14, c2 = 13 e d1 = d2 = 1 quando os dados das regioes CE e SW

sao utilizadas. Os hiperparametros das distribuicoes a priori das quantidades φi,i, i = 1, 2

sao os mesmos dos modelos 2 e 3 em ambos os conjuntos de dados, isto e, o produzido

pelas regioes NE e CE e o produzido pelas regioes CE e SW. Nos tambem temos que em

66

ambos os casos, os hiperparametros das distribuicoes a priori de µi sao ei = −3 e fi = 10,

i = 1, 2. O parametro ψ0 tem uma distribuicao a priori Normal N(0, 10) e ψ1 tem uma

distribuicao a priori Beta(20, 1) e os hiperparametros da distribuicao a priori de σ2ρ sao

c3 = d3 = 1 quando as regioes NE e CE ou as regioes CE e SW sao levadas em conta. Na

Tabela 4, as estimativas das quantidades de interesse sao apresentadas.


Regiao Parametro Media D.P. Intervalo de Credibilidade (95%)NE-CE τ ρ 4,13 1,20 (2,39; 7,19)

τ η1 17,05 3,61 (10,46; 24,39)τ η2 11,61 2,89 (6,65; 18,94)µ1 -3,02 0,06 (-3,13; -2,90)µ2 -2,89 0,07 (-3,02; -2,76)φ1,1 0,72 0,13 (0,38; 0,87)φ2,2 0,72 0,08 (0,55; 0,85)ψ1 0,85 0,04 (0,76; 0,92)ψ0 2,13 0,15 (1,84; 2,43)

CE-SW τ ρ 6,79 2,05 (3,16; 11,05)τ η1 15,41 2,97 (10,66; 22,13)τ η2 16,43 3,49 (10,66; 24,15)µ1 -2,79 0,05 (-2,89; -2,69)µ2 -2,84 0,06 (-2,96; -2,74)φ1,1 0,59 0,12 (0,30; 0,77)φ2,2 0,55 0,13 (0,25; 0,77)ψ1 0,93 0,03 (0,87; 0,97)ψ0 2,99 0,24 (2,49; 3,42)

A Figura 7 apresenta graficos similares aos dados pelas Figuras 4, 5 e 6, mas agora

nos utilizamos o modelo 4.

67

(a) NE (b) CE

(c) CE (d) SW


Os graficos do modelo 4 apresentam resultados bem diferentes dos obtidos pelos mo-

delos anteriores. Aqui observa-se uma pequena variacao nos nıveis de volatilidade. Porem

ainda e possıvel observar uma pequena diferenca entre as regioes NE e CE. Estas apre-

sentam valores bem similares de volatilidade, assim como foi observado nos modelos 2 e

3, porem e possıvel observar em alguns perıodos uma maior variacao na regiao CE.

Modelo 5 No caso do modelo 5, isto e, quando a causalidade de Granger esta presente

em h2(t), t = 1, 2, ..., N, nos temos que os hiperparametros das distribuicoes a priori

das quantidades φi,i, i = 1, 2 e σ2ρ sao os mesmos do modelo 4 e os hiperparametros

da distribuicao a priori de φ2,1 sao os mesmos do modelo 3 para ambos os conjuntos de

dados, isto e, o produzido pelas regioes NE e CE e o produzido pelas regioes CE e SW. Os

hiperparametros das distribuicoes a priori de σ2η1

e σ2η2

sao c1 = 17, c2 = 12 e d1 = d2 = 1

68

quando as regioes NE e CE sao consideradas e sao c1 = 15, c2 = 16 e d1 = d2 = 1

para as regioes NE e SW, ψ0 tem uma distribuicao a priori Gama cujos hiperparametros

sao (2, 1) e (3, 1) quando as regioes NE e CE e as regioes CE e SW sao consideradas,

respectivamente. A quantidade ψ1 tambem tem uma distribuicao a priori Gama com

hiperparametros (1, 1) para quando as regioes NE-CE e CE-SW sao consideradas. (Aqui

nos estamos considerando uma distribuicao Gama(α, β) com media α/β e variancia α/β2.)

A Tabela 5 traz o sumario das quantidades de interesse.


Regiao Parametro Media D.P. Intervalo de Credibilidade (95%)NE-CE τ ρ 3,03 1,02 (1,59; 5,96)

τ η1 19,99 4,34 (12,27; 30,15)τ η2 13,10 2,62 (8,77; 19,48)µ1 -3,04 0,07 (-3,18; -2,90)µ2 -2,90 0,07 (-3,03; -2,77)φ1,1 0,80 0,11 (0,44; 0,91)φ2,2 0,68 0,11 (0,40; 0,83)φ2,1 0,21 0,15 (0,02; 0,60)ψ1 0,78 0,06 (0,65; 0,89)ψ0 2,11 0,13 (1,86; 2,36)

CE-SW τ ρ 6,24 2,30 (2,32; 11,18)τ η1 16,72 3,25 (10,78; 24,01)τ η2 18,96 3,74 (12,71; 27,17)µ1 -2,84 0,07 (-2,97; -2,69)µ2 -2,88 0,07 (-3,02; -2,74)φ1,1 0,82 0,05 (0,71; 0,91)φ2,2 0,42 0,12 (0,20; 0,65)φ2,1 0,43 0,11 (0,25; 0,69)ψ1 0,91 0,04 (0,82; 0,96)ψ0 2,95 0,21 (2,54; 3,34)

Atraves dos resultados do modelo 5 nos podemos notar que as estimativas neste modelo

sao mais precisas do que as do modelo 4. Outra vantagem deste modelo para com os demais

apresentados anteriormente, e que este incorpora todos os parametros introduzidos por

estes. Isto mostra que este e um modelo mais completo e que pode trazer resultados mais

sofisticados.

Os graficos da volatilidade estimada produzida pelo modelo 5 sao dados na Figura 8.

69

(a) NE (b) CE

(c) CE (d) SW


Podemos observar atraves da Figura 8, que o modelo 5 nao apresenta valores tao

pequenos de volatilidade como o modelo 4. E novamente mostra nao existir muita dife-

renca em cada par de regioes. Apenas, mais uma vez, e possivel obervar que a regiao CE

apresenta uma maior variacao em alguns perıodos quando comparada a regiao NE.

Modelo 6 Quando o modelo 6 e levando em conta, nos temos que as quantidades φi,i,

µi, i = 1, 2 e ρε tem distribuicoes a priori com os mesmos hiperparametros do modelo

4 em ambos os casos, isto e, quando as regioes NE e CE sao consideradas e quando as

regioes CE e SW sao levadas em conta. Os hiperparametros das distribuicoes a priori de

σ2η1

e σ2η2

sao c1 = 20, c2 = 13 e d1 = d2 = 1 quando as regioes NE e CE sao consideradas

e sao c1 = 17, c2 = 19 e d1 = d2 = 1 para as regioes NE e SW. Os hiperparametros da

distribuicao a priori da quantidade λw sao o = p = 1 quando dos dados das regioes NE

70

e CE sao considerados e sao o = 1 e p = 2 quando os dados das regioes CE e SW sao

utilizados. Na Tabela 6 nos apresentamos o sumario das quantidades de interesse.


Regiao Parametro Media D.P. Intervalo de Credibilidade (95%)NE-CE λ 6,03 1,33 (3,92; 9,11)

τ η1 23,13 4,63 (14,53; 32,98)τ η2 13,71 2,95 (8,62; 19,93)µ1 -3,16 0,07 (-3,28; -3,01)µ2 -3,03 0,07 (-3,17; -2,90)φ1,1 0,78 0,08 (0,58; 0,91)φ2,2 0,71 0,11 (0,40; 0,85)ρε 0,75 0,02 (0,71; 0,78)

CE-SW λ 4,39 0,81 (3,10; 6,26)τ η1 19,41 3,72 (13,08; 27,61)τ η2 21,45 4,87 (13,54; 33,20)µ1 -2,99 0,07 (-3,13; -2,86)µ2 -3,04 0,07 (-3,17; -2,91)φ1,1 0,62 0,14 (0,29; 0,81)φ2,2 0,55 0,21 (0,06; 0,81)ρε 0,88 0,01 (0,86; 0,90)

Pela tabela acima vemos que o modelo 6 apresenta boas estimativas, porem con-

siderando as estimativas do modelo 2, que pode ser considerado um modelo similar a este,

exceto pelo fato de que o modelo 6 considera uma distribuicao t-Student para os com-

ponentes do erro, nos nao observamos um ganho consideravel com relacao a estimativas

mais precisas quando consideramos o modelo 6.

Na Figura 9, nos temos graficos similares aqueles apresentados nas Figuras 4, 5, 6, 7

e 8, mas usando o modelo 6.

71

(a) NE (b) CE

(c) CE (d) SW


O modelo 6 apresenta menores valores de volatilidade quando comparado ao modelo

2. E assim, como os demais modelos nos mostra nao existir grande diferenca em cada par

de regioes.

4.1.3 Selecao do Modelo

Como criterio de selecao de modelos foi utilizado o Deviance Information Criterion

(DIC). A Tabela 7 apresenta os valores estimados de DIC para cada modelo e cada con-

junto de dados.

72

Tabela 7: Deviance Information Criterion para os Modelos 1, 2, 3, 4, 5 e 6 quando ospares de regioes NE-CE e CE-SW sao considerados

DICRegioes Modelo 1 Modelo 2 Modelo 3 Modelo 4 Modelo 5 Modelo 6NE-CE -128,71 –745,87 -753,38 -877,87 -882,2 -740,97CE-SW 4,38 -1090,84 -1118,85 -1251,66 -1268,99 -1080,25

Considerando os valores de DIC, nos podemos notar que para ambos os pares de

regioes, o modelo 1 e muito inferior aos demais modelos. Nota-se tambem que os modelos

2 e 6, que sao similares, exceto pela distribuicao considerada para os componentes dos

erros, apresentam valores praticamente iguais de DIC. E por fim, temos que os modelos

4 e 5 se destacam com relacao aos demais, onde o modelo 5 apresenta uma valor de DIC

um pouco menor do que o do modelo 4.

4.1.4 Discussao dos Resultados Obtidos

Atraves dos resultados obtidos pelos 6 modelos pode-se observar que nao existe muita

diferenca com relacao aos valores de volatilidade em cada par de regioes estudadas. E

importante salientar tambem, que para ambos os pares de regioes temos uma volatilidade

menor e mais estavel para as medidas semanais de ozonio apos a quadragesima semana

(proximo ao ano de 1998).

Notar que nos temos um grafico mais estavel dos nıveis de volatilidade, com raras

excecoes, apos a semana 600. Esta semana corresponde a um perıodo proximo do fim do

ano 2001 e o do comeco do ano de 2002. E interessante chamar a atencao para o fato de

que no ano 2000 nos temos uma serie de medidas que foram tomadas pelas autoridades

ambientais do Mexico desde 1990 com o objetivo de reduzir os nıveis de ozonio em grandes

cidades do paıs e em particular na Cidade do Mexico. A serie de medidas podem ser

aproximadamente descritas como se segue. Desde 1990, os carros que circulam na Area

Metropolitana da Cidade do Mexico devem fazer inspecoes periodicas de suas condicoes

mecanicas. Alem disso, restricoes foram impostas na circulacao de carros de acordo com

o numero final de suas placas. Em 1997, restricoes adicionais foram implementadas em

73

termos de permitir carros ”limpos”de circular livremente. Em 1999, os fabricantes foram

encorajados a produzir carros com tecnologia moderna e limpa. A producao desses carros

foi feita de maneira compulsoria em 2001.

E possıvel observar que muitas medidas tomadas pelas autoridades ambientais pro-

duziram um efeito positivo em direcao ao objetivo de diminuir e estabilizar os nıveis das

medicoes de ozonio na Cidade do Mexico e em modificar seu comportamento a longo

prazo.

4.2 Analise dos Dados de Poluicao da Cidade de Sao

Paulo

4.2.1 Descricao dos Dados

Nesta secao, nos aplicamos os modelos de volatilidade estocastica multivariados para o

caso das medias semanais das medidas de poluicao na cidade de Sao Paulo no perıodo entre

Maio de 1996 e Dezembro de 2006. As medidas diarias dos poluentes (SO2, NO2, PM10,

O3 e CO) foram coletadas pelo Instituto de Astronomia, Geofısica e Ciencias Atmosfericas

da USP (http://www.iag.usp.br/).

Um dos grandes interesses das autoridades ligadas as polıticas ambientais da cidade

de Sao Paulo e analisar e modelar o grau de variabilidade desses diferentes poluentes e

a relacao deles com as medidas de intervencao tomadas pelas autoridades autoridades

ambientais, com o objetivo de melhorar a qualidade do ar, nos ultimos anos.

Na Figura 2, nos temos as series temporais para os poluentes (SO2, NO2, PM10, O3

e CO) para o perıodo que vai de Maio de 1996 a Dezembro de 2006. Desses graficos

da Figura 2, nos observamos que existe um decrescimo consistente nas medias semanais

durante o perıodo de 1996 a 2006 para os poluentes CO, SO2 e PM10. Nos tambem

observamos uma forma cıclica para as tres series temporais, um possıvel indicador de que

74

as medicoes de poluicao mudam entre as temporadas de verao e inverno. O poluente NO2

tambem apresenta uma tendencia decrescente ao longo do perıodo de estudo, exceto para

o perıodo que vai da 300a ate a 400a semana, onde observamos um aumento das medias

semanais do poluente. Ja para o poluente O3 nos nao observamos um claro decrescimo

para as medias semanais durante o perıodo de observacao.

E importante citar que no ano de 1996, as autoridades da cidade de Sao Paulo decre-

taram uma lei para restringir o uso de veıculos em dias uteis de acordo com o seu numero

de licenca. Apos a implementacao desta lei, nos observamos um decrescimo na concen-

tracao dos poluentes, especialmente o CO e o NO2 apos a 100a semana, isto e, em meados

de 1998. Nos tambem observamos da Figura 2, que apos a 400a semana, que representa

o perıodo proximo ao comeco do ano de 2004, que ha uma estabilizacao em um menor

nıvel dos poluentes, especialmente para o CO, NO2, SO2 e PM10. Observar que no ano

de 2004, houve a implementacao de outra lei estabelecida pelas autoridades da cidade de

Sao Paulo: a resolucao CONAMA 297/02 que estabelece o Programa de Poluicao do Ar

por motocicletas e veıculos similares.

Mesmo com a implementacao dessas leis, nos nao observamos um claro decrescimo

para os nıveis de ozonio. Uma possıvel explicacao para esse fato e que o ozonio, diferente

dos outros poluentes aqui estudados, e formado na atmosfera pela reacao quımica entre

poluentes primarios, como oxidos de nitrogenio, e componentes da atmosfera, onde a

simples reducao na emissao de alguns poluentes nao e suficiente para reduzir os nıveis de

emissao de ozonio.

Na Figura 10, nos temos os graficos dos log-retornos yi(t), t = 1, 2, ..., N ; i = 1, 2, 3, 4, 5,

centrados em suas respectivas medias (de Maio de 1996 a Dezembro de 2006).

75

(a) CO (b) NO2

(c) PM10 (d) SO2

(e) O3

Figura 10: Log-retornos centrados em suas medias para os poluentes de Maio de 1996 aDezembro de 2006.

4.2.2 Uso dos Modelos e Resultados

Modelo 1 Como foi citado na secao 3.4, devido ao grande custo computacional encon-

trado ao utilizar os modelos de volatilidade estocastica multivariados, foram aplicados

apenas os modelos 1 e 2 para os dados de poluicao do ar da cidade de Sao Paulo.

76

Para obter a convergencia do algoritmo foram assumidas as distribuicoes a priori

N(ai, bi) e N(gi, hi) para φ e ρ, respectivamente. Assim sendo, assumindo o modelo

1, os hiperparametros das distribuicoes a priori de φ,σ2η e µ, onde φ = (φ11, φ22, ..., φkk) ,

σ2η =

(σ2ηi, σ2

η2, ..., σ2

ηk

)e µ = (µ1, µ2, ..., µk), sao aii = 0; bii = 10; ci = di = 1; ei = 0 e

fi = 10, i = 1, 2, ..., k. A escolha desses hiperparametros foi feita com o objetivo de termos

distribuicoes a priori nao-informativas e tambem que possibilitassem a convergencia do

algoritmo. Na Tabela 8, temos o sumario das quantidades estimadas quando o modelo 1

e considerado.

Tabela 8: Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de credibili-dade (95%) para as quantidades de interesse quando o modelo 1 e utilizado.

Poluente Parametro Media D. P. Intevalo de credibilidade 95%µ1 -2,17 0,08 (-2,33; -2,03)

O3 φ1 0,41 0,20 (0,01; 0,76)τ η1 4,97 1,47 (2,71; 8,50)µ2 -2,78 0,11 (-3,00; -2,58)

CO φ2 0,70 0,10 (0,46; 0,85)τ η2 4,19 1,30 (2,13; 7,26)µ3 -2,71 0,09 (-2,89; -2,55)

SO2 φ3 0,46 0,13 (0,19; 0,68)τ η3 2,90 0,79 (1,67; 4,80)µ4 -2,65 0,10 (-2,84; -2,46)

PM10 φ4 0,66 0,13 (0,36; 0,84)τ η4 4,36 1,39 (2,08; 7,57)µ5 -2,81 0,08 (-2,96; -2,66)

NO2 φ5 0,31 0,17 (-0,03; 0,61)τ η5 4,26 1,34 (2,29; 7,37)

Na Figura 11 temos os graficos das estimativas das raızes quadradas da volatilidade

quando o modelo 1 e utilizado.

77

(a) CO (b) NO2

(c) PM10 (d) SO2

(e) O3

Figura 11: Raızes quadradas das volatilidades para os poluentes quando o modelo 1 eutilizado.

Observando os graficos da Figura 11, nota-se que os valores de volatilidade encontrados

para o NO2 e do O3 apresentam uma variacao bem menor quando comparados aos obtidos

para os demais poluentes. Tambem podemos observar que, para os demais poluentes, nao

existe uma estabilizacao dos nıveis de volatilidade ao longo do perıodo estudado.

78

Modelo 2 No caso do modelo 2, os hiperparametros das distribuicoes a priori de φ,σ2η,µ

e ρ, onde ρ =(ρ12, ρ13, ..., ρk−1,k

), sao (a1 = 0, 41; a2 = 0, 70; a3 = 0, 46; a4 = 0, 66;

a5 = 0, 31); (b1 = 25, 15; b2 = 103, 53; b3 = 63, 09; b4 = 62, 39; b5 = 35, 86); (c1 = 11, 42;

c2 = 10, 38; c3 = 13, 38; c4 = 9, 79; c5 = 10, 18); (d1 = 2, 30; d2 = 2, 48; d3 = 4, 61;

d4 = 2, 24; d5 = 2, 39); (e1 = −2, 17; e2 = −2, 78; e3 = −2, 71; e4 = −2, 65; e5 = −2, 81);

(f1 = 167, 79; f2 = 88, 83; f3 = 129, 13; f4 = 107, 01; f5 = 173, 36); gi = 0 e hi = 10. E

importante chamar atencao para o fato que a mudanca nos hiperparametros de φ,σ2η e µ

e porque foram utilizadas as informacoes providas pelos resultados do modelo 1. Portanto,

estamos usando uma aproximacao Bayesiana empırica para analisar o problema (CARLIN

e LOUIS, 2000). Na Tabela 9 apresentamos o sumario das quantidades estimadas.

Tabela 9: Estimativa da media a posteriori, desvio-padrao (D.P.) e intervalo de credibili-dade (95%) para as quantidades de interesse quando o modelo 2 e utilizado.

Poluente Parametro Media D. P. Intevalo de credibilidade 95%µ1 -2,18 0,06 (-2,30; -2,06)

O3 φ1 0,56 0,10 (0,36; 0,74)τ η1 1,55 0,26 (1,12; 2,16)µ2 -2,76 0,07 (-2,90; -2,64)

CO φ2 0,66 0,07 (0,52; 0,79)τ η2 5,48 1,15 (3,57; 7,92)µ3 -2,72 0,07 (-2,85; -2,59)

SO2 φ3 0,53 0,08 (0,37; 0,67)τ η3 1,31 0,19 (0,98; 1,76)µ4 -2,63 0,07 (-2,77; -2,50)

PM10 φ4 0,62 0,11 (0,39; 0,81)τ η4 2,01 0,40 (1,40; 2,94)µ5 -2,82 0,05 (-2,92; -2,72)

NO2 φ5 0,40 0,12 (0,16; 0,62)τ η5 5,58 1,20 (3,56; 8,26)ρ12 0,24 0,04 (0,15; 0,32)ρ13 0,59 0,03 (0,53; 0,64)ρ14 0,76 0,02 (0,72; 0,80)ρ15 0,69 0,02 (0,64; 0,73)ρ23 0,53 0,03 (0,46; 0,59)ρ24 0,66 0,03 (0,60; 0,71)ρ25 0,54 0,03 (0,47; 0,60)ρ34 0,79 0,02 (0,75; 0,83)ρ35 0,68 0,02 (0,63; 0,73)ρ45 0,72 0,02 (0,68; 0,76)

79

Atraves dos resultados da Tabela 9 nota-se que o modelo 2 aparenta estimar melhor a

volatilidade do que o modelo 1. Ele apresenta estimativas mais precisas do que as obtidas

pelo modelo 1. Alem disso ele incorpora os parametros de covariancia dos componentes

dos erros.

A Figura 12 apresenta os graficos das estimativas das raızes quadradas da volatilidade

quando o modelo 2 e utilizado.

(a) CO (b) NO2

(c) PM10 (d) SO2

(e) O3

Figura 12: Raızes quadradas das volatilidades para os poluentes quando o modelo 2 eutilizado.

80

Observando os graficos obtidos pelo modelo 2, nota-se que os valores de volatilidade

para o NO2 sao semelhantes aos obtidos pelo modelo 1. Porem para os demais poluentes

observa-se uma grande diferenca nos resultados obtidos, principalmente para o ozonio, que

agora, apresenta altos valores de volatilidade. O que se pode observar com base nestes

resultados, exceto para o NO2, que os nıveis de volatilidade estao muito elevados, e nao

apresentam uma tendencia de estarem diminuindo ou estabilizando.

4.2.3 Selecao do Modelo

O criterio DIC foi utilizado para discriminar os dois modelos propostos. Os valores de

DIC assumindo os modelos 1 e 2 sao apresentados na Tabela 10.

Tabela 10: Deviance Information Criterion para os Modelos 1 e 2.

DICModelo 1 Modelo 2

928,30 -873,54

Observando os resultados da Tabela 10, nos verificamos que o modelo 2, atraves do

criterio DIC, seria o melhor modelo para descrever a volatilidade da serie de dados de

poluicao da cidade de Sao Paulo.

4.2.4 Discussao dos Resultados Obtidos

E importante observar atraves dos resultados obtidos pelos dois modelos aplicados

que nao ha uma evidente diminuicao e estabilizacao dos nıveis de volatilidade para os

poluentes SO2, PM10, O3 e CO ao longo do perıodo de estudo. Apenas para o poluente

NO2 foi possıvel observar uma estabilizacao nos nıveis de volatilidade.

Como citado anteriormente, em 1996 , foi decretada uma lei para restringir o uso

de veıculos em dias uteis de acordo com o seu numero de licenca. E em 2004, houve a

81

implementacao de outra lei, a resolucao CONAMA 297/02, que estabeleceu o Programa

de Poluicao do Ar por motocicletas e veıculos similares.

Portanto, esses resultados aqui obtidos demonstram que as leis implementadas pelas

autoridades ambientais da cidade de Sao Paulo, aparentemente, nao tem resultado em

uma diminuicao dos nıveis de volatilidade, principalmente para o SO2, PM10, O3 e CO

. Isso e um indıcio de que as medidas que estao sendo tomadas nao tem surtido grande

efeito, e portanto e necessario tentar implementar novas leis para se obter resultados mais

significativos.

A Figura 13 apresenta os graficos de volatilidade dos cinco poluentes, obtidos pelo mo-

delo 2, para o perıodo que vai de Janeiro de 1998 a Dezembro de 2006 e os graficos contendo

o numero de internacoes por doencas do aparelho circulatorio e respiratorio. Os dados de

internacoes foram obtidos atraves do DATASUS (http://www2.datasus.gov.br/DATASUS/

index.php). Este sistema possui informacoes sobre essas internacoes a partir de 1998.

82

(a) CO (b) NO2 (c) PM10

(d) SO2 (e) O3

(f) Aparelho circulatorio (g) Aparelho respiratorio

Figura 13: Raızes quadradas das volatilidades para os poluentes quando o modelo 2 eutilizado e numero de internacoes por doencas do aparelho circulatorio e respiratorio.

Pode-se observar atraves da Figura 13 que o numero de internacoes por doencas do

aparelho respiratorio e circulatorio tende a aumentar no perıodo observado, apresentando

grandes saltos a partir de 2002. Nota-se tambem que a volatilidade tambem tende a apre-

sentar maiores valores a partir desse ano, indicando assim que, possivelmente, a grande

variacao encontrada para os poluentes no perıodo aqui estudado tem contribuido para o

aumento do numero de internacoes para essas doencas. E importante citar que o numero

83

de internacoes apresentado representa o total de internacoes em toda a cidade de Sao

Paulo, enquanto as medidas de poluicao foram obtidas em apenas uma regiao de Sao

Paulo. Portanto essa relacao observada entre a volatilidade e o numero de internacoes e

um pouco subjetiva, ja que ambos foram coletados em proporcoes diferentes.

84

5 Conclusoes e Perspectivas Futuras

A aplicacao dos modelos de volatilidade estocastica nos conjuntos de dados apresen-

tados se mostrou uma otima ferramenta para estudar a volatilidade em dados de poluicao

do ar. Os resultados obtidos nos permitiram observar com clareza o comportamento da

volatilidade ao longo dos perıodos de estudo, permitindo assim avaliar se as leis imple-

mentadas para o combate a poluicao estao tendo um efeito positivo.

Foi possıvel observar que a medida que aplicamos modelos mais complexos, com in-

clusao de mais parametros, as estimativas obtidas tendem a ser melhores do que as dos

modelos com menos parametros.

A analise dos dados da Cidade do Mexico nos mostrou que as medidas tomadas pelo

governo para diminuir os problemas causados pelo ozonio tem sido eficaz e isso pode

ser observado com clareza com as estimativas obtidas da volatilidade. A cidade de Sao

Paulo tambem apresentou um decrescimo dos nıveis de volatilidade dos poluentes apos

a implementacao de algumas leis, porem nao com a mesma clareza observada na Cidade

do Mexico. Indicando assim que novas solucoes para os problemas de poluicao devem ser

buscadas.

E importante ressaltar que existem outros modelos que podem ser aplicados para

analisar este tipo de dados. Alem dos modelos multivariados aqui apresentados, pode-se

aplicar, por exemplo, os modelos de volatilidade estocastica com a presenca de outro efeito

aleatorio, dado por:

y (t) = σ (t) ε (t)

onde σ (t) = exp

12h (t) + ω (t)

e ω (t) e um efeito aleatorio com distribuicao Normal

N (0, σ2ω).

Esses modelos, ao inves de utilizar a estrutura multivariada aqui apresentada, fazem

uso de um efeito aleatorio para capturar a correlacao entre as series de dados. Sera

85

importante aplicar esses modelos e compara-los aos modelos multivariados, com o objetivo

de verificar qual apresenta melhores estimativas da volatilidade.

Tambem e interessante estudar mais sobre a relacao entre o numero de internacoes e

a volatilidade. Alem da observacao do comportamento desses dados ao longo do tempo,

e importante incorporar as informacoes do numero de internacoes nos modelos de volati-

lidade, com o objetivo de termos resultados mais consistentes.

Alem de novos modelos, e interessante aplicar outros criterios de selecao de modelos,

como o criterio BIC (SCHWARZ, 1978) e o PCO (GELFAND et al., 1992).

86

Referencias Bibliograficas

ACHCAR, J. A., FERNANDEZ-BREMAUNTZ, A. A., RODRIGUES, E. R. e TZIN-

TZUN, G. Estimating the number of ozone peaks in Mexico City using a non-homogeneous

Poisson model. Environmetrics, v.19, p.469-485, 2007.

ACHCAR, J. A.; ZOZOLOTTO, H. C. e RODRIGUES, E. R. Bivariate stochastic

volatility models applied to Mexico City ozone pollution data. Em: Gaia C.

Romano e Alice G. Conti. (Org.). Air Quality in the 21st Century, 1a ed., Nova Iorque:

Nova Publishers, v.01, p.0001-0030, 2008a.

ACHCAR, J. A.; RODRIGUES, E. R. e ZOZOLOTTO, H. C. Bivariate stochastic vo-

latility models applied to air pollution data. Revista Brasileira de Biometria, v.26,

p.67-81, 2008b.

ALVAREZ, L. J., FERNANDEZ-BREMAUNTZ, A. A., RODRIGUES, E. R. e TZIN-

TZUN, G. Maximum a posteriori estimation of the daily ozone peaks in Mexico City.

Journal of Agricultural, Biological, and Environmental Statistics, v.10, p.276-

290, 2005.

AUSTIN, J. e TRAN, H. A characterization of the weekday-weekend behavior

of ambient ozone concentrations in California. Em: Air Pollution VII. UK: WIT

Press, p.645-661, 1999.

BELL, M. L., McDERMONTT, A., ZEGER, S. L., SAMET, J. M. e DOMINICI, F.

Ozone and short-term mortality in 95 U.S. urban communities. Journal of the Ameri-

can Medical Society, v.292, p.2372-2378, 2004.

BELL, M. L., PENG, R. e DOMINICI, F. The exposure-response curve for ozone and

risk of mortality and the adequacy of current ozone regulations. Environmental Health

Perspectives, v.114, p.532-536, 2005.

87

BELL, M. L., GOLDBERG, R., ROGREFE, C., KINNEY, P. L., KNOWLTON, K.,

LYNN, B., ROSENTHAL, J., ROSENZWEIG, C. e PATZ, J. A. Climate change, ambi-

ent ozone, and health in 50 US cities. Climate change, v.82, p.61-76, 2007.

BERG, A., MEYER R. e YU, J. Deviance information criterion for computing volatility

models. Journal of Business and Economic Statistics, v.22, p.107-120, 2004.

BERNARDO, J. M. e SMITH, A. F. M. Bayesian theory. New York: John Wiley &

Sons, 1995.

BERHANE, K. e MOLITOR, N-T. A Bayesian approach to functional-based multilevel

modeling of longitudinal data: applications to environmental epidemiology. Biostatis-

tics, v.9, n.4, p.686-699, 2008.

BESAG, J. E. Discussion on the paper by Grenander and Miller. Journal of the Royal

Statistical Society, v.56, p.591-592, 1994.

BOLLERSLEV, T. Generalized autoregressive conditional heterocedasticity. Journal of

Econometrics, v.31, p.307-327, 1986.

BRAGA, A. L. F., ZANOBETTI, A. e SCHWARTZ, J. The effect of weather on respira-

tory and cardiovascular deaths in 12 U.S. cities. Enviromental Health Perspectives,

v.110, p.859-863, 2002.

BURSTYN, I., CHERRY, N. M., YASUI, Y e KIM, H. Relative performance of different

exposure modeling approaches for sulfur dioxide concentrations in the air in rural western

Canada. BMC Medical Research Methodology, v.8, n.43, 2008.

CARLIN, B. P. e LOUIS, T. A. Bayes and empirical Bayes methods for data ana-

lysis. 2a ed., Boca Raton: Chapman and Hall Press, 2000.

CENDON, S., PEREIRA, L. A. A., BRAGA, A. L. F., CONCEICAO, G. M. S., JUNIOR,

A. C., ROMALDINI, H., LOES, A. C. e SALDIVA, P. H. N. Efeitos da poluicao do ar no

88

infarto do miocardio. Rev. Saude Publica, v.40, p.414-419, 2006.

CONCEICAO, G. M. S., SALDIVA, P. H. S. e SINGER, J. M. Modelos MLG e e MAG

para analise de associacao entre poluicao atmosferica e marcadores de morbi-mortalidade:

uma introducao baseada em dados da cidade de Sao Paulo. Revista Brasileira de Epi-

demiologia, v.4, n.3, p.206-219, 2001.

DANIELSSON, J. Stochastic volatility in asset prices: estimation with simulated maxi-

mum likelihood. Journal of Econometrics, v.64, p.375-400, 1994.

DOMINIC, F., WANG, C., CRAINICEANU, C. e PARMIGIANI, G. Model selection

and health effect estimation in environmental epidemiology. Epidemiology, v.19, n.4,

p.558-560, 2008.

ENGLE, R. F. Autoregressive conditional heteroscedasticity with estimates of the vari-

ance of United Kingdon inflation. Econometrica, v.50, p.987-1007, 1982.

FARHAT, S. C. L., PAULO, R. L. P., SHIMODA, L. M., CONCEICAO, G. M. S., LIN,

C. A., BRAGA, A. L. F., WARTH, M. P. N. e SALDIVA, P. H. N. Effect of air pollu-

tion on pediatric respiratory emergency room visits and hospital admissions. Brazilian

Journal of Medical and Biological Research, v.38, p.227-235, 2005.

FREITAS, C., BREMNER, S. A., GOUVEIA, N., PEREIRA, L. A. A., SALDIVA, P. H.

N. Internacoes e obitos e sua relacao com a poluicao atmosferica em Sao Paulo, 1993 a

1997. Rev. Saude Publica, v.38, p.751-757, 2004.

GAMERMAN, D. Markov Chain Monte Carlo: stochastic simulation for Baye-

sian inference. 1a ed., Londres: Chapman & Hall, 245p, 1997.

GELFAND, A. E., DEY, D. K. e CHANG, H. Model determination using predictive

distributions with implementation sampling-based method (with discussion).

In Bayesian Statistics 4, BERNARDO, J. M., BERGER, J. O., DAWID, A. P. e SMITH,

A. F. P. , editores. Oxford: Oxford University Pres, 1992.

89

GELFAND, A. E. e SMITH, A. F. M. Sampling-based approaches to calculating marginal

densities. Journal of the American Statistical Association, v.85, p.398-409, 1990.

GEMAN, S. e GEMAN, D. Stochastic relaxation, Gibbs distribution and the Bayesian

restoration of images. IEEE Transactions on Pattern Analysis and Machine In-

telligence, v.6, p.721-741, 1984.

GEWEKE, J. Evaluating the accuracy of sampling-based approaches to the cal-

culation of posterior moments Em: J.O. Berger, J.M. Bernardo, A.P. Dawid, e A.F.M.

Smith. (Eds.). Bayesian Statistics 4, Oxford: Oxford University Press, p.169-194, 1992.

GHYSELS, E., HARVEY, A. C. e RENAULT, E. Stochastic volatility. Em: RAO, C.

R. e MADALLA, G. S. (Eds.). Statistical Models in Finance, Amsterdam:North-Holland,

1996.

GELMAN, A. e RUBIN, B. D. Inference from iterative simulation using multiple se-

quences. Statistical Science, v.7, p.457-511, 1992.

GOUVEIA, N., FREITAS, C. U., MARTINS, L. C. e MARCILIO, I. O. Hospitaliza-

coes por causas respiratorias associadas a contaminacao atmosferica no municıpio de Sao

Paulo, Brasil. Cad. Saude Publica, Rio de Janeiro, v.22, p.2669-2677, 2006.

GUARDANI, R., NASCIMENTO, C. A. O., GUARDANI, M. L. G., MARTINS, M. H.

R. B. e ROMANO, J. Study of atmospheric ozone formation by means of a neural network

based model. J. Air and Waste Management Assoc., v.49, p.316-323, 1999.

GUARDANI, R., AGUIAR, J. L., NASCIMENTO, C. A. O., LACAVA, C. I. V. e

YANAGI, Y. Ground-level ozone mapping in large urban areas using multivariate analy-

sis: application to the Sao Paulo Metropolitan Area. J. Air and Waste Management

Assoc., v.53, p.553-559, 2003.

HOROWITZ, J. Extreme values from a nonstationary stochastic process: an application

to air quality analysis. Technometrics, v.22, p.469-482, 1980.

90

JAVITS, J. S. Statistical interdependencies in the ozone national ambient air quality stan-

dard. Journal of the Air Pollution Control Association, v.30, p.58-59, 1980.

KIM, S., SHEPARD, N. e CHIB, S. Stochastic volatility: likelihood inference and com-

parison with ARCH models. Review of Economic Studies, v.65, p.361-393, 1998.

LANFREDI, M. e MACCHIATO, M. Searching for dimensionality in air pollution time

series. Europhysics Lett., v.40, p.589-594, 1997.

LARSEN, L. C., BRADLEY, R. A. e HONCOOP, G. L. A new method of characte-

rizing the variability of air quality-related indicators. Em: Air and Waste Mana-

gement Association’s International Speciality Conference of Tropospheric Ozone and the

Environment. California. USA, 1990.

LEADBETTER, M. R. On a basis for ”peak over threshold”modeling. Statistics and

Probability Letters, v.12, p.357-362, 1991.

LIN, S., BELL, E. M., LIU, W., WALKER, R. J., KIM, N. K. e HWANG, S-A. Ambient

ozone concentration and hospital admissions due to childhood respiratory diseases in New

York State, 1991– 2001. Environmental Research, v.108, p.42-47, 2008.

LOOMIS, D. P., BORJA-ARBUTO, V. H., BANGDIWALLA, S. I. e SHY, C. M. Ozone

exposure and daily mortality in Mexico City: a time series analysis. Health Effects

Institute Research Report, v.75, p.1-46, 1996.

MARCILIO, I., GOUVEIA, N. Quantificacao do impacto da poluicao atmosferica sobre a

populacao urbana brasileira. Cad. Saude Publica, Rio de Janeiro, v.23, p.529-536,

2007.

MARTINS, L. C., PEREIRA, L. A. A., LIN, C. A., SANTOS, U. P., PRIOLI, G., LUIZ,

O. C., SALDIVA, P. H. N. e BRAGA, A. L. F. The effects of air pollution on cardiovas-

cular diseases: lag structures. Rev. Saude Publica, v.40, n.4, p.677-683, 2006.

91

MARTINS, L. C., LATORRE, M. R. D. O., CARDOSO, M. R. A., GONCALVES, F.

L. T., SALDIVA, P. H. N. e BRAGA, A. L. F. Poluicao atmosferica e atendimentos por

pneumonia e gripe em Sao Paulo, Brasil. Rev. Saude Publica, v.36, p.88-94, 2002.

MEDEIROS, A. e GOUVEIA, N. Relacao entre baixo peso ao nascer e a poluicao do ar

no municıpio de Sao Paulo. Rev. Saude Publica, v.39, n.6, p.965-972, 2005.

MEYER, R. e YU, J. BUGS for a Bayesian analysis of stochastic volatility models.

Econometrics Journal, v.3, p.198-215, 2000.

MORETTIN, P. A. Econometria financeira. Sao Paulo: Editora Edgard Blucher, v.1,

340p, 2008.

MOURA, M., JUNGER, W. L., MENDONCA, G. A. S. e DE LEON, A. P. Qualidade

do ar e transtornos respiratorios agudos em criancas. Rev. Saude Publica, v.42, n.3,

p.503-511, 2008.

NASCIMENTO, L. F. C., PEREIRA, L. A. A., BRAGA, A. L. F., MODOLO, M. C. C.

e JR, J. A. C. Efeitos da poluicao atmosferica na saude infantil em Sao Jose dos Campos,

SP. Rev. Saude Publica, v.40, n.1, p.77-82, 2006.

PAN, J.-N. e CHEN, S.-T. Monitoring long-memory air quality data using ARFIMA

model. Environmetrics, (http://www.interscience.wiley.com/10.1002/env.882).

PAULINO, C., TURKMAN, M. e MURTEIRA, B. Estatıstica Bayesiana. Lisboa, Por-

tugal: Fundacao Calouste Gulbenkian, 446 p., 2002.

RAFTERY, A. E. Are ozone exceedance rate decreasing?, Comment of the paper ”Ex-

treme value analysis of environmental time series: an application to trend detection in

ground-level ozone” por R. L. Smith. Statistical Sciences, v.4, p.378-381, 1989.

ROBERTS, E. M. Review of statistics extreme values with applications to air quality

data. Part I. Review. Journal of the Air Pollution Control Association, v.29,

92

p.632-637, 1979a.

ROBERTS, E. M. Review of statistics extreme values with applications to air quality

data. Part II. Applications. Journal of the Air Pollution Control Association,

v.29, p.733-740, 1979b.

SCHWARZ, G. Estimating the dimension of a model. The Annals of Statistics, v.6,

p.461-464, 1978.

SMITH, R. L. Extreme values analysis of environmental time series: an application to

trend detection in ground-level ozone. Statistical Sciences, v.4, p.367-393, 1989.

SMITH, A. F. M. e ROBERTS, G. O. Bayesian Computation via the Gibbs samples and

related Markov Chain Monte Carlo methods (with discussion). Journal of the Royal

Statistical Society, serie B, v.55, p.3-23, 1993.

SPIEGELHALTER, D. J., BEST, N. G., CARLIN, B. P. e VAN DER LINDE, A. Baye-

sian measures of model complexity and fit. Journal of the Royal Statistical Society,

serie B, v.64, p.583-639, 2002.

SPIEGELHALTER, D. J., THOMAS, A., BEST, N. G. e LUNN, D. Winbugs user

manual. Cambridge, United Kingdon: MRC Biostatistics Unit, 1999.

WELTY, L. J., PENG, R. D., ZEGER, S. L. e DOMINIC, F. Bayesian distributed lag

models: estimating effects of particulate matter air pollution on daily mortality. Bio-

metrics, v.65, p.282-291, 2009.

WILSON, R., COLOME, S. D., SPENGLER, J. D. e WILSON, D. G. Health effects in

fossil fuel burning: assessment and mitigation. Cambridge, USA: Ballinger, 1980.

YU, J. Forecasting volatility in the New Zeland stock market. Applied Financial Eco-

nomics, v.12, p.193-202, 2002.

93

YU, J. e MEYER, R. Multivariate stochastic volatility models: Bayesian estimation and

model comparison. Econometric Reviews, v.25, p.361-384, 2006.

ZOLGHADRI, A. e HENRY, D. Minmax statistical models for air pollution time series.

Application to ozone time series data measured in Bordeaux. Environmental Moni-

toring and Assessment , v.98, p.275-294, 2004.

94

A Programas

Esse apendice apresenta alguns programas computacionais utilizados nesta dissertacao

de mestrado. Os programas utilizados em cada aplicacao desta dissertacao estao apresen-

tados a seguir.

A.1 Cidade do Mexico

A.1.1 Modelo 1

Listagem 1: Programa Desenvolvido no software Winbugs (Modelo 1).

1

2 model

3 4 f o r ( i in 1 :N)

5 6 Y1var [ i ] <− exp ( th1 [ i ] )

7 Y1tau [ i ] <− 1/Y1var [ i ]

8 Y1 [ i ] ˜ dnorm (0 , Y1tau [ i ] )

9 SdY1 [ i ] <− s q r t ( Y1var [ i ] )

10 Y2var [ i ] <− exp ( th2 [ i ] )

11 Y2tau [ i ] <− 1/Y2var [ i ]

12 Y2 [ i ] ˜ dnorm (0 , Y2tau [ i ] )

13 SdY2 [ i ] <− s q r t ( Y2var [ i ] )

14 15 thm1 [ 1 ] <− mu1

16 thm2 [ 1 ] <− mu2

17 th1 [ 1 ] ˜ dnorm(thm1 [ 1 ] , sigma . eta1 )


19 f o r ( i in 2 :N)

20 21 thm1 [ i ] <− mu1 + phi1 *( th1 [ i −1] − mu1)

22 th1 [ i ] ˜ dnorm(thm1 [ i ] , sigma . eta1 )

23 thm2 [ i ] <− mu2 + phi2 *( th2 [ i −1] − mu2)

95


25 26 phi1 ˜ dbeta (1 , 1 )

27 sigma . eta1 ˜ dgamma(1 , 1 )

28 mu1 ˜ dnorm ( 0 , 0 . 0 1 )

29 phi2 ˜ dbeta (1 , 1 )


31 mu2 ˜ dnorm ( 0 , 0 . 0 1 )

32

A.1.2 Modelo 2


1

2 model

3 4 f o r ( i in 1 :N)

5 6 ysdet [ i ] <− exp ( th [ i , 1 ] + th [ i , 2 ] ) * ( 1 − rhoep* rhoep )

7 Yis2 [ i , 1 , 1 ] <− exp ( th [ i , 2 ] ) / ysdet [ i ]


9 Yis2 [ i , 1 , 2 ] <− −rhoep*exp (0 . 5* th [ i , 1 ] + 0 .5* th [ i , 2 ] ) / ysdet [ i ]

10 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]

11 Y[ i , 1 : 2 ] ˜ dmnorm(muy [ ] , Yis2 [ i , , ] )

12 Yis2a [ i ] <− exp(−th [ i , 1 ] )

13 Yis2b [ i ] <− exp(−th [ i , 2 ] )

14 sdY1 [ i ] <− 1/ s q r t ( Yis2a [ i ] )

15 sdY2 [ i ] <− 1/ s q r t ( Yis2b [ i ] )

16 17 muy [ 1 ] <− 0

18 muy [ 2 ] <− 0

19 thm [ 1 , 1 ] <− mu1

20 thm [ 1 , 2 ] <− mu2

21 th [ 1 , 1 ] ˜ dnorm(thm [ 1 , 1 ] , i taua2 )

22 th [ 1 , 2 ] ˜ dnorm(thm [ 1 , 2 ] , i taub2 )

23 f o r ( i in 2 :N)

24

96

25 thm [ i , 1 ] <− mu1 + phi1 *( th [ i −1 ,1] − mu1)

26 thm [ i , 2 ] <− mu2 + phi2 *( th [ i −1 ,2] − mu2)

27 th [ i , 1 ] ˜ dnorm(thm [ i , 1 ] , i taua2 )

28 th [ i , 2 ] ˜ dnorm(thm [ i , 2 ] , i taub2 )

29 30 phi1 ˜ dbeta (1 , 1 )

31 phi2 ˜ dbeta (1 , 1 )

32 i t aua2 ˜ dgamma(5 , 1 )

33 i taub2 ˜ dgamma(4 , 1 )

34 mu1 ˜ dnorm(−3 ,0.01)

35 mu2 ˜ dnorm(−3 ,0.01)

36 rhoep ˜ dun i f (−1 ,1)

37

A.1.3 Modelo 3


1

2 model

3 4 f o r ( i in 1 :N)

5 6 ysdet [ i ] <− exp ( th [ i , 1 ] + th [ i , 2 ] ) * ( 1 − rhoep* rhoep )




10 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]


12 Yis2a [ i ] <− exp(−th [ i , 1 ] )

13 Yis2b [ i ] <− exp(−th [ i , 2 ] )

14 sdY1 [ i ] <− 1/ s q r t ( Yis2a [ i ] )

15 sdY2 [ i ] <− 1/ s q r t ( Yis2b [ i ] )

16 17 muy [ 1 ] <− 0

18 muy [ 2 ] <− 0

19 thm [ 1 , 1 ] <− mu1

20 thm [ 1 , 2 ] <− mu2

97



23 f o r ( i in 2 :N)

24 25 thm [ i , 1 ] <− mu1 + phi11 *( th [ i −1,1]−mu1)

26 thm [ i , 2 ] <− mu2 + phi21 *( th [ i −1,1]−mu1) + phi22 *( th [ i −1,2]−mu2)



29 30 phi11 ˜ dbeta (1 , 1 )

31 phi21 ˜ dbeta (1 , 1 )

32 phi22 ˜ dbeta (1 , 1 )

33 i t aua2 ˜ dgamma(9 , 1 )


35 mu1 ˜ dnorm(−3 ,0.01)

36 mu2 ˜ dnorm(−3 ,0.01)

37 rhoep ˜ dun i f (−1 ,1)

38

A.1.4 Modelo 4


1

2 model

3 4 f o r ( i in 1 :N)

5 6 ysdet [ i ] <− exp ( th [ i , 1 ] + th [ i , 2 ] ) * ( 1 − rhoep [ i ]* rhoep [ i ] )



9 Yis2 [ i , 1 , 2 ] <− −rhoep [ i ]* exp (0 . 5* th [ i , 1 ]+0 .5* th [ i , 2 ] ) / ysdet [ i ]

10 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]


12 Yis2a [ i ] <− exp(−th [ i , 1 ] )

13 Yis2b [ i ] <− exp(−th [ i , 2 ] )

14 sdY1 [ i ] <− 1/ s q r t ( Yis2a [ i ] )

15 sdY2 [ i ] <− 1/ s q r t ( Yis2b [ i ] )

98

16 17 muy [ 1 ] <− 0

18 muy [ 2 ] <− 0

19 thm [ 1 , 1 ] <− mu1

20 thm [ 1 , 2 ] <− mu2



23 q [ 1 ] ˜ dnorm( psi0 , i t au2 )

24 rhoep [ 1 ] <− ( exp ( q [ 1 ] ) − 1)/( exp ( q [ 1 ] ) + 1)

25 f o r ( i in 2 :N)

26 27 thm [ i , 1 ] <− mu1 + phi1 *( th [ i −1 ,1] − mu1)

28 thm [ i , 2 ] <− mu2 + phi2 *( th [ i −1 ,2] − mu2)



31 qmean [ i ] <− ps i 0 + p s i *( q [ i −1] − ps i 0 )

32 q [ i ] ˜ dnorm(qmean [ i ] , i t au2 )

33 rhoep [ i ] <− ( exp ( q [ i ] ) − 1)/( exp ( q [ i ] ) + 1)

34 35 phi1 ˜ dbeta (1 , 1 )

36 phi2 ˜ dbeta (1 , 1 )

37 i t aua2 ˜ dgamma(13 ,1 )


39 mu1 ˜ dnorm(−3 ,0.1)

40 mu2 ˜ dnorm(−3 ,0.1)

41 ps i 0 ˜ dnorm ( 0 , 0 . 1 )

42 p s i ˜ dbeta (20 ,1 )

43 i t au2 ˜ dgamma(1 , 1 )

44

A.1.5 Modelo 5


1

2 model

3 4 f o r ( i in 1 :N)

99

5 6 ysdet [ i ] <− exp ( th [ i , 1 ] + th [ i , 2 ] ) * ( 1 − rhoep [ i ]* rhoep [ i ] )



9 Yis2 [ i , 1 , 2 ] <− −rhoep [ i ]* exp (0 . 5* th [ i , 1 ]+0 .5* th [ i , 2 ] ) / ysdet [ i ]

10 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]


12 Yis2a [ i ] <− exp(−th [ i , 1 ] )

13 Yis2b [ i ] <− exp(−th [ i , 2 ] )

14 sdY1 [ i ] <− 1/ s q r t ( Yis2a [ i ] )

15 sdY2 [ i ] <− 1/ s q r t ( Yis2b [ i ] )

16 17 muy [ 1 ] <− 0

18 muy [ 2 ] <− 0

19 thm [ 1 , 1 ] <− mu1

20 thm [ 1 , 2 ] <− mu2



23 q [ 1 ] ˜ dnorm( ps i0 , i t au2 )

24 rhoep [ 1 ] <− ( exp ( q [ 1 ] ) − 1)/( exp ( q [ 1 ] ) + 1)

25 f o r ( i in 2 :N)

26 27 thm [ i , 1 ] <− mu1 + phi11 *( th [ i −1,1]−mu1)

28 thm [ i , 2 ] <− mu2 + phi21 *( th [ i −1,1]−mu1) + phi22 *( th [ i −1,2]−mu2)



31 qmean [ i ] <− ps i 0 + p s i *( q [ i −1] − ps i 0 )

32 q [ i ] ˜ dnorm(qmean [ i ] , i t au2 )

33 rhoep [ i ] <− ( exp ( q [ i ] ) − 1)/( exp ( q [ i ] ) + 1)

34 35 phi11 ˜ dbeta (1 , 1 )

36 phi21 ˜ dbeta (1 , 1 )

37 phi22 ˜ dbeta (1 , 1 )

38 i t aua2 ˜ dgamma(17 ,1 )

39 i taub2 ˜ dgamma(12 ,1 )

40 mu1 ˜ dnorm(−3 ,0.1)

41 mu2 ˜ dnorm(−3 ,0.1)

42 ps i 0 ˜ dgamma(2 , 1 )

43 p s i ˜ dgamma(1 , 1 )

100

44 i t au2 ˜ dgamma(1 , 1 )

45

A.1.6 Modelo 6


1

2 model

3 4 f o r ( i in 1 :N)

5 6 ysdet [ i ] <− exp ( th [ i , 1 ] + th [ i , 2 ] ) * ( 1 − rhoep* rhoep )/w[ i ]




10 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]

11 w[ i ] ˜ dgamma( dstar , d s ta r )


13 Yis2a [ i ] <− exp(−th [ i , 1 ] )

14 Yis2b [ i ] <− exp(−th [ i , 2 ] )

15 sdY1 [ i ] <− 1/ s q r t ( Yis2a [ i ] )

16 sdY2 [ i ] <− 1/ s q r t ( Yis2b [ i ] )

17 18 muy [ 1 ] <− 0

19 muy [ 2 ] <− 0

20 thm [ 1 , 1 ] <− mu1

21 thm [ 1 , 2 ] <− mu2



24 f o r ( i in 2 :N)

25 26 thm [ i , 1 ] <− mu1 + phi1 *( th [ i −1 ,1] − mu1)

27 thm [ i , 2 ] <− mu2 + phi2 *( th [ i −1 ,2] − mu2)



30 31 phi1 ˜ dbeta (1 , 1 )

101

32 phi2 ˜ dbeta (1 , 1 )

33 i t aua2 ˜ dgamma(20 ,1 )

34 i taub2 ˜ dgamma(13 ,1 )

35 mu1 ˜ dnorm(−3 ,0.1)

36 mu2 ˜ dnorm(−3 ,0.1)

37 dstar ˜ dgamma(1 , 1 )

38 rhoep ˜ dun i f (−1 ,1)

39

A.2 Sao Paulo

A.2.1 Modelo 1


1

2 model

3 4 f o r ( i in 1 :N)

5 6 Y1var [ i ] <− exp ( th1 [ i ] )

7 Y1tau [ i ] <− 1/Y1var [ i ]

8 Y1 [ i ] ˜ dnorm (0 , Y1tau [ i ] )

9 SdY1 [ i ] <− s q r t ( Y1var [ i ] )

10 Y2var [ i ] <− exp ( th2 [ i ] )

11 Y2tau [ i ] <− 1/Y2var [ i ]

12 Y2 [ i ] ˜ dnorm (0 , Y2tau [ i ] )

13 SdY2 [ i ] <− s q r t ( Y2var [ i ] )

14 Y3var [ i ] <− exp ( th3 [ i ] )

15 Y3tau [ i ] <− 1/Y3var [ i ]

16 Y3 [ i ] ˜ dnorm (0 , Y3tau [ i ] )

17 SdY3 [ i ] <− s q r t ( Y3var [ i ] )

18 Y4var [ i ] <− exp ( th4 [ i ] )

19 Y4tau [ i ] <− 1/Y4var [ i ]

20 Y4 [ i ] ˜ dnorm (0 , Y4tau [ i ] )

21 SdY4 [ i ] <− s q r t ( Y4var [ i ] )

22 Y5var [ i ] <− exp ( th5 [ i ] )

102

23 Y5tau [ i ] <− 1/Y5var [ i ]

24 Y5 [ i ] ˜ dnorm (0 , Y5tau [ i ] )

25 SdY5 [ i ] <− s q r t ( Y5var [ i ] )

26 27 thm1 [ 1 ] <− mu1

28 thm2 [ 1 ] <− mu2

29 thm3 [ 1 ] <− mu3

30 thm4 [ 1 ] <− mu4

31 thm5 [ 1 ] <− mu5






37 f o r ( i in 2 :N)

38 39 thm1 [ i ] <− mu1 + phi1 *( th1 [ i −1] − mu1)


41 thm2 [ i ] <− mu2 + phi2 *( th2 [ i −1] − mu2)


43 thm3 [ i ] <− mu3 + phi3 *( th3 [ i −1] − mu3)


45 thm4 [ i ] <− mu4 + phi4 *( th4 [ i −1] − mu4)


47 thm5 [ i ] <− mu5 + phi5 *( th5 [ i −1] − mu5)


49 50 phi1 ˜ dnorm (0 ,10 )


52 mu1 ˜ dnorm ( 0 , 0 . 0 1 )

53 phi2 ˜ dnorm (0 ,10 )


55 mu2 ˜ dnorm ( 0 , 0 . 0 1 )

56 phi3 ˜ dnorm (0 ,10 )


58 mu3 ˜ dnorm ( 0 , 0 . 0 1 )

59 phi4 ˜ dnorm (0 ,10 )


61 mu4 ˜ dnorm ( 0 , 0 . 0 1 )

103

62 phi5 ˜ dnorm (0 ,10 )


64 mu5 ˜ dnorm ( 0 , 0 . 0 1 )

65

A.2.2 Modelo 2


1

2 model

3 4 f o r ( i in 1 :N)

5 6 Y[ i , 1 : 5 ] ˜ dmnorm(muy [ ] , Yis22 [ i , 1 : 5 , 1 : 5 ] )

7 Yis2 [ i , 1 , 1 ] <− exp ( th [ i , 1 ] )

8 Yis2 [ i , 2 , 2 ] <− exp ( th [ i , 2 ] )

9 Yis2 [ i , 3 , 3 ] <− exp ( th [ i , 3 ] )

10 Yis2 [ i , 4 , 4 ] <− exp ( th [ i , 4 ] )

11 Yis2 [ i , 5 , 5 ] <− exp ( th [ i , 5 ] )

12 Yis2 [ i , 1 , 2 ] <− rhoep12*exp (0 . 5* th [ i , 1 ] + 0 .5* th [ i , 2 ] )

13 Yis2 [ i , 2 , 1 ] <− Yis2 [ i , 1 , 2 ]


15 Yis2 [ i , 3 , 1 ] <− Yis2 [ i , 1 , 3 ]


17 Yis2 [ i , 4 , 1 ] <− Yis2 [ i , 1 , 4 ]


19 Yis2 [ i , 5 , 1 ] <− Yis2 [ i , 1 , 5 ]


21 Yis2 [ i , 3 , 2 ] <− Yis2 [ i , 2 , 3 ]


23 Yis2 [ i , 4 , 2 ] <− Yis2 [ i , 2 , 4 ]


25 Yis2 [ i , 5 , 2 ] <− Yis2 [ i , 2 , 5 ]


27 Yis2 [ i , 4 , 3 ] <− Yis2 [ i , 3 , 4 ]


29 Yis2 [ i , 5 , 3 ] <− Yis2 [ i , 3 , 5 ]

104


31 Yis2 [ i , 5 , 4 ] <− Yis2 [ i , 4 , 5 ]

32 sdY1 [ i ] <− 1/ s q r t ( Yis22 [ i , 1 , 1 ] )

33 sdY2 [ i ] <− 1/ s q r t ( Yis22 [ i , 2 , 2 ] )

34 sdY3 [ i ] <− 1/ s q r t ( Yis22 [ i , 3 , 3 ] )

35 sdY4 [ i ] <− 1/ s q r t ( Yis22 [ i , 4 , 4 ] )

36 sdY5 [ i ] <− 1/ s q r t ( Yis22 [ i , 5 , 5 ] )

37 sdY11 [ i ] <− s q r t ( Yis2 [ i , 1 , 1 ] )

38 sdY22 [ i ] <− s q r t ( Yis2 [ i , 2 , 2 ] )

39 sdY33 [ i ] <− s q r t ( Yis2 [ i , 3 , 3 ] )

40 sdY44 [ i ] <− s q r t ( Yis2 [ i , 4 , 4 ] )

41 sdY55 [ i ] <− s q r t ( Yis2 [ i , 5 , 5 ] )

42 Yis22 [ i , 1 : 5 , 1 : 5 ] <− i n v e r s e ( Yis2 [ i , 1 : 5 , 1 : 5 ] )

43 44 muy [ 1 ] <− 0

45 muy [ 2 ] <− 0

46 muy [ 3 ] <− 0

47 muy [ 4 ] <− 0

48 muy [ 5 ] <− 0

49 thm [ 1 , 1 ] <− mu1

50 thm [ 1 , 2 ] <− mu2

51 thm [ 1 , 3 ] <− mu3

52 thm [ 1 , 4 ] <− mu4

53 thm [ 1 , 5 ] <− mu5



56 th [ 1 , 3 ] ˜ dnorm(thm [ 1 , 3 ] , i t auc2 )

57 th [ 1 , 4 ] ˜ dnorm(thm [ 1 , 4 ] , i taud2 )

58 th [ 1 , 5 ] ˜ dnorm(thm [ 1 , 5 ] , i t aue2 )

59 f o r ( i in 2 :N)

60 61 thm [ i , 1 ] <− mu1 + phi1 *( th [ i −1 ,1] − mu1)

62 thm [ i , 2 ] <− mu2 + phi2 *( th [ i −1 ,2] − mu2)

63 thm [ i , 3 ] <− mu3 + phi3 *( th [ i −1 ,3] − mu3)

64 thm [ i , 4 ] <− mu4 + phi4 *( th [ i −1 ,4] − mu4)

65 thm [ i , 5 ] <− mu5 + phi5 *( th [ i −1 ,5] − mu5)



68 th [ i , 3 ] ˜ dnorm(thm [ i , 3 ] , i t auc2 )

105

69 th [ i , 4 ] ˜ dnorm(thm [ i , 4 ] , i taud2 )

70 th [ i , 5 ] ˜ dnorm(thm [ i , 5 ] , i t aue2 )

71 72 phi1 ˜ dnorm ( 0 . 4 1 0 2 , 2 5 . 1 5 )

73 phi2 ˜ dnorm ( 0 . 7 0 3 7 , 1 0 3 . 5 3 )

74 phi3 ˜ dnorm ( 0 . 4 5 8 2 , 6 3 . 0 9 )

75 phi4 ˜ dnorm ( 0 . 6 5 7 4 , 6 2 . 3 9 )

76 phi5 ˜ dnorm ( 0 . 3 1 4 8 , 3 5 . 8 6 )

77 i t aua2 ˜ dgamma( 1 1 . 4 2 , 2 . 3 0 )

78 i taub2 ˜ dgamma( 1 0 . 3 8 , 2 . 4 8 )

79 i t auc2 ˜ dgamma( 1 3 . 3 8 , 4 . 6 1 )

80 i taud2 ˜ dgamma( 9 . 7 9 , 2 . 2 4 )

81 i t aue2 ˜ dgamma( 1 0 . 1 8 , 2 . 3 9 )

82 mu1 ˜ dnorm (−2.173 ,167 .79)

83 mu2 ˜ dnorm (−2 .783 ,88 .83)

84 mu3 ˜ dnorm (−2.713 ,129 .13)

85 mu4 ˜ dnorm (−2.647 ,107 .01)

86 mu5 ˜ dnorm (−2.809 ,173 .36)

87 rhoep12 ˜ dnorm (0 ,10 )

88 rhoep13 ˜ dnorm (0 ,10 )

89 rhoep14 ˜ dnorm (0 ,10 )

90 rhoep15 ˜ dnorm (0 ,10 )

91 rhoep23 ˜ dnorm (0 ,10 )

92 rhoep24 ˜ dnorm (0 ,10 )

93 rhoep25 ˜ dnorm (0 ,10 )

94 rhoep34 ˜ dnorm (0 ,10 )

95 rhoep35 ˜ dnorm (0 ,10 )

96 rhoep45 ˜ dnorm (0 ,10 )

97

Documents

Aplica˘c~ao de Modelos de Volatilidade Estoc astica em ... · Social da Faculdade de Medicina de Ribeir~ao Preto da Universidade de Sao Paulo para a obtenc~ao do t tulo de Mestre