Universidade Federal do Rio de Janeiro Instituto de Matem ... · Universidade Federal do Rio de Janeiro Instituto de Matem atica Modelos de convolu˘c~ao para dados espa˘co-temporais

Universidade Federal do Rio de JaneiroInstituto de Matematica

Modelos de convolucao para dadosespaco-temporais

Geraldo Marcelo da Cunha

2009



Tese de Doutorado submetida ao programa de

Pos-graduacao em Estatıstica do Instituto de

Matematica da Universidade Federal do Rio de

Janeiro como parte dos requisitos necessarios a

obtencao do grau de Doutor em Estatıstica.

Orientador: Dani Gamerman

Rio de Janeiro, julho de 2009




Tese de Doutorado submetida ao programa de Pos-graduacao em Estatıstica do

Instituto de Matematica da Universidade Federal do Rio de Janeiro como parte dos

requisitos necessarios a obtencao do grau de Doutor em Estatıstica.

Aprovada por:

Presidente Prof. Dani Gamerman Profa. Esther Salazar

IM–UFRJ IM–UFRJ

Prof. Hedibert Freitas Lopes Prof. Ronaldo Dias

University of Chicago UNICAMP

Prof. Josemar Rodrigues

UFSCAR

Rio de Janeiro, Julho de 2009

Cunha, Geraldo Marcelo

Modelos de convolucao para dados espaco-temporais/

Geraldo Marcelo da Cunha. – Rio de Janeiro: UFRJ/IM,

2009.

xvii, 133 f. : il. ; 31cm.

Tese (Doutorado) – UFRJ/IM. Programa de Pos-

Graduacao em Estatıstica, 2009.


Referencias bibliograficas: p.101–106.

1. Estatıstica Matematica - Tese. I. Gamerman, Dani.

II. Universidade Federal do Rio de Janeiro. Instituto de

Matematica. III. Tıtulo.

Modelos de convolucao para dados espaco-temporais 1

Resumo

Nos recentes anos, uma serie de modelos para dados espacialmente correlacionados

surgiram no intuito de relaxar ou desconsiderar a suposicao de estacionaridade na

estrutura de covariancia. Alguns desses modelos utilizam a ideia de convolucao de

processos por funcoes nucleo. Partindo de modelos conhecidos na literatura (Fuentes

e Smith (2001); Fuentes (2002); Banerjee, Gelfand, Knight, e Sirmans (2004b)) um

novo modelo que considera estrutura de covariancia nao-estacionaria e apresentado. O

modelo e definido como uma mistura de processos estacionarios latentes ponderados por

componentes de misturas. Funcoes nucleo contınuas sao utilizadas para definir estas

componentes. Nossa abordagem consiste em a partir dos dados, estimar os parametros

de suavizacao das funcoes nucleo consideradas.

Seguindo adiante no trabalho, o modelo proposto para processos espaciais nao-

estacionarios e inserido no contexto de dados que variam no tempo e espaco. Dois

diferentes modelos dinamicos espaco-temporais sao propostos. O primeiro, considera

uma mesma estrutura espacial nao-estacionaria em todos os tempos. O segundo, ge-

neraliza o primeiro, permitindo tambem, que esta estrutura espacial evolua no tempo

atraves de seus parametros de suavizacao.

No ultimo capıtulo, dados de temperatura mınima mensal observados no estado do

Rio de Janeiro de 1961 a 2000 sao analisados considerando um modelo hierarquico que

incorpora nossa abordagem para a estrutura espacial nao-estacionaria.


Abstract

In recent years, different models for correlated spatial data emerged in order to relax

or disregard the stationarity assumption of covariance structure. Some of these mod-

els are built around kernel convolution of stationary processes. Starting from known

models in the literature (Fuentes e Smith (2001); Fuentes (2002); Banerjee, Gelfand,

Knight, e Sirmans (2004b)) a new model considering non-stationary covariance struc-

ture is presented. The model is defined as a weighted combination of latent stationary

processes by mixture components. Continuous kernel functions are used to define these

componentes. Our approach consists in estimation of the bandwidths in the kernels,

from the data.

Following on the work, the proposed model for spatial non-stationary processes is

inserted to account for data varying in time and space. Two different spatio-temporal

dynamic models are proposed. The first one considers the same non-stationary spatial

structure for all times. The second one generalizes the first one by allowing this spatial

structure to vary over time through their bandwidths parameters.

In the last chapter monthly minimum temperature data observed in Rio de Janeiro

state from 1961 to 2000 are analyzed considering a hierarchical model that incorporates

our non-stationary approach.


AGRADECIMENTOS

Inicialmente gostaria de agradecer ao meu orientador Dani Gamerman, pessoa pela

qual passei a admirar nos ultimos anos como pesquisador e ser humano. Este tra-

balho e fruto de muitos dos seus conselhos e intuicoes, ditos ou escritos muitas vezes

em pequenos rascunhos em pedacos de papel, ou em seu pequeno quadro branco, o

que sempre me surpreendia. Eu espero que esta relacao forneca frutos de amizade e

trabalhos futuros.

Gostaria tambem de agradecer pelas outras pessoas maravilhosas que conheci e que

participaram deste processo, todos professores do programa e meus amigos e colegas da

pos: Valmaria, Mario, Adelmo, Fidel, Esther e Fernando. Meus amigos da FIOCRUZ

que sempre me apoiaram e possibilitaram minha visita a Universidade Estadual da

Carolina do Norte.

Gostaria de agradecer a minha famılia e amigos distantes que sao o meu suporte

e dao razao ao meu viver. Meu pais, Mario e Ana, meus irmaos, Rosana e Marcio e

outros do coracao, Buda, Ivan, Alan, Alex primo, Didier, Alex Erikson, Edicleia, Tio

Romulo. A Lane que esteve comigo em boa parte deste perıodo. E finalmente, a Thati,

minha companheira, amiga e mulher que sempre me incentivou e incentiva em tudo o

que faco.

Durante este doutorado, outras coisas alegres e tristes ocorreram em minha vida,

sou grato a todas elas. Citando meu preferido poeta, Fernando Pessoa:

“...

Quem quere passar alem do Bojador

Tem que passar alem da dor.

Deus ao mar o perigo e o abysmo deu,

Mas nelle e que espelhou o ceu.”

Sumario

1 Processos espaciais via convolucao de processos 7

1.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.2 Modelos para processos espaciais nao-estacionarios . . . . 10

1.3 Processos espaciais via convolucao de processos . . . . . 12

1.3.1 Processos estacionarios via convolucao de processos 12

1.3.2 Processos nao-estacionarios via convolucao de pro-

cessos com nucleos variando no espaco . . . . . . 16

1.3.3 Processos nao-estacionarios via convolucao de pro-

cessos localmente estacionarios . . . . . . . . . . . 17

1.4 Justificativa deste trabalho . . . . . . . . . . . . . . . . . 20

1.5 Organizacao da tese . . . . . . . . . . . . . . . . . . . . . 20

2 Processos nao-estacionarios via convolucao de processos

estacionarios latentes 24

2.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.2 Definicao do modelo e suas propriedades . . . . . . . . . 25

2.3 Inferencia . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3


2.3.1 Introducao . . . . . . . . . . . . . . . . . . . . . . 28

2.3.2 Verossimilhanca . . . . . . . . . . . . . . . . . . . 29

2.3.3 Distribuicoes a Priori e a Posteriori . . . . . . . . 29

2.4 Aspectos computacionais . . . . . . . . . . . . . . . . . . 32

2.4.1 Interpolacao . . . . . . . . . . . . . . . . . . . . . 35

2.5 Simulacoes . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.5.1 Estudo 1 . . . . . . . . . . . . . . . . . . . . . . 36

2.5.2 Estudo 2 . . . . . . . . . . . . . . . . . . . . . . 44

3 Modelo dinamico com estrutura de covariancia espacial

nao-estacionaria 51

3.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.2 Modelo Proposto . . . . . . . . . . . . . . . . . . . . . . 53

3.3 Inferencia . . . . . . . . . . . . . . . . . . . . . . . . . . 55


3.3.2 Distribuicao a Priori . . . . . . . . . . . . . . . . 55

3.3.3 Distribuicao a Posteriori . . . . . . . . . . . . . . 56


3.5 Simulacao . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4 Modelo dinamico com estrutura de covariancia espacial

dinamica nao-estacionaria 66

4.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.2 Modelo Geral . . . . . . . . . . . . . . . . . . . . . . . . 67


4.3 Inferencia . . . . . . . . . . . . . . . . . . . . . . . . . . 72


4.3.2 Distribuicao a Priori . . . . . . . . . . . . . . . . 72

4.3.3 Distribuicao a Posteriori . . . . . . . . . . . . . . 73


4.5 Simulacao . . . . . . . . . . . . . . . . . . . . . . . . . . 77

5 Aplicacao 81

5.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . 81

5.2 Modelo para variacoes na temperatura . . . . . . . . . . 82

5.2.1 Modelando a estrutura nao-estacionaria . . . . . . 83

5.2.2 Inferencia dos parametros do modelo . . . . . . . 85

5.2.3 Aspectos computacionais . . . . . . . . . . . . . . 86

5.2.4 Dados faltantes e interpolacao . . . . . . . . . . . 86

5.3 Exemplo simulado . . . . . . . . . . . . . . . . . . . . . . 87

5.4 Resultados da analise . . . . . . . . . . . . . . . . . . . . 88

5.4.1 Descricao dos dados e modelo . . . . . . . . . . . 88

5.4.2 Comparacao dos modelos . . . . . . . . . . . . . . 90

5.4.3 Estimacao do termo espacial . . . . . . . . . . . . 93

5.4.4 Estimacao da estrutura media . . . . . . . . . . . 93

5.5 Conclusao . . . . . . . . . . . . . . . . . . . . . . . . . . 98

6 Consideracoes finais 100


Referencias Bibliograficas 103

Capıtulo 1

Processos espaciais via convolucao

de processos

1.1 Introducao

Nesta secao fazemos uma revisao breve e simples de modelos para dados espacialmente

correlacionados. O objetivo e apresentar algumas definicoes da teoria classica e chamar

atencao de como algumas propriedades desses modelos foram construıdas de modo a

facilitar a estimacao e interpretacao de seus parametros.

Considere uma regiao espacialmente contınua1 D ⊂ <2 na qual, somente para um

conjunto n de posicoes fixas s1, . . . , sn, sao conhecidas medidas de interesse Y =

(Y (s1), . . . , Y (sn))′. O objetivo final da analise estatıstica espacial e ser capaz de

fornecer a qualquer localizacao s0 ∈ D, onde nao se conhece o valor da medida de

interesse, a melhor estimativa Y (s0), a partir dos dados observados Y .

A funcao aleatoria Y (.) e chamada campo aleatorio. Uma abordagem usual na

modelagem de valores observados de um campo aleatorio e feita ao considerar Y de-

composto nas seguintes componentes:

1Pode-se considerar uma regiao definida em um espaco de dimensao L qualquer, onde L ≥ 1 e L e

um numero inteiro.

7


Y = (media)+(componente espacial)+(erro de medida) (1.1)

= µ+Z + ε,

onde µ representa uma media que pode, por exemplo, ter forma linear, isto e, µ = Xβ,

onde X e uma matriz (n × p) que acomoda p covariaveis associadas a media e β e

um vetor (p × 1) dos parametros associados a essas covariaveis. A componente Z

corresponde a um vetor (n× 1) de realizacao de um campo aleatorio espacial de media

0 e ε e um vetor (n×1) de erros de medida independentes e identicamente distribuıdos,

tal que ε e Z sao independentes. ε e comumente chamado de efeito pepita (Cressie,

1993).

Em muitas aplicacoes, o campo aleatorio espacial Z(.) e assumido Gaussiano, ou

diz-se simplesmente que Z(.) segue um processo Gaussiano.

De um modo geral, Z(.) segue um processo Gaussiano de media 0, se para qualquer

n ≥ 1, o conjunto de observacoes Z = (Z(s1), . . . , Z(sn))′ tem distribuicao normal

multivariada com media 0 e matriz de covariancia Ω. Vamos considerar a seguinte

notacao,

Z(.) ∼ PG(0,Ω).

No contexto da teoria espacial classica, e utilizada a simplificacao Ω = σ2Σ, onde

Σ e a matriz de correlacao dos dados e σ2 e a variancia, igual para todos os Z(s).

Assumindo que Y e condicionalmente independente dado os parametros e que o

processo Z(.) segue um processo Gaussiano podemos escrever:

Y | β,Z, σ2ε ∼ N

(Xβ +Z, σ2

ε I), (1.2)

onde σ2ε e a variancia do erro de medida ε.

Da expressao acima, podemos marginalizar em Z, obtendo

Y | β, σ2,Σ, σ2ε ∼ N

(Xβ, σ2Σ + σ2

ε I). (1.3)


Outra simplificacao util ocorre pela especificacao da matriz Σ definida a partir de

uma unica funcao de correlacao ρ entre medidas realizadas em dois pontos quaisquer

si e sj,

ρij = ρ(Z(si), Z(sj)).

Como esperado, ρij deve ser definido como uma funcao que decresca com o aumento

da distancia que separa as localizacoes si e sj. A geoestatıstica classica (Cressie, 1993)

apresenta uma classe parametrica de funcoes de correlacao ρ que dependem somente

do vetor diferenca entre os pontos si e sj.

Sendo mais especıfico, se ν = si− sj e o vetor da diferenca entre duas localizacoes

quaisquer si e sj e ρ e funcao somente dessa diferenca, ρ = ρ(ν), Z(.) e estacionario,

significando que a correlacao entre quaisquer dois pontos em D depende da orientacao

do vetor da diferenca entre os pontos e do modulo dessa diferenca. Caso ρ dependa

somente da distancia entre as localizacoes, ρ = ρ(|ν|), Z(.) e isotropico, e portanto, ρ

passa a nao depender mais da orientacao do vetor diferenca. Sendo assim, segue que

todo processo isotropico e tambem estacionario.

Como exemplo, temos abaixo a funcao de correlacao estacionaria Matern:

ρ =1

2ξ−1Γ(ξ)(2ξ1/2 | ν | φ)νκξ(2ξ

1/2φ | ν |),

onde κν e a funcao modificada de Bessel de terceira ordem. Neste caso, a matriz

de covariancias e totalmente especificada se conhecemos os parametros σ2, φ e ξ. O

parametro φ e o inverso da amplitude r (range) que mede como a correlacao decai

com a distancia. O parametro ξ e responsavel pela suavidade do processo sendo este

mais suave na medida em que ξ aumenta. O parametro σ2 e a variancia do processo.

Outras importantes funcoes de correlacao sao derivadas a partir da funcao de correlacao

Matern. Por exemplo, fazendo ξ → ∞ obtemos a funcao de correlacao gaussiana.

Fixando ξ = 12, obtemos a funcao de correlacao exponencial:

ρ = exp −φ | ν | .


Apesar de atraentes do ponto de vista da interpretacao e facil implementacao,

os modelos que consideram processos espaciais estacionarios/isotropicos sao limitados

dadas as suas especificacoes. Essas propriedades limitam uma possıvel complexidade

real na estrutura espacial que dados reais possam vir a ter. Recorre-se portanto a mode-

los espaciais mais elaborados, onde por exemplo, seja permitido a funcao de covariancia

variar sobre a regiao de estudo.

1.2 Modelos para processos espaciais nao-estacionarios

Devido ao avanco e maior utilizacao de procedimentos computacionalmente intensivos

nos recentes anos, uma serie de procedimentos surgiram no intuito de relaxar ou des-

considerar a suposicao de estacionaridade para dados espaciais.

Uma serie de procedimentos que obtiveram notorio reconhecimento na literatura

parte do artigo inicial de Sampson e Guttorp (1992) que sugerem um modelo de de-

formacao espacial para dados espaco-temporais onde a estacionaridade e assumida no

processo temporal enquanto e permitido ao processo espacial ser nao-estacionario. Eles

utilizam tecnicas de escalonamento multidimensional para transformar o espaco onde

se encontram as observacoes, em um espaco latente estacionario no qual procedimentos

de geoestatıstica classica podem ser aplicados. Seguindo essa linha, Schmidt e O’Hagan

(2003) propoem uma abordagem Bayesiana para este modelo.

Hass (1995) propoe uma abordagem baseada em janelas moveis para estruturas de

covariancia nao-estacionarias. Entretanto, uma desvantagem deste modelo e que ele

pressupoe isotropia nas janelas que incidem sobre as localizacoes.

Uma outra abordagem proposta por Nychka e Saltzman (1998) utiliza funcoes or-

togonais empıricas para modelar dados com estruturas de covariancia mais complexas.

Entretanto, o inconveniente deste metodo e a necessidade de que as localizacoes, onde

os dados foram coletados, sejam um subconjunto da grade onde e feita a interpolacao

dos dados.

Outras duas diferentes abordagens surgem da convolucao de processos por funcoes

nucleo (kernel). A primeira abordagem (Higdon, Swall, e Kern (1999); Swall (1999))


permite que nucleos normais bivariados possam variar no espaco e considera os proces-

sos como sendo um mesmo ruıdo branco ou um mesmo processo espacial estacionario

definido por um mesmo conjunto de parametros. Paciorek e Schervish (2004) e Pa-

ciorek e Schervish (2006) mostraram ser possıvel, a partir de outros nucleos, estabelecer

diferentes estruturas de covariancias nao-estacionarias. A segunda abordagem (Fuentes

e Smith (2001); Fuentes (2002)) utiliza um mesmo nucleo para diferentes centros,

mas processos espaciais estacionarios definidos por diferentes conjuntos de parametros.

Banerjee et al. (2004b) estende este modelo ao substituir os nucleos por funcoes ponde-

radas dos nucleos. Todos os trabalhos que consideram esta segunda abordagem fixam

os parametros de suavizacao em valores conhecidos, como sera discutido na proxima

secao.

Gelfand, Kottas, e MacEachern (2005b) definem processos nao estacionarios nao-

parametricos a partir da convolucao de processos de Dirichlet. Kottas, Duan, e Gelfand

(2008) utilizam uma abordagem equivalente a partir de um processo de Dirichlet cen-

trado em torno de uma normal multivariada. Ainda neste artigo, o modelo e estendido

para analise de dados espaco-temporais atraves de uma formulacao dinamica (West

e Harrison, 1997) para os efeitos aleatorios. Fuentes e Reich (2009) permitem aos

parametros de suavizacao que compoe a convolucao de processos de Dirichlet vari-

arem no espaco, para caracterizar a falta de estacionaridade na dependencia espacial

e dependencia cruzada, quando sao considerados processos espaciais multivariados.

Tambem, Fuentes, Henry, e Reich (2009) permitem ao parametro de suavizacao do

nucleo ser funcao espacial para explicar falta de estacionaridade em distribuicoes de

extremos, utilizando uma mistura nao-parametrica para dados de extremos de tempe-

ratura.

Para processos nao-estacionarios multivariados Schmidt e Gelfand (2003) e Gelfand,

Schmidt, Banerjee, e Sirmans (2004) utilizam modelos de coregionalizacao. Tambem

Calder (2003) e Calder (2007) propoem um modelo nao-estacionario com nucleos va-

riando espacialmente para dados espaco-temporais multivariados considerando modelos

dinamicos.

Ha ainda uma serie de modelos propostos na literatura para processos espaciais nao-


estacionarios. Uma revisao dos principais modelos pode ser encontrada em Banerjee,

Carlin, e Gelfand (2004a) e Smith (2001). Em portugues, a revisao de alguns desses

modelos pode ser encontrada em Schmidt e Sanso (2006).

A proxima secao e dedicada a uma melhor compreensao de modelos espaciais nao-

estacionarios envolvendo convolucao de processos.

1.3 Processos espaciais via convolucao de processos

1.3.1 Processos estacionarios via convolucao de processos

Considere uma grade de M pontos u1, ...,uM contidos em uma regiao D, como na

Figura 1.1 e suponha localizacoes (s1, . . . , sn) onde sao observadas as realizacoes de

um processo de interesse Z = Z(s1), . . . , Z(sn). Vejamos como isto pode ser obtido

a partir da convolucao de processos, de modo que Z represente a realizacao de um

processo espacial estacionario.

Inicialmente a cada um, m = 1, . . . ,M sao atribuıdas variaveis aleatorias ω(um)

independentes e identicamente distribuıdas,

ω(um) ∼ N(0, σ2

).

Deste modo, ω(um) representa um processo ruıdo branco.

Seja k(si − um;hm) o valor de uma funcao nucleo de centro um avaliada em si,

i = 1, . . . , n. O parametro hm e chamado de parametro de suavizacao. Quanto maiores

os valores de hm, m = 1, . . . ,M , mais suave e o processo final observado. De agora

em diante, a seguinte notacao sera considerada. Para um nucleo contınuo qualquer

k(.), onde u varia continuamente no espaco, k(si − u) = ku(si) representa o nucleo

centrado em u e avaliado em si. No caso em que u varia discretamente no espaco,

k(si − um) = km(si) representa o nucleo centrado em um e avaliado em si. Alem

disto, subentende-se um parametro de suavizacao para definir cada uma destas funcoes

nucleo.

Considerando ainda o caso discreto, se todos os parametros de suavizacao hm sao


0 2 4 6 8 10

02

46

810

u1 u2 u3 ……

…… u99 u100

s1

s2

s3

s4

s5

s6

s7

s8

s9

s10

s11s12

s13

s14

s15s1

s2

s3

s4

s5

s6

s7

s8

s9

s10

s11s12

s13

s14

s15

Figura 1.1: Grade formada de um pontos, m = 1, . . . , 100 em D = [0, 10] × [0, 10]. A cada um

e atribuıdo um processo ruıdo branco ω(um). A figura apresenta tambem o que poderiam ser as

possıveis localizacoes si, i = 1, . . . , 15 das medidas de interesse Z(si).

iguais para m = 1, . . . ,M , a realizacao de um processo estacionario Z(si) pode ser

obtida por convolver os processos ω(um) por km(si) atraves da soma finita,

Z(si) =M∑m=1

km(si)ω(um). (1.4)

Segue entao que Z e a realizacao de um processo Gaussiano no qual,

E[Z(si)] = E

[M∑m=1

km(si)ω(um)

]

=M∑i=m

km(si)E[ω(um)]

= 0,


V ar[Z(si)] = E

( M∑m=1

km(si)ω(um)

)2

= σ2

M∑m=1

km(si)2,

e:

Cov[Z(si), Z(sj)] = E[Z(si)Z(sj)]− E[Z(si)]E[Z(sj)]

= σ2

M∑m=1

km(si)km(sj)− 0

= σ2

M∑m=1

km(si)km(sj).

Na Figura 1.2 o processo estacionario foi gerado sobre a grade de pontos repre-

sentada pela Figura 1.1 considerando uma funcao nucleo normal padrao, km(s) =

12π

exp−1

2s′s

, ou seja, hm = 1 para m = 1, . . . ,M e um processo ruıdo branco de

variancia σ2 = 0.01.

A versao contınua de convolucao de processos, segue diretamente do limite da grade

de pontos formada por um cada vez mais densa,

Z(s) =

∫ku(s)ω(u)du. (1.5)

A variancia e covariancia de Z(si) sao dadas por,

V ar[Z(si)] = σ2

∫(ku(si))

2du, (1.6)

Cov[Z(si), Z(sj)] = σ2

∫ku(si)ku(sj)du. (1.7)

Uma demonstracao interessante da aproximacao discreta para a contınua de con-

volucao de processos pode ser encontrada em Smith (2001).

Na equacao (1.7), se δ = sj − u,


Figura 1.2: A realizacao de um processo estacionario e obtida a partir da convolucao do nucleo

normal padrao km(si) e um processo ruıdo branco ω(um) de variancia 0.01.


∫k(si − u)k(sj − u)du

= σ2

∫k(si − sj + δ)k(δ)dδ

que depende somente do vetor da diferenca entre as localizacoes ν = si−sj e portanto,

Z(.) e estacionario.

A equacao (1.4) pode ser utilizada para gerar processos estacionarios cujas funcoes

de correlacao podem ser obtidas na sua forma explıcita, utilizando a relacao entre

a transformada de Fourier de Cov[si − sj ] e a transformada de Fourier do nucleo

k(si − sj). Uma discussao sobre este tema pode ser encontrada em Kern (2000).

Cabe ressaltar que os ruıdos brancos ω(um) podem ser substituıdos por realizacoes

de um mesmo processo estacionario definido por um conjunto de parametros fixados e

ainda assim, o processo Z(.) sera estacionario.

Na pratica, a versao discreta do modelo (1.4) e utilizada pois a integracao das es-

truturas de covariancias estabelecidas por (1.5) e de difıcil manipulacao por metodos


numericos. Uma das vantagens atribuıdas ao modelo (1.4) e a reducao da dimensionali-

dade no problema, quando M e pequeno em relacao ao numero de localizacoes (Higdon

(1998); Kern (2000)).

Calder (2003) estende o modelo discreto de convolucao de processos (1.4) para

modelos dinamicos (West e Harrison, 1997). Isso e feito ao definir componentes de

ruıdo ωt(um), m = 1, . . . ,M , t = 1, . . . , T , evoluindo no tempo.

1.3.2 Processos nao-estacionarios via convolucao de processos

com nucleos variando no espaco

No intuito de obter processos nao-estacionarios, duas diferentes abordagens sao moti-

vadas a partir de (1.4) e (1.5).

Uma primeira abordagem para processos nao-estacionarios proposta por Higdon

et al. (1999) e Swall (1999) e obtida ao se manterem os processos ruıdo branco ω(u),

mas serem definidos nucleos ksi(.) que variam de acordo com a localizacao si. O modelo

e entao representado por,

Z(si) =

∫ksi(u)ω(u)du, (1.8)

onde ksi(.) e um nucleo de centro si.

De modo que a funcao de covariancia do processo e dada por,


∫ksi(u)ksj(u)du. (1.9)

Uma vez que Cov[Z(si), Z(sj)] depende tambem das localizacoes de si e sj, o

processo Z(.) e nao-estacionario.

Higdon et al. (1999) e Swall (1999) escolhem propositalmente a forma do nucleo,

ks(.), como a de uma normal bivariada com matriz de covariancias Λ(s). Como as

curvas de nıvel deste nucleo tem forma de elipse, eles utilizam elementos da equacao de

uma elipse (os focos, a area e um fator de expansao) como um modo de parametrizar

os elementos da matriz Λ(s). Alem disso, eles estabelecem que os focos dessa elipse


sigam um processo Gaussiano. A variacao dos focos permite que o modelo se ajuste a

estruturas complexas dos dados podendo eles exibir um comportamento suave ou nao.

Novamente na pratica, a integracao em (1.8) e aproximada pela soma finita,

Z(si) =M∑m=1

ksi(um)ω(um), (1.10)

sendo portanto a funcao de covariancia


M∑m=1

ksi(um)ksj(um). (1.11)

1.3.3 Processos nao-estacionarios via convolucao de processos

localmente estacionarios

Nosso trabalho tem como ponto de partida uma outra abordagem baseada na con-

volucao de processos, proposta inicialmente por Fuentes e Smith (2001) e Fuentes

(2002). Esta abordagem difere da anterior por manter fixo o nucleo ku(s) e substituir

o unico processo ruıdo branco ω(u) por processos espaciais localmente estacionarios

Wη(u)(s) de media 0 e conjunto de parametros da funcao de covariancia η(u). Por

exemplo, se Wη(u)(.) representa um processo estacionario com funcao de correlacao

exponencial de parametros η(u) = (σ2u, φu), de modo que,

Wη(u)(.) ∼ PG(0, σ2

u exp −φu | si − sj |).

O modelo nao-estacionario e entao representado por,

Z(s) =

∫ku(s)Wη(u)(s)du. (1.12)

Fuentes e Smith (2001) resolvem a integral em (1.12) cobrindo a regiao de estudo

D por uma grade fina de pontos um, m = 1, . . . ,M . Definindo para cada um um

processo Wη(um)(.), os autores assumem que funcoes (por exemplo, logarıtmica) de cada

um dos parametros em η(um) variam no espaco segundo processos Gaussianos. Eles

impoem que cada nucleo ku(.) decresca rapidamente, a partir de seus parametros de


suavizacao de modo que para pontos si e sj proximos de um, Z(si) e Z(sj) sejam como

realizacoes de um mesmo processo estacionario de parametros η(um). Deste modo, o

processo observado global e nao-estacionario, mas retem localmente caracterısticas de

um processo estacionario.

Apesar das interessantes consideracoes teoricas e propriedades deste modelo descrito

em Fuentes e Smith (2001), parecem haver problemas praticos na convergencia dos

parametros por considerar a aproximacao da integral por uma grade fina de pontos

um (Barber (1999); Banerjee et al. (2004a)). Deste modo, a integral em (1.12) e de

fato aproximada por uma grade de pontos relativamente muito menor ao numero de

localizacoes (Fuentes (2002); Fuentes, Chen, Davis, e Lackmann (2005); Banerjee et al.

(2004b)).

Define-se portanto processos localmente estacionarios Wη(um)(si) = Wm(si), m =

1, . . . ,M , satisfazendo Cov[Wm(si),Wm′(sj)] = 0 para m 6= m′ de modo que assume-se

a versao discreta de (1.12),

Z(si) =M∑m=1

km(si)Wm(si). (1.13)

E facil ver que E[Z(si)] = 0, desde que os processos Wm(si) tem media 0 para

m = 1, . . . ,M . Pode-se mostrar tambem que,

Cov[Z(si), Z(sj)] =M∑m=1

km(si)km(sj)Cov[Wm(si),Wm(sj)], (1.14)

e que portanto o processo geral e nao-estacionario.

Nesta abordagem, os parametros de suavizacao hm dos nucleos km(.), m = 1, . . . ,M,

sao escolhidas e fixadas nos menores valores possıveis de modo que a estacionaridade

e aproximadamente obtida dentro de cada sub-regiao. Para isto, as amplitudes dos

parametros de suavizacao hm sao escolhidas em seus menores valores possıveis, satis-

fazendo km(si) > 0 para todo si em D . A proposta apresentada por Banerjee et al.

(2004b) estende o modelo ao substituir os nucleos km(si) por pesos ponderados dos

nucleos,


γm(si) =km(si)∑M

m′=1 km′(si),ou ainda, γm(si) = km(si)√∑M

m′=1 k2m′ (si)

; (1.15)

e o modelo (1.13) e entao reescrito como,

Z(si) =M∑m=1

γm(si)Wm(si). (1.16)

Estas duas abordagens apresentadas nos modelo (1.13) e (1.16) que fixam os parame-

tros de suavizacao hm, fazem uma particao de D similar a utilizada na Tesselagem de

Voronoi. Em outras palavras, o processo nao-estacionario e basicamente obtido pela

composicao de processos estacionarios ao inves de uma mistura destes processos. Como

um exemplo, considere uma regiao D = [0, 100] × [0, 100], com M = 4 centros como

dispostos na Figura 1.3. Considere tambem 40 localizacoes (s1, . . . , s40) e pesos pon-

derados (1.15, primeira equacao) definidos a partir de nucleos normais padrao, k(si) =

12π

exp−1

2(si)

′(si)

, ou seja, nucleos com parametros de suavizacao hm = 1. Tomando

como exemplo, a localizacao s37, temos2, k1(s37) = O(10−400), k2(s37) = O(10−311),

k3(s37) = O(10−280), k4(s37) = O(10−139), de modo que γ1(s37).= 0, γ2(s37)

.= 0,

γ3(s37).= 0, γ4(s37)

.= 1, e portanto, Z(s37) = W4(s37), independente dos valores ob-

servados em cada um dos processos W1(s37), W2(s37) e W3(s37). O mesmo ocorre para

cada uma das outras localizacoes que tem somente um peso relevante associado, de

modo que sao estabelecidas 4 sub-regioes descritas por 4 processos estacionarios. As-

sim, o processo resultante e globalmente nao-estacionario mas retendo uma estrutura

local estacionaria.

Esta abordagem leva a inversao de matrizes esparsas da dimensao dos dados. A

proposta adotada por Kim, Mallick, e Holmes (2006) define processos estacionarios

independentes em sub-regioes obtidas pela Tesselagem de Voronoi. Neste caso, ma-

trizes de dimensao do numero de localizacoes em cada uma das sub-regioes e que sao

invertidas, tornando esta abordagem computacionalmente mais atrativa.

2Se limn→∞f(n)g(n) = L, onde 0 ≤ L <∞, entao f(n) = O(g(n)).


1.4 Justificativa deste trabalho

Este trabalho e focado na extensao do modelo (1.16) por tratar cada um dos parametros

de suavizacao hm dos nucleos km(.) que definem γm(.), m = 1, . . . ,M , como parametros

do modelo. Essa abordagem e original uma vez que nao estamos propondo um modelo

de misturas de processos localmente estacionarios (1.13) e (1.16). Nosso modelo assume

processos estacionarios latentes que se misturam entre si, ponderados por componentes

de misturas a serem definidas no proximo capıtulo. Nosso modelo tambem se difere

dos modelos (1.13) e (1.16) por ser nao-estacionario nao so atraves dos processos esta-

cionarios latentes mas tambem atraves das componentes de mistura, uma vez que e

permitido a elas variarem no espaco.

Os procedimentos de inferencia e amostragem adotados para os parametros pos-

sui a vantagem de selecionar automaticamente o numero de componentes de misturas

envolvidas no problema. Este procedimento evita o uso de algoritmos de saltos re-

versıveis (Green, 1985) tornando o procedimento computacional relativamente simples.

Exemplos simulados do proximo capıtulo sugerem que mesmo considerando diferentes

criterios de mistura dos processos estacionarios envolvidos, os procedimentos de in-

ferencia e amostragem adotados sao eficazes na recuperacao dos parametros utilizados

na simulacao dos dados.

Este trabalho tambem se justifica por extender esta nova proposta para processos

espaciais nao-estacionarios em modelos espaco-temporais atraves de uma abordagem

de modelos dinamicos (West e Harrison, 1997) na qual duas propostas sao apresen-

tadas. A primeira extensao considera uma estrutura dinamica na media na qual a

estrutura espacial nao-estacionaria dos dados permanece constante ao longo tempo

enquanto, a segunda extensao, considera uma estrutura dinamica na media e uma

estrutura dinamica para a matriz de covariancia espacial nao-estacionaria dos dados.

1.5 Organizacao da tese

Esta tese e organizada da forma a seguir.


No Capıtulo 2 o modelo proposto para processos espaciais nao-estacionarios e

definido e suas propriedades sao estabelecidas. Partindo de uma abordagem Bayesiana,

o procedimento de inferencia e esquema de amostragem sao apresentados para os

parametros que compoem a funcao de covariancia do modelo. Dois exemplos simu-

lados sao apresentados. Estes estudos ilustram como o procedimento de inferencia e

de amostragem adotados sao capazes de identificar corretamente o numero de com-

ponentes do modelo, descartando aquelas nao existentes. Nestes estudos, o modelo

proposto e tambem comparado ao modelo de mistura de processos localmente esta-

cionarios no que diz respeito a capacidade preditiva destes modelos .

No Capıtulo 3, a funcao de covariancia do modelo proposto no Capıtulo 2 e inserida

no contexto de dados espaco-temporais. Para isso, e utilizada uma abordagem de mode-

los dinamicos para a media do processo (West e Harrison, 1997) na qual e sugerida a

matriz de covariancias das observacoes ser nao-estacionaria. O modelo proposto neste

capıtulo considera que esta estrutura de covariancia nao evolua no tempo. Todo o

procedimento de inferencia e esquemas de amostragem sao apresentados para todos os

parametros do modelo. Ao final, um exemplo simulado de regressao dinamica e uti-

lizado para exemplificacao. Para os parametros de variancia dos coeficientes dinamicos

da regressao duas prioris sao sugeridas, a Gamma Invertida e a Semi-Cauchy. Alem

disso, o efeito delas na estimativa a posteriori dos parametros e analisada.

O Capıtulo 4 estende o modelo estabelecido no Capıtulo 3 ao permitir que estrutura

de covariancia nao-estacionaria evolua no tempo. Isso e feito por definir uma estrutura

dinamica para os parametros de suavizacao. Por serem positivos, a distribuicao log-

normal e utilizada na construcao da estrutura dinamica desses parametros. Essas

distribuicoes sao parametrizadas de modo que a cada evolucao no tempo o parametro

de suavizacao tenha media igual ao valor do parametro de suavizacao no tempo anterior.

No final do capıtulo um exemplo simulado e utilizado para exemplificar este modelo

dinamico de estrutura de covariancia dinamica.

No Capıtulo 5 e apresentado um estudo das alteracoes climaticas sofridas no estado

do Rio de Janeiro a partir da analise dos dados de temperatura mınima mensal de

1961 a 2000. O modelo considerado e baseado em uma especificacao hierarquica onde


os dados sao modelados atraves de uma estrutura de regressao para a media e um

termo de ruıdo espacialmente estruturado com uma representacao nao-estacionaria.

Ao final, este modelo e comparado a outros que consideram para o termo espacial, um

modelo de mistura de processos localmente estacionarios e um unico processo espacial

estacionario.

O Capıtulo 6 traca os plano de trabalhos futuros e possıveis extensoes para os

modelos.


+

+

+

+

+

+

+

+

+

+

+

+

+

++

+

+

+

+

+

+

+

+

+

+

+

+

+

+

++

++

+

++

+

+

+

+

0 20 40 60 80 100

020

4060

8010

0

longitude

latit

ude 37

u1 u2

u3 u4

0 20 40 60 80 100

020

4060

8010

0

longitude

latit

ude

u1 u2

u3 u4

Figura 1.3: Localizacoes de 40 pontos (s1, . . . , s40) de medidas de interesse Z(si) e centros dos

processos localmente estacionarios (u1,u2,u3,u4) (acima) e como esses pontos sao alocados a cada

uma das sub-regioes (abaixo) a partir dos seus respectivos pesos γm(si), m = 1, . . . , 4. Por exemplo,

o ponto s37 e alocado a sub-regiao 4, pois s37 esta mais proximo de u4 que dos outros centros.

Capıtulo 2

Processos nao-estacionarios via

convolucao de processos

estacionarios latentes

2.1 Introducao

Neste capıtulo apresentamos uma nova representacao nao-estacionaria para processos

espaciais. O modelo e formalmente definido e suas propriedades sao estabelecidas.

Partindo de uma abordagem Bayesiana, o procedimento de inferencia e esquema de

amostragem sao apresentados para o conjunto de parametros que compoem o modelo.

Um estudo simulado e utilizado para mostrar como as estrategias de inferencia

e de amostragem adotadas selecionam o numero correto de misturas envolvidas na

construcao do processo espacial nao-estacionario. Neste estudo, sao comparadas as

capacidades preditivas do modelo proposto e modelo de mistura de processos localmente

estacionarios. Outro estudo e utilizado para mostrar que o modelo proposto e capaz

de identificar um processo nao-estacionario definido a partir da mistura de processos

localmente estacionarios.

24


2.2 Definicao do modelo e suas propriedades

O modelo proposto e introduzido aqui. Considere M sub-regioes em D de centros

um, m = 1, . . . ,M e para cada uma dessas sub-regioes defina processos estacionarios

latentes Wm(.) seguindo um processo Gaussiano de media 0 e funcao de covariancia

definida por parametros η(um). Os processos Wm(.) sao independentes, isto e, para

todo i, j = 1, · · · , n, Wm(si) e Wm′(sj) sao independentes para todo m 6= m′ e si,

sj ∈ D. Seja tambem,

γ(si − um) = γm(si)

a componente de mistura associada ao processo Wm(.).

A restricao de idenficabilidadeM∑m=1

γm(si) = 1, e considerada da mesma forma como

em Banerjee et al. (2004b).

Se km(si −um) = km(si) representa uma funcao nucleo, esta restricao e facilmente

satisfeita se as componentes de mistura γm(si) sao definidas como pesos relativos,

γm(si) =km(si)∑M

m′=1 km′(si).

Note que em nossa abordagem, os nucleos km(.) dependem de seus respectivos

parametros de suavizacao desconhecidos hm. Quanto maior o valor de hm maior a

influencia da componente m sobre o processo observado na regiao.

O processo nao-estacionario e obtido pela mistura (convolucao) discreta,

Z(s) =M∑m=1

γm(s)Wm(s), (2.1)

para cada s ∈ D.

A equacao (2.1) fornece uma media ponderada dos processos e pode ser vista como

uma generalizacao espacial da estimacao dos parametros de suavizacao de Nadaraya

(1964) e Watson (1964) em modelos de regressao.


Alem disso, e esperado que com uma escolha adequada de componentes, esta re-

presentacao finita forneca uma aproximacao representativa da mistura contınua subja-

cente,

∫γu(s)Wu(s)du, (2.2)

como sugerido em Fuentes e Smith (2001). Note que os centros e suas correspondentes

regioes de vizinhanca sao meramente uma aproximacao para a forma de convolucao

infinita. Deste modo, eles nao precisam necessariamente estar relacionados a alguma

caracterıstica local especıfica da regiao. De fato, eles sao um artefato para a apro-

ximacao finita e as componentes do modelo podem nao ter uma interpretacao espacial

associadas as suas localizacoes especıficas.

Dada estas especificacoes do modelo segue que a esperanca, variancia e covariancia

do processo obtido sao,

E [Z(si)] = E

[M∑m=1

γm(si)Wm(si)

]= 0,

V ar[Z(si)] = V ar

[M∑m=1

γm(si)Wm(si)

]

= E

[M∑m=1

γm(si)Wm(si)

]2

=M∑m=1

[γm(si)]2 V ar[Wm(si)] e


Cov[Z(si), Z(sj)] = Cov

[M∑m=1

γm(si)Wm(si),M∑

m′=1

γm′(sj)Wm′(sj)

]

= E

[(M∑m=1

γm(si)Wm(si)

)(M∑

m′=1

γm′(sj)Wm′(sj)

)]

=M∑m=1

M∑m′=1

γm(si)γm′(sj)E [Wm(si)Wm′(sj)]

=M∑m=1

γm(si)γm(sj)Cov [Wm(si),Wm(sj)] .

A nao-estacionaridade e obtida porque os pesos γm(si) podem variar no espaco

atraves de seus parametros de suavizacao hm e os processos Wm(.) podem tambem

variar no espaco atraves de seus parametros η(um) que os definem. Banerjee et al.

(2004a) observam que esta abordagem apresenta algumas similaridades qualitativas

com o modelo proposto por Higdon et al. (1999).

O modelo (1.16) proposto por Banerjee et al. (2004b) pode ser visto como um caso

especial do modelo (2.1) quando os valores de hm, m = 1, . . . ,M , sao pequenos o

suficiente para evitar uma interacao entre as componentes do modelo, definindo assim

processos localmente estacionarios. Neste caso, se uma localizacao si e mais proxima

do centro um que dos demais centros entao γm(si) ≈ 1 e γm′(si) ≈ 0 para todo

m′ 6= m. Entao Z(si) ≈ Wm(si). O mesmo raciocınio pode ser aplicado a todas

as localizacoes em D, implicando uma representacao de processo aproximadamente

estacionario para todas as localizacoes, mas variando no espaco uma vez que diferentes

localizacoes podem estar mais proximas de diferentes centros. Note tambem que no

caso em que M = 1 e h1 e suficientemente grande, γ1(si) = 1 para todo si ∈ D, e

portanto Z(si) = W1(si) define um unico processo estacionario em toda regiao.

Uma possıvel interpretacao para o modelo (2.1) e a seguinte: para todo ponto

s ∈ D, o processo observado Z(s) e uma mistura de processos estacionarios latentes

Wm(s), ponderados por componentes de mistura γm(s). Naturalmente, os parametros

η(um) das funcoes de covariancia Wm(s) estimados nao podem ser comparados com os

parametros estimados da funcao de covariancia de outros modelos, por exemplo, com


o modelo de mistura de processos localmente estacionarios (1.16).

E bastante conhecido o fato de que os parametros de suavizacao desempenham um

papel mais importante do que a forma do nucleo em modelos de mistura e suavizadores

via nucleo (ver por exemplo, Epanechnikov (1969)). O mesmo efeito e valido aqui.

Algumas suposicoes a respeito do nucleo e funcao de covariancia sao agora assumidas

para serem consideradas neste capıtulo e capıtulos subsequentes.

Vamos considerar o seguinte nucleo Gaussiano,

km(si − um) =1

2π|Λm|−

12 exp

−1

2(si − um)′Λ−1

m (si − um)

, (2.3)

onde Λm =

h2m 0

0 h2m

e que os processos estacionarios de mistura sejam isotropicos

e Gaussianos com variancia σ2m e funcao de correlacao espacial de amplitude 1/φm,

para m = 1, . . . ,M . Ou seja,

Wm(.) ∼ PG(0, σ2

m exp(−φm|si − sj|)), (2.4)

de modo que a funcao de covariancia do modelo proposto (2.1) sera dada por,

Cov[Z(si), Z(sj)] =∑M

m=1 γm(si)γm(sj)σ2m exp(−φm|si − sj|).

Note que outras possibilidades mais elaboradas na definicao do nucleo sao possıveis,

por exemplo, a utilizacao de nucleos Gaussianos anisotropicos correlacionados como

em Higdon et al. (1999) ou nucleos mais gerais como os estabelecidos em Paciorek e

Schervish (2006). Note tambem que outras diferentes funcoes de correlacao para os

processos estacionarios latentes poderiam ser consideradas.

2.3 Inferencia

2.3.1 Introducao

Sejam as observacoes ao longo do tempo denotadas por Z = (Z ′1, . . . ,Z′T )′, onde

Zt = (Zt(s1), . . . , Zt(sn))′. Dadas as especificacoes em (2.3) e (2.4), σ = (σ1, . . . , σM)


e φ = (φ1, . . . , φM) representam os parametros das funcoes de covariancia dos processos

Wm(.) e h = (h1, . . . , hM) os parametros de suavizacao que compoem as componentes

de mistura γm(.) O conjunto de parametros a serem estimados e Θ = (σ,φ,h). Vamos

supor inicialmente que as observacoes Zt, t = 1, . . . T sao independentes para t 6= t′ e

que a estrutura espacial nao-estacionaria seja a mesma para todos os tempos,

Z|Θ ∼ N(0,Ω⊗ IT ), (2.5)

onde Ω =∑M

m=1σ2mAmR(φm)Am, Am = diag(γm(s1), . . . , γm(sn)) , [R(φm)]i,j =

exp(−φm|si − sj|) e In e a matriz identidade n-dimensional. Esta especificacao e

tal que (2.5) pode ser facilmente incorporado a um modelo mais geral, como um ruıdo

espacialmente estruturado nao-estacionario.

A partir de cada uma das amostras observadas Zt = (Zt(s1), . . . , Zt(sn))′ o objetivo

e interpolar Zt(.) nao observado em localizacoes (sn+1, . . . , sn+k). As respostas nao-

observadas nessas localizacoes, Znot = (Zt(sn+1), . . . , Zt(sn+k))

′ no tempo t sao obtidas

atraves de E[Znot |Zt].

2.3.2 Verossimilhanca

Dado o conjunto de observacoes Z = (Z ′1, . . . ,Z′T )′ e Θ = (σ,φ,h) o conjunto de

todos os parametros do modelo, a funcao de verossimilhanca e dada por,

p(Z | Θ) ∝T∏t=1

| Ω |−12 exp

−1

2Z ′tΩ

−1Zt

(2.6)

∝ | Ω |−T2 exp

−1

2

T∑t=1

Z ′tΩ−1Zt

.

2.3.3 Distribuicoes a Priori e a Posteriori

A distribuicao a priori dos parametros e assumida na forma


p(Θ) = p(σ)p(φ)p(h)

=M∏m=1

p(σm)p(φm)p(hm).

Embora pareca razoavel, na ausencia de informacao a priori relevante, assumir pri-

oris Gama e Gama-Invertida pouco informativas para aos parametros de suavizacao,

variancias e amplitudes, estudos pilotos simulados mostraram que estas distribuicoes

tiveram impacto na posteriori desses parametros. Berger, Oliveira, e Sanso (2001) con-

sideram o uso de prioris de referencia para estes parametros quando se tem um unico

processo Gaussiano isotropico definido para toda regiao. Em nosso caso, o exercıcio de

obter estas prioris e complicado pois estamos definindo um processo nao-estacionario

como mistura de diferentes processos Gaussianos. Alem disso, esta mistura tambem

considera parametros de suavizacao dos nucleos a serem estimadas. Algumas outras

distribuicoes alternativas foram testadas, mas bons resultados foram geralmente obti-

dos com distribuicoes a priori semi-Cauchy (half-Cauchy) sugeridas por Gelman (2006).

Se κ e um dado parametro com distribuicao semi-Cauchy, entao sua densidade e pro-

porcional a

p(κ) ∝(

1 +(κc

)2)−1

, (2.7)

se κ > 0 e 0, caso contrario.

A Figura 2.1 mostra como essa priori se comporta para diferentes valores do parametro

de escala c.

Gelman (2006) sugere o uso dessas prioris para o parametro desvio-padrao e nao

para a variancia para evitar distribuicoes a posteriori improprias.

Um cuidado especial deve ser tomado para a especificacao a priori dos parametros

de suavizacao hm.


Seja s− um = (xm, ym). Entao, para todo m,

γm(s) =km(s)∑M

m′=1 km′(s)

=

12πh2

mexp

− 1

2h2m

[x2m + y2

m]

∑Mm′=1

12πh2

m′exp

− 1

2h2m′

[x2m′ + y2

m′ ]

=

12πh2

mexp

− 1

2h2m

[x2m + y2

m]

12πh2

mexp

− 1

2h2m

[x2m + y2

m]

+ c

Segue diretamente que limhm→∞ γm(s) = 0 e que limhm→0 γm(s) = 0. Esses calculos

mostram que γm(s) → 0 quando hm → ∞ e γm(s) → 0 quando hm → 0. Assim, a

verossimilhanca perfilada de hm (que depende de hm somente atraves de Ω, que depende

de hm somente atraves de γm(s)) converge para o mesmo valor se hm converge para 0

ou infinito. Ou seja, se L(hm;Y ,Θ−hm) representa a verossimilhanca perfilada de hm,

limhm→∞

L(hm;Y ,Θ−hm) = limhm→0

L(hm;Y ,Θ−hm). (2.8)

Deste modo, os parametros de suavizacao devem ser adicionalmente restritos a um

suporte finito a fim de se verificar a possibilidade de sua convergencia para 0. Se isto

acontece, esta componente de mistura e desnecessaria e pode ser removida do modelo.

Esta observacao fornece um procedimento simples para determinacao do numero de

componentes. Um exercıcio simulado na proxima secao e outro no capıtulo 5 fornecerao

maiores evidencias da adequacao desta abordagem. Neste caso, outras alternativas

complexas para estimar o numero desconhecido de componentes sao evitadas.

Os centros das componentes de mistura sao alocados em torno da regiao de interesse

a fim de cobri-la adequadamente. Consideracoes especıficas sobre as caracterısticas

locais das sub-regioes podem ser utilizadas tambem. Elas podem, mas nao necessaria-

mente precisam ser interpretadas em associacao com os aspectos fısicos locais. De fato,

eles sao utilizados na intencao de se fazer uma aproximacao adequada da representacao

espacial infinitamente dimensional. Os centros sao assumidos fixados neste trabalho

mas podem tambem ser estimados como em Fuentes, Chaudhuri, e Holland (2007).


0 50 100 150

0.0

0.2

0.4

0.6

0.8

1.0

σσm

Den

sida

de

c=25c=10c=5

Figura 2.1: Densidades proporcionais de uma distribuicao Semi-Cauchy para tres diferentes valores

da constante c = 5, 10, 25. A distribuicao apresenta cauda pesada e sugere prioris pouco informativas.

A combinacao das especificacoes a priori acima com a verossimilhanca (2.6) fornece

a distribuicao a posteriori pelo teorema de Bayes,

p(Θ | Z) ∝ | Ω |−T2 exp

−1

2

T∑t=1

Z ′tΩ−1Zt

M∏m=1

p(σm)p(φm)p(hm). (2.9)

Como usualmente ocorre em modelos com este grau de complexidade, e difıcil ex-

trair informacoes analiticamente e algum mecanismo de aproximacao deve ser uti-

lizado. Dentre as poucas alternativas disponıveis, um algoritmo MCMC hıbrido com

amostrador de Gibbs (ver capıtulo 5), passos de Metropolis-Hastings e slice sampling foi

escolhido. No restante do texto que segue, Θ−ξ denota o vetor completo de parametros

sem seu componente ξ.

2.4 Aspectos computacionais

As condicionais completas de σm, φm e hm sao respectivamente proporcionais a,


p(σm | Θ−σm) ∝ | Ω |−T2 exp

−1

2

T∑t=1

Z ′tΩ−1Zt

p(σm). (2.10)

p(φm | Θ−φm) ∝ | Ω |−T2 exp

−1

2

T∑t=1

Z ′tΩ−1Zt

p(φm). (2.11)

p(hm | Θ−hm) ∝ | Ω |−T2 exp

−1

2

T∑t=1

Z ′tΩ−1Zt

p(hm). (2.12)

Cada um dos parametros σm, φm, m = 1, . . . ,M e amostrado individualmente por

passos de Metropolis. Por aproveitarem da mesma estrutura da funcao de verossimi-

lhanca e terem distribuicao a priori semi-Cauchy, os passos de amostragem de cada um

dos parametros σm e φm, m = 1, . . . ,M no algoritmo MCMC sao parecidos. O quadro

abaixo apresenta como exemplo, como um elemento σm pode ser amostrado.

Algoritmo: Metropolis-Hastings para σm, m = 1, . . . ,M .

A cada passo da iteracao i+ 1 do algoritmo de Metropolis-Hastings,

• Gere cada componente σ(p)m proposto, a partir de uma proposta de transicao

Gama-Invertida (GI) com media σ(i)m e variancia σ

2(i)m

∆σm,

σ(p)m ∼ GI

[∆σm + 2, σ(i)

m (∆σm + 1)].

• Aceite σ(p)m com probabilidade,

B(σ(i)m , σ

(p)m ) = min

1,p(σ

(p)m | Θ(i)

−σm)q(σ(p)m → σ

(i)m )

p(σ(i)m | Θ(i)

−σm)q(σ(i)m → σ

(p)m )

onde q(.) e a densidade da proposta de transicao e p(.) e fornecida pela equacao

(2.10).

Amostrar os parametros de suavizacao hm eficientemente requer utilizar outra es-

trategia mais especıfica. E importante ser capaz de reconhecer quando a convergencia


deles vai para 0 o que implicaria na remocao deste componente do modelo. O algoritmo

slice sampling proposto por Neal (2003) e utilizado com este fim. Agarwal e Gelfand

(2005) utilizam este algoritmo no contexto espacial e relatam as vantagens sobre os

esquemas de reamostragem e Metropolis.

O quadro abaixo apresenta como exemplo, como um elemento hm pode ser amostrado.

Algoritmo: Slice Sampling para hm, m = 1, . . . ,M .

• Amostre aleatoriamente h(1)m da distribuicao uniforme [0, R] e um ponto S da

distribuicao uniforme[0, L(Θh

(1)m ;Y )

], onde L(Θh

(1)m ;Y ) e a verossimilhanca de

h(1)m e valores correntes dos outros parametros.

• O ponto S define um slice horizontal na funcao de verossimilhanca. Um novo

ponto h(2)m e amostrado da distribuicao condicional completa de hm ate que ele

esteja dentro do slice.

Maiores detalhes tecnicos que tornam o uso deste algoritmo viavel como a utilizacao

de um esquema de amostragem com retracao (shrinkage sampling scheme) e o uso do

algoritmo na escala logarıtmica podem ser encontrados em Neal (2003) e Agarwal e

Gelfand (2005).

Nossa abordagem consiste em definir suportes iniciais [0, R] de amplitudes largas

em relacao a regiao de estudo para todos os parametros de suavizacao hm. Se um dado

parametro de suavizacao nao converge, indo para R, o algoritmo e reinicializado para

um novo e menor suporte [0, R1], ondeR1 < R, para este parametro de suavizacao. Esta

estrategia fornece bons resultados em exercıcios simulados, como descrito na proxima

secao, provendo indicacoes claras da relevancia de um dado componente adicionado

ao modelo. Ela tambem evita o uso de algoritmos de saltos reversıveis (Green, 1985)

para determinar o numero de componentes tornando o procedimento computacional

relativamente simples.


2.4.1 Interpolacao

As respostas nao-observadasZnot = (Zt(sn+1), . . . , Zt(sn+k))

′ para localizacoes (sl, . . . , sl+k)

no tempo t sao obtidas atraves de E[Znot |Zt], onde Zno

t |Zt representa a distribuicao

preditiva, cuja densidade e dada por

p(Znot |Zt) =

∫p(Zno

t ,Θ|Zt)dΘ (2.13)

=

∫p(Zno

t |Zt,Θ)p(Θ|Zt)dΘ

≈1

G

G∑g=1

p(Znot |Zt,Θ

(g))

onde p(Znot |Zt,Θ) tem distribuicao condicional normal surgindo da distribuicao normal

conjunta de Znot e Zt, e por ultimo Θ(1), . . . ,Θ(G) sao amostras da distribuicao a

posteriori p(Θ|Zt) obtidas pelo algoritmo MCMC descrito. Cabe observar tambem

que E[Znot |Zt] e V ar[Zno

t |Zt] podem ser obtidos pelas aproximacoes:

E(Znot |Zt) = E [E(Zno

t |Zt,Θ)]

≈1

G

G∑g=1

E(Znot |Zt,Θ

(g))

V ar(Znot |Zt) = E [V ar(Zno

t |Zt,Θ)] + V ar [E(Znot |Zt,Θ)]

≈1

G

G∑g=1

V ar(Znot |Zt,Θ

(g))

+

1

G

[E(Znot |Zt,Θ

(g))− E (Zno

t |Zt)] [E(Znot |Zt,Θ

(g))− E (Zno

t |Zt)]′.

2.5 Simulacoes

Dois diferentes estudos simulados sao considerados nesta secao. Estes estudos servem

para ilustrar as vantagens obtidas quando os parametros de suavizacao das com-

ponentes de mistura sao estimadas no modelo. Em ambos os estudos foram con-

sideradas n = 30 localizacoes (s1, . . . , s30) em D = [0, 100] × [0, 100] e 5 centros


(u1, u2, u3, u4, u5). Os valores dos parametros da funcao de covariancia dos proces-

sos estacionarios de media 0 foram, σ21 = 1.0, σ2

2 = 5.0, σ23 = 2, σ2

4 = 3, σ25 = 0.5 e

φ1 = 0.040, φ2 = 0.015, φ3 = 0.025, φ4 = 0.020, φ5 = 0.030). A partir desse conjunto

de parametros e valores fixados, T = 200 realizacoes independentes foram geradas,

(Z ′1, . . . ,Z′200)′, onde Zt = (Zt(s1), . . . , Zt(s30))′.

Os algoritmos de estimacao dos parametros utilizado em todos os exemplos e si-

mulacoes deste trabalho, foram implementados no programa Ox; as simulacoes e analise

dos dados feitas no R.

2.5.1 Estudo 1

Neste primeiro estudo, supomos observacoes de um processo nao-estacionario definido

a partir da mistura de cinco processos estacionarios. Entretanto, para um destes pro-

cessos, a componente de mistura e gerada a partir de parametro de suavizacao proximo

de 0. Isto significa na realidade que apenas quatro destes processos sao utilizados para

gerar os dados. Nossos interesses neste estudo sao: 1) mostrar como o procedimento de

amostragem para o modelo proposto e capaz de diferenciar o parametro de suavizacao

proximo de 0. Neste caso, isso implica que esta componente se faz desnecessaria no

modelo e deve ser removida; 2) Comparar a capacidade preditiva em termos da media

e variancia do modelo proposto (modelo I) com o modelo de misturas de processos

localmente estacionarios (modelo II) considerando os mesmos cinco centros e com o

modelo de misturas de processos localmente estacionarios (modelo III) considerando

os quatro centros que de fato possuem pesos associados.

Os dados foram gerados a partir dos parametros de suavizacao h1 = 0.01, h2 =

30, h3 = 20, h4 = 40, h5 = 15. A Figura 2.2 apresenta as curvas de nıvel das compo-

nentes de mistura γ1(.), γ2(.), γ3(.), γ4(.) e γ5(.) que sao observadas ao considerarmos

estas componentes para o vetor h. Esta especificacao faz com que a correlacao obser-

vada nos dados varie de 0.004 a 0.7989.

Modelo I

No algoritmo de estimacao foi inicialmente assumido um suporte relativamente largo


0 20 40 60 80 100

020

4060

8010

0

γγ1

longitude

latit

ude

u1

u2 u3

u4 u5

0 20 40 60 80 100

020

4060

8010

0

γγ2

longitude

latit

ude

u1

u2 u3

u4 u5

0.2 0.4 0.6

0.8

0 20 40 60 80 100

020

4060

8010

0

γγ3

longitude

latit

ude

u1

u2 u3

u4 u5

0.2

0.4

0.6

0.8

0 20 40 60 80 100

020

4060

8010

0

γγ4

longitude

latit

ude

u1

u2 u3

u4 u5

0.2

0.4 0.6

0.8

0 20 40 60 80 100

020

4060

8010

0

γγ5

longitude

latit

ude

u1

u2 u3

u4 u5

0.2 0.4

0.6 0.8

Figura 2.2: Curvas de nıvel de cada uma das componentes de mistura γ1(.), γ2(.), γ3(.),γ4(.) e

γ5(.) em D. As curvas fornecem uma ideia da influencia de cada componente nas observacoes Zt(si).

de 0 a 100 para cada um dos 5 parametros de suavizacao hm considerados. Este e

um limite bastante elevado, uma vez que maxm hm = 40. As distribuicoes a priori

Semi-Cauchy foram assumidas com c = 30. A Figura 2.3 apresenta a verossimilhanca

perfilada para cada hm, onde a concentracao em torno de 0 somente e observada para a

componente no. 1. A Figura 2.4 apresenta os tracos para os 5 parametros de suavizacao.

Ela mostra que embora os parametros que influenciam no processo observado sejam

bem estimados a distribuicao para h1 se acumula no limite superior do intervalo. Ou-

tras exemplos simulados mostram que nem sempre os parametros de suavizacao que

influenciam o processo observado sao bem estimadas quando um destes parametros

nao converge (veja o exemplo simulado do capıtulo 5). De qualquer maneira ha uma

indicacao de problemas de convergencia e a discussao da secao anterior ja indicava que

isto poderia indicar h1 concentrado em torno de 0. O algoritmo de amostragem foi

reiniciado com um suporte mais curto [0, 15], somente para h1. A Figura 2.5 mostra

que todas as distribuicoes dos parametros de suavizacao ficam concentradas em torno

de seus verdadeiros valores, incluindo h1 que se concentra em torno de 0. Isto indica

que esta componente e desnecessaria e nao deve ser incluıda no modelo.


0 50 100 150 200

−15

50−

1450

−13

50−

1250

h1

log(

Ver

ossi

milh

ança

)

0 50 100 150 200

−15

50−

1450

−13

50−

1250

h2

log(

Ver

ossi

milh

ança

)

0 50 100 150 200

−15

00−

1400

−13

00

h3

log(

Ver

ossi

milh

ança

)

0 50 100 150 200

−15

50−

1450

−13

50−

1250

h4

log(

Ver

ossi

milh

ança

)

0 50 100 150 200

−17

00−

1600

−15

00−

1400

−13

00

h5

log(

Ver

ossi

milh

ança

)

Figura 2.3: Log das funcoes de verossimilhanca perfilada para h1, h2, h3, h4 e h5.

iteração

h 1

0 1000 2000 3000 4000 5000

4050

6070

8090

100

iteração

h 2

0 1000 2000 3000 4000 5000

3040

5060

70

iteração

h 3

0 1000 2000 3000 4000 5000

2040

6080

iteração

h 4

0 1000 2000 3000 4000 5000

3040

5060

7080

9010

0

iteração

h 5

0 1000 2000 3000 4000 5000

2030

4050

6070

80

Figura 2.4: Tracos do slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha horizontal)

considerando o suporte [0, 100] for h1.

O Algoritmo e reiniciado excluindo a componente 1. A Tabela 2.1 confirma o achado

visual da Figura 2.5 para os outros parametros do modelo. Quando os parametros de


iteração

h 1

0 1000 2000 3000 4000 5000

02

46

8

iteração

h 2

0 1000 2000 3000 4000 5000

2030

4050

6070

80

iteração

h 3

0 1000 2000 3000 4000 5000

2030

4050

6070

80

iteração

h 4

0 1000 2000 3000 4000 5000

3040

5060

7080

iteração

h 5

0 1000 2000 3000 4000 5000

2040

6080

100

Figura 2.5: Tracos do slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha horizontal)

considerando o suporte [0, 15] para h1.

suavizacao e numero de componentes sao corretamente estimados, todos os parametros

podem ser identificados e eles sao bem estimados.

Modelo II

Para obter as estimativas do modelo II cada um dos valores de hm e fixado como

igual a 1. Este valor gera automaticamente 5 processos que sao localmente esta-

cionarios. A Tabela 2.2 fornece as estimativas obtidas para σ2m e φm, m = 1, . . . , 5.

Modelo III

O modelo III difere do modelo II no fato de que somente as componentes 2, 3, 4 e

5, que de fato possuem pesos associados ao processo observado, sao consideradas. Este

novo esquema gera 4 processos que sao localmente estacionarios. A Tabela 2.3 fornece

as estimativas obtidas para σ2m e φm, m = 2, . . . , 5.

Comparacao dos modelos

As estimativas de σ2m e φm obtidas pelo modelo I (Tabela 2.1) nao podem ser

comparadas com as estimativas obtidas pelo modelo II (2.2) e modelo III (2.3). As


Tabela 2.1: Valores estimados de parametros a posteriori para o Modelo I e seus respectivos

intervalos de credibilidade de 95%.

Parametros Verdadeiros valores Medias a posteriori I.C. 95 %

σ22 5.0 4.7 (4.1, 5.2)

σ23 2.0 2.2 (2.0, 2.5)

σ24 3.0 3.1 (2.5, 3.6)

σ25 0.5 0.5 (0.4, 0.6)

φ2 0.015 0.015 (0.013, 0.018)

φ3 0.025 0.021 (0.018, 0.025)

φ4 0.020 0.019 (0.016, 0.024)

φ5 0.030 0.031 (0.026, 0.038)

Tabela 2.2: Valores estimados de parametros a posteriori para o Modelo II e seus respectivos


Parametros Medias a posteriori I.C. 95 %

σ21 1.2 (1.0, 1.4)

σ22 2.4 (2.1, 2.7)

σ23 1.6 (1.4, 1.9)

σ24 1.8 (1.6, 2.1)

σ25 0.4 (0.4, 0.5)

φ1 0.025 (0.0203, 0.030)

φ2 0.018 (0.015, 0.022)

φ3 0.023 (0.018, 0.027)

φ4 0.033 (0.028, 0.039)

φ5 0.033 (0.028, 0.039)


Tabela 2.3: Valores estimados de parametros a posteriori para o Modelo III e seus respectivos


Parametros Medias a posteriori I.C. 95 %

σ22 2.3 (2.0, 2.6)

σ23 1.6 (1.4, 1.8)

σ24 1.7 (1.5, 1.9)

σ25 0.4 (0.4, 0.5)

φ2 0.017 (0.015, 0.020)

φ3 0.024 (0.020, 0.028)

φ4 0.021 (0.018, 0.025)

φ5 0.033 (0.029, 0.038)

estimativas destes parametros, obtidas pelos modelos II e III revelam o grau de nao-

estacionaridade do processo global e podem ser interpretadas, pois cada uma delas

corresponde a uma distinta sub-regiao de estacionaridade. Os parametros estimados

pelo modelo I nao podem ser interpretados dada a misturas dos processos pelas com-

ponentes de mistura.

Uma vez que estamos tratando de um estudo simulado, conhecemos a verdadeira dis-

tribuicao condicional [Znot |Zt] e podemos comparar diretamente a esperanca e variancia

condicionais E[Znot |Zt] e V ar[Zno

t |Zt] com suas respectivas estimativas obtidas pelos

modelos I, II e III. Para isto, observacoes geradas em 256 pontos de uma grade regular

foram geradas e guardadas para comparacao.

A Figura 2.6 apresenta os verdadeiros valores de E[Znot |Zt] no tempo t = 10 e os

valores interpolados para estes pontos obtidos a partir das equacoes descritas na sub-

secao 2.4.1. E nıtido observar que as estimativas obtidas pelo Modelo I sao superiores,

as obtidas pelos modelos II e III. De fato, a soma dos resıduos quadrados sao 0.019

(modelo I), 3.733 (modelo II) e 7.539 (modelo III). O mesmo cenario ocorre para outros

diferentes valores de t.

As verdadeiras variancias V ar[Znot |Zt] sao as mesmas para todo t. A Figura 2.7


−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

Figura 2.6: Verdadeira esperanca condicional E[Zno10 |Z10] (acima e a esquerda) no tempo t = 10

e suas estimativas fornecidas pelo Modelo I (acima e a direita), Modelo II (abaixo e a esquerda) e

Modelo III (abaixo e a direita).

apresenta estes verdadeiros valores e os valores interpolados para estes pontos obtidos

no tempo t = 10. Os graficos revelam que o modelo I recupera estas variancias,

enquanto que os outros modelos fornecem estimativas que sobreestimam ou subestimam

estes valores em algumas regioes. A soma dos resıduos quadrados das variancias sao

0.02 (modelo I), 5.64 (modelo II) e 5.94 (modelo III).

Estas estimativas ruins dos valores estimados E[Znot |Zt] e V ar[Zno

t |Zt] obtidos


0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

0.0

0.2

0.4

0.6

0.8

1.0

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

Figura 2.7: Verdadeira variancia condicional V ar[Zno|Z] (acima e a esquerda) e suas estimativas

fornecidas pelo Modelo I (acima e a direita), Modelo II (abaixo e a esquerda) e Modelo III (abaixo e

a direita) no tempo t = 10.

pelos modelos II e III acabam por influenciar na capacidade de interpolacao real dos

dados. Por exemplo, considere o ponto se = (53.3, 46.7). A Figura 2.8 apresenta os

valores verdadeiros observados de Zt(se) para t = 1, . . . , 200 (em ordem crescente),

contra as estimativas da media e respectivos intervalos de credibilidade de 95%. A

qualidade dos intervalos de confianca e valores preditos obtidos pelo modelo I sao

notavelmente superiores aqueles obtidos pelos modelos II e III. O tamanho medio

dos intervalos para o Modelo I (4.39) e menor quando comparado ao tamanho dos


intervalos obtidos pelos modelos II e III (6.28, 5.45 respectivamente). Alem disso,

o erro quadratico medio, EQM =∑T

t=1(E[Zt(se)|Zt]−Zt(se))2

Te menor para o modelo I

(0.84) comparado aos modelos II e III (1.49, 1.30 respectivamente).

Este exemplo simulado mostra que quando um processo nao-estacionario subjacente

e de fato dado pela mistura de processos estacionarios latentes como no modelo (2.1),

parte da capacidade preditiva e perdida ao se considerar os modelos de misturas de

processos localmente estacionarios. Este estudo tambem revela que o mecanismo de

amostragem adotado e capaz de identificar o numero correto de misturas, bem como

quais devem ser as componentes a serem consideradas no modelo. Este esquema e

robusto no sentido de recuperar os parametros utilizados no processo de simulacao.

No capıtulo 5 dados reais de temperatura mınima mensal para o estado do Rio de

Janeiro sao modelados. Esta aplicacao mostra a importancia cientıfica de considerar

abordagens mais elaboradas para misturas de processos, envolvendo a estimacao dos

parametros de suavizacao dos nucleos.

2.5.2 Estudo 2

Neste segundo estudo, supomos observacoes de um processo nao-estacionario global

definido a partir da mistura de cinco processos localmente estacionarios. Portanto

todas os parametros de suavizacao utilizados no processo de estimacao tem seus valores

relativamente pequenos, iguais a 1. Nosso interesse neste estudo e mostrar que o

procedimento de amostragem para o modelo proposto e capaz de estimar parametros

de suavizacao suficientemente pequenos que possam identificar os processos localmente

estacionarios em suas respectivas sub-regioes. O modelo proposto sera chamado de

modelo I e o modelo de mistura de processos localmente estacionario, modelo II.

Os dados foram gerados a partir de parametros de suavizacao h1 = 1, h2 = 1, h3 =

1, h4 = 1, h5 = 1. A Figura 2.9 apresenta as curvas de nıvel dos componentes de mistura

γ1(.), γ2(.), γ3(.), γ4(.) e γ5(.) que sao observadas ao considerarmos estas componentes

para o vetor h. Esta especificacao faz com que a correlacao observada nos dados varie

de aproximadamente 0 a 0.8280.


0 50 100 150 200

−4

−2

02

4

t

inte

rval

os

0 50 100 150 200

−4

−2

02

4

t

inte

rval

os

0 50 100 150 200

−4

−2

02

4

t

inte

rval

os

Figura 2.8: Valores interpolados (-) de Zt(.) para a localizacao se e seus respectivos intervalos

de credibilidade de 95% para o modelo proposto (acima) e para o modelo de mistura de processos

localmente estacionarios (abaixo). Os verdadeiro valores simulados Zt(se) sao representados por ∆.


0 20 40 60 80 100

020

4060

8010

0

γγ1

longitude

latit

ude

u1

u2 u3

u4 u5

0.2

0.4

0.6

0.8 1

0 20 40 60 80 100

020

4060

8010

0

γγ2

longitude

latit

ude

u1

u2 u3

u4 u5

0.2 0.4

0.6

0.8 1

0 20 40 60 80 100

020

4060

8010

0

γγ3

longitude

latit

ude

u1

u2 u3

u4 u5

0.2

0.4 0.6

0.8

1

0 20 40 60 80 100

020

4060

8010

0

γγ4

longitude

latit

ude

u1

u2 u3

u4 u5

0.2 0.4 0.6 0.8 1

0 20 40 60 80 100

020

4060

8010

0

γγ5

longitude

latit

ude

u1

u2 u3

u4 u5

0.2 0.4

0.6

0.8 1

Figura 2.9: Curvas de nıvel de cada uma das componentes de mistura γ1(.), γ2(.), γ3(.),γ4(.) e

γ5(.) em D. As curvas fornecem uma ideia da influencia de cada componente e respectivo processo

estacionario m nas observacoes Zt(si).

Modelo I

No algoritmo de estimacao foi inicialmente assumido um suporte de 0 a 100 para

cada um dos 5 parametros de suavizacao hm considerados. Este e um limite bastante

elevado, uma vez que maxm hm = 1. As distribuicoes a priori semi-Cauchy foram as-

sumidas com c = 30. A Figura 2.10 mostra os tracos para os 5 parametros. Inicialmente

o algoritmo sugere que h3 e h4 estejam proximos de 0. O algoritmo e reiniciado com um

suporte [0, 3] para estes 2 parametros e os tracos passam a indicar h1 tambem proximo

de 0 (Figura 2.11). Este procedimento continua e apos certo cuidado com as taxas

de aceitacao dos outros parametros nos passos de Metropolis, chega-se a um resultado

final indicado pelos tracos da Figura 2.12. Neste caso eles indicam todas parametros de

suavizacao com amplitudes proximas de 1. Estes valores sao suficientemente pequenos

para definir processos localmente estacionarios.

Modelo II

Para obter as estimativas do modelo de misturas de processos localmente esta-


iteração

h 1

0 1000 2000 3000 4000 5000

2040

6080

100

iteração

h 2

0 1000 2000 3000 4000 5000

1020

3040

50

iteração

h 3

0 1000 2000 3000 4000 5000

5060

7080

9010

0

iteração

h 4

0 1000 2000 3000 4000 5000

2040

6080

100

iteração

h 5

0 1000 2000 3000 4000 5000

2040

6080

100

Figura 2.10: Tracos do slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha horizon-

tal) considerando o suporte [0, 100].

iteração

h 1

0 1000 2000 3000 4000 5000

6070

8090

100

iteração

h 2

0 1000 2000 3000 4000 5000

3035

4045

5055

60

iteração

h 3

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

iteração

h 4

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

iteração

h 5

0 1000 2000 3000 4000 5000

2040

6080

100

Figura 2.11: Tracos do slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha horizon-

tal) considerando o suporte [0, 3] para h3 e h4.


iteração

h 1

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

2.0

2.5

iteração

h 2

0 1000 2000 3000 4000 5000

0.5

1.0

1.5

2.0

2.5

iteração

h 3

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

2.0

2.5

3.0

iteração

h 4

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

2.0

2.5

3.0

iteração

h 5

0 1000 2000 3000 4000 5000

0.0

0.5

1.0

1.5

2.0

2.5

3.0

Figura 2.12: Tracos do slice sampling para hm, m = 1, . . . , 5 para dois diferentes pontos de partida

e verdadeiros valores (linha horizontal) considerando o suporte [0, 3].

cionarios cada um dos parametros hm teve seu valor fixado em 1. Os outros parametros

do modelo foram estimados.

Comparacao dos modelos

A Tabela 2.4 compara as estimativas obtidas para o modelo I e modelo II. Os

resultados indicam nao haver diferencas significativas nas estimativas dos parametros

σ2m e φm considerando estes modelos.

As estimativas de σ2m e φm, m = 1, . . . , 5, do modelo I ficam muito proximas das

estimativas obtidas pelo modelo II. Alem disto, os valores de cada um dos hm obtidos

pelo modelo I sao pequenos, tais que para cada si, um dos valores de γm(si) e igual a

1 para algum m e aproximadamente 0 para os outros valores de m. Por causa disto,

as predicoes obtidas pelo modelo I sao proximas das predicoes obtidas pelo modelo II,

que sao proximas dos dados simulados.

A Figura 2.13 apresenta os verdadeiros valores de E[Znot |Zt] no tempo t = 10 e os

valores interpolados para os modelos I e II. A figura revela que as estimativas obtidas

por estes modelos sao proximas da verdadeira esperanca. Alem disso, nao existem


Tabela 2.4: Valores estimados dos parametros a posteriori do estudo simulado 2 e seus respectivos

intervalos de credibilidade de 95%. Dois modelos sao comparados: o modelo proposto (modelo I) que

estima os os parametros hm e o modelo de misturas de processos localmente estacionario (modelo II)

para o qual os valores de hm foram fixados em 1.

Modelo I Modelo II

verdadeiro estimado 2.5% 97.5% estimado 2.5% 97.5%

σ21 1.0 0.8 0.7 0.9 0.8 0.8 1.0

σ22 5.0 4.7 4.1 5.6 4.6 4.1 5.2

σ23 2.0 2.0 1.8 2.3 2.0 1.8 2.4

σ24 3.0 2.7 2.4 3.1 2.8 2.4 3.3

σ25 0.5 0.5 0.4 0.5 0.5 0.4 0.5

φ1 0.040 0.045 0.038 0.052 0.046 0.038 0.053

φ2 0.015 0.015 0.013 0.018 0.015 0.013 0.018

φ3 0.025 0.025 0.020 0.028 0.024 0.020 0.029

φ4 0.020 0.021 0.018 0.026 0.021 0.017 0.024

φ5 0.030 0.032 0.028 0.037 0.031 0.026 0.036

h1 1 2.4 2.0 2.5 - - -

h2 1 2.3 2.2 2.3 - - -

h3 1 2.3 2.2 2.3 - - -

h4 1 2.3 1.6 2.9 - - -

h5 1 2.4 1.9 2.9 - - -


diferencas significativas entre as estimativas das esperancas obtidas por estes modelos.

−2

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

−2

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

−2

−1

0

1

2

3

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

Figura 2.13: Da esquerda para a direita: Verdadeira esperanca condicional E[Zno10 |Z10] no

tempo t = 10, estimativa fornecida pelo Modelo I e estimativa fornecida pelo modelo II.

A Figura 2.14 apresenta os verdadeiros valores de V ar[Znot |Zt] no tempo t = 10 e os

valores interpolados para os modelos I e II. A figura revela que as estimativas obtidas

por estes modelos sao proximas da verdadeira variancia. Alem disso, nao existem

diferencas significativas entre as estimativas das variancias obtidas por estes modelos.

0.0

0.5

1.0

1.5

2.0

2.5

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

0.0

0.5

1.0

1.5

2.0

2.5

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

0.0

0.5

1.0

1.5

2.0

2.5

0 20 40 60 80 100

0

20

40

60

80

100

longitude

latit

ude

Figura 2.14: Da esquerda para a direita: Verdadeira Variancia condicional V ar[Zno10 |Z10] no

tempo t = 10, estimativa fornecida pelo Modelo I e estimativa fornecida pelo modelo II.

Mostramos neste exemplo que o procedimento de amostragem para o modelo pro-

posto e capaz de identificar os processos localmente estacionarios e seus respectivos

parametros. Este exemplo reforca a ideia de se considerar a estimacao dos parametros

de suavizacao dos nucleos em modelos de misturas de processos estacionarios. As-

sim, nao se perde nada ao considerar estruturas mais gerais para a correlacao espacial

mesmo que a estrutura de dependencia dos dados seja mais simples.

Capıtulo 3

Modelo dinamico com estrutura de

covariancia espacial

nao-estacionaria

3.1 Introducao

Em muitas situacoes reais, observacoes que sao correlacionadas espacialmente sao ob-

servadas durante o tempo, fazendo com elas tambem estejam correlacionadas no tempo.

Neste capıtulo a estrutura de covariancias nao-estacionaria do modelo proposto (2.1)

do capıtulo anterior e inserida neste contexto. Uma abordagem de modelos dinamicos

(West e Harrison, 1997) para dados espaco-temporais que considera as observacoes

contınuas no espaco e discretas no tempo e utilizada. Diferentes modelos para dados

espaco-temporais na literatura, sao baseados em modelos dinamicos (Gelfand, Baner-

jee, e Gamerman (2005a); Paez, Gamerman, e Oliveira (2005a); Sanso, Schmidt, e

Nobre (2008); Paez, Gamerman, Landim, e Salazar (2008); Lopes, Salazar, e Gamer-

man (2008)).

Considerando convolucao de processos, a utilizacao de modelos dinamicos para da-

dos espaco-temporais univariados pode ser encontrada em Higdon (2002). O seguinte

modelo como extensao do modelo estacionario via convolucao de processos e proposto,

51


Yt(s) =M∑m=1

k(s− um)ωt(um) + µ+ εs,t

ωt(um) = ωt−1(um) + vm,t

onde ωt(um) ∼ N(0, σ2) representa um ruıdo branco no tempo t. Alem disso, εs,t e vm,t

sao erros independentes e identicamente distribuıdos, εs,t ∼ N(0, σ2ε ) e vm,t ∼ N(0, σ2

v).

Calder (2003) tambem propoe a extensao deste modelo para dados espaco-temporais

multivariados.

Utilizando uma abordagem parecida ao modelo nao-estacionario proposto por Hig-

don et al. (1999), Stroud, Muller, e Sanso (2001) propoem um modelo dinamico espaco-

temporal na qual a superfıcie media do processo e modelada como uma mistura pon-

derada local de regressores lineares. No caso, a ponderacao e feita atraves de nucleos

variando no espaco.

Partindo de uma abordagem diferente da considerada em modelos de espaco de

estados, Fuentes et al. (2005) fazem uma extensao espaco-temporal da funcao de co-

variancia dos modelos de mistura de processos localmente estacionarios (1.15). No

artigo, padroes de correntes de vento observados no tempo sao modelados a partir do

seguinte modelo hierarquico,

Yt(s) =k∑i=1

βi,t(s)fi,t(s) + εt(s)

βi,t(s) = ξi,t + ξi,t(s)

ξi,t = ξi,t−1 + wt

ξi,t(s) = ξi,t−1(s) + ηt(s),

onde wt e ηt(s) sao processos ruıdo branco independentes e εt(s) tem funcao de co-

variancia nao-separavel e nao-estacionaria, definida por uma mistura de processos es-

paciais Wm(.) localmente estacionarios ,

Cov[εti(si), εtj(sj)] =M∑m=1

k(si)k(sj)Cov1[Wm(si),Wm(sj)]Cov2[ti − tj],


onde Cov1(.) e uma funcao de covariancia espacial estacionaria e Cov2(.) e uma funcao

de covariancia temporal estacionaria.

3.2 Modelo Proposto

Seja Y t = (Yt(s1), . . . , Yt(sn))′ o processo observado em localizacoes s1, . . . , sn, no

tempo t, t = 1, . . . , T . Vamos considerar o seguinte argumento heurıstico para o

modelo dinamico a ser proposto. Seja o modelo espacial geral proposto em (1.1) sem

o erro de medida,

Y t = µt +Zt,

de modo que o processo observado no tempo t seja decomposto em uma media que

evolui no tempo, mais a realizacao de um campo aleatorio espacial Zt de media 0

que evolui, ou nao, no tempo. Supondo que Zt nao evolua no tempo, seja Wm =

(Wm(s1) . . . ,Wm(sn))′ o vetor com as realizacoes de um m-esimo processo estacionario

Wm,t(.), m = 1, . . . ,M, em s1, . . . , sn e Am = diag(γm(s1 − um), . . . , γm(sn − um)),

a matriz diagonal formada pelas componentes de mistura relacionadas a este m-esimo

processo. Suponha tambem para Wm(.) uma funcao de covariancia estacionaria expo-

nencial Cov [Wm(si),Wm(sj)] = σ2m exp(−φm|si − sj|). Entao segue que,

Zt =M∑m=1

AmWm,t.

De modo que, se Wm e Wm′ sao independentes para todo m 6= m′,

E [Zt] = 0,

V ar [Y t] =∑M

m=1σ2mAmR(φm)Am,

onde [R(φm)]i,j = exp(−φm|si − sj|).

Denotando Cov [Zt] = Ω, segue que Y t ∼ Nn(µt,Ω) e este resultado e utilizado

para definir uma equacao para as observacoes como indicado em (3.1) .


Considerando que µt = F ′tθt, a equacao de evolucao do sistema (3.2) e obtida a

partir de θt e definimos o seguinte modelo dinamico espaco-temporal:

Y t = F ′tθt + ωt, ωt ∼ Nn(0,Ω) (3.1)

θt = Gtθt−1 + νt, νt ∼ Np(0,V ) (3.2)

com distribuicao a priori θ0 ∼ Np(m0, c0) e V ∼ WIn0(S0). Alem disso, definimos

θt = (θt,1, . . . , θt,p)′, F t e Gt sao matrizes de elementos conhecidos e WI representa

a distribuicao Wishart Invertida. Supomos tambem que ωt e νt sao independentes e

mutuamente independentes.

Cada elemento de Ω e dado por,

Ωi,j = Cov[Yt(si), Yt(sj)] (3.3)

=M∑m=1

γm(si)γm(sj)Cov[Wm(si),Wm(sj)], (3.4)

onde Wm(.) , m = 1, . . . ,M , sao processos estacionarios com Cov[Wm(si),Wm′(sj)] = 0

para m 6= m′. Como no capıtulo anterior, sao aqui considerados processos Gaus-

sianos com variancia σ2m e funcao de correlacao espacial de amplitude 1/φm, para

m = 1, . . . ,M .

Cada uma das componentes de mistura γm(si) e definida a partir de pesos relativos

de nucleos contınuos km(si),

γm(si) =km(si)∑M

m′=1 km′(si), (3.5)

e portanto, satisfazem a restricaoM∑m=1

γm(si) = 1. Como no capıtulo anterior iremos

considerar nucleos Gaussianos.

Comparando este modelo aos modelos descritos na sub-secao 3.1, destacamos as

seguintes diferencas. Nosso modelo se diferencia do modelo dinamico espaco-temporal


proposto por Higdon (2002) e Calder (2003) por considerar nao-estacionaridade para es-

trutura espacial dos dados. Enquanto que os modelos dinamicos espaco-temporais pro-

postos por Higdon et al. (1999) e Stroud et al. (2001) permitem a nao-estacionaridade

espacial dos dados somente atraves dos parametros de suavizacao dos nucleos que po-

dem variar espacialmente, nosso modelo considera que tanto estes parametros, quanto

os parametros dos processos estacionarios latentes podem variar espacialmente. Em-

bora Fuentes et al. (2005) proponham um modelo espaco-temporal que considera uma

funcao de covariancia espaco-temporal, nao-estacionaria espacialmente e nao-separavel,

sua abordagem nao e inserida no contexto de modelos dinamicos, conhecidos por suas

inumeras vantagens no que diz respeito a estimacao e implementacao do modelo.

3.3 Inferencia


Sejam θ = (θ1, . . . ,θT ), σ = (σ1, . . . , σM), φ = (φ1, . . . , φM), h = (h1, . . . , hM), e

defina Θ = (θ,θ0,σ,φ,h,V ) o conjunto de todos os parametros do modelo. Dado

um conjunto de observacoes Y = (Y 1, . . . ,Y T )′, temos a seguinte funcao de verossi-

milhanca,

p(Y | Θ) ∝T∏t=1

| Ω |−12 exp

−1

2(Y t − F ′tθt)′Ω−1(Y t − F ′tθt)

(3.6)

∝ | Ω |−T2 exp

−1

2

T∑t=1

(Y t − F ′tθt)′Ω−1(Y t − F ′tθt)

.

3.3.2 Distribuicao a Priori

A distribuicao a priori dos parametros e dada por,

p(Θ) =

[T∏t=1

p(θt|θt−1,V )

][M∏m=1

p(σm)p(φm)p(hm)

]p(θ0)p(V ). (3.7)


Como no capıtulo 2, prioris semi-Cauchy (Gelman, 2006) sao utilizadas para cada

um dos parametros σm, φm e hm, m = 1, . . . ,M (ver sub-secao 2.3). As prioris para

θt, t = 1, . . . , T , seguem uma estrutura hierarquica definida pela equacao (3.2). Sao

estabelecidas distribuicoes a priori pouco informativas para θ0 ∼ Np[m0, c0], com c0

grande e V ∼ WIn0(S0), com n0 pequeno.

3.3.3 Distribuicao a Posteriori

A distribuicao a posteriori e obtida combinando a funcao de verossimilhanca e in-

formacao a priori,

p(Θ | Y ) ∝ | Ω |−T2 exp

−1

2

T∑t=1

(Y t − F ′tθt)′Ω−1(Y t − F ′tθt)

(3.8)[

T∏t=1

p(θt|θt−1,V )

][M∏m=1

p(σm)p(φm)p(hm)

]p(θ0)p(V ).

Como a expressao acima nao possui forma fechada, o algoritmo MCMC e utilizado

na amostragem de todos os parametros.


O esquema de amostragem a partir do algoritmo MCMC se da a partir das condicionais

completas de todos os parametros (ou blocos de parametros). Utilizamos amostrador

de Gibbs caso o parametro apresente a condicional completa com forma conhecida ou

caso contrario, utilizamos passos de Metropolis - Hastings ou slice sampling.


Amostrando θ

Os elementos do vetor de estados latentes θ = (θ1, . . . ,θT ) sao amostrados uti-

lizando o algoritmo FFBS ( forward filtering backward sampler) proposto por Fruhwirth-

Schnater (1994) e Carter e Kohn (1994) . Dada a estrutura markoviana da equacao de

sistema de um modelo linear dinamico, esse algoritmo fornece amostras do conjunto

completo de vetores θ = (θ1, . . . ,θT ). Esse metodo de simulacao e mais eficiente que

simular estado por estado, ou seja, a partir das distribuicoes condicionais completas de

cada um dos θt, t = 1, . . . , T .

Para utilizar o FFBS, aproveitamos as equacoes de atualizacao do filtro de Kalman

(West e Harrison, 1997).

Filtro de Kalman

Seja Dt−1 = (y1, . . . ,yt−1) o conjunto de observacoes ate o tempo t− 1 e suponha

conhecida a informacao no tempo t = 0,

(θ0 |D0) ∼ N(m0,C0)

.

Deste modo as equacoes de atualizacao, a cada instante t sao:

• Posteriori em t− 1:

Para alguma media mt−1 e matriz de covariancias Ct−1,

(θt−1 |Dt−1) ∼ N(mt−1,Ct−1)

.

• Posteriori em t,

(θt |Dt−1) ∼ N(at,Rt),

com


at = Gtmt−1 e Rt = GtCt−1G′t + V .

• Previsao um passo a frente,

(yt |Dt−1) ∼ N(f t,Qt),

com

ft = F ′tat e Qt = F ′tRtF′t + Ω.

• Posteriori em t,

(θt |Dt) ∼ N(mt,Ct),

com

mt = at +At, Ct = Rt −AtQtA′t

At = RtF tQ−t 1 e et = yt − f t.

Todas as medias mt e variancias Ct obtidas pelo filtro de Kalman sao armazenadas

para serem entao utilizadas no algoritmo FFBS descrito a seguir.

Algoritmo: FFBS

• A partir do filtro de Kalman calculamos e guardamos as medias e variancias mt

e Ct, t = 1, . . . , T .

• Amostramos um vetor de (θT | DT ). Neste passo, utilizamos a equacao da

posteriori em T do filtro de Kalman.

• Recursivamente amostramos de (θt−1 | θt,DT ) para t = T − 1, . . . , 2, da sua

distribuicao condicional dada por


(θt−1 | θt,DT ) ∼ N[(G′tV

−1Gt +C−1t−1)−1(G′tV

−1θt +C−1t−1mt−1), (G′tV

−1Gt +C−1t−1)−1

].

Amostrando θ0

A distribuicao condicional completa de θ0 tem forma conhecida,

(θ0 | Θ−θ0) ∼ N[(C−1

0 +G′1V−1G1

)−1 (C−1

0 m0 +G′1V−1θ1

),(C−1

0 +G′1V−1G1

)−1],

obtida a partir do FFBS.

Amostrando os elementos σm, φm, hm,Os passos de amostragem de cada um dos parametros σm, φm e hm, m = 1, . . .M, no

algoritmo MCMC sao os mesmos que foram apresentados na sub-secao 2.3, adaptando-

se apenas a forma das condicionais completas. Por exemplo, a condicional completa

de σm e dada por,

p(σm | Θσm) ∝| Ω |−T2 exp

−1

2

T∑t=1

(yt − F ′tθt)′Ω−1(yt − F ′tθt)

p(σm) (3.9)

Como as condicionais completas desses parametros nao tem forma conhecida, passos

de Metropolis sao utilizados para σm e φm. A amostragem dos parametros de suavizacao

hm e feita por slice sampling, como descrito no capıtulo 2.

Amostrando V

A distribuicao condicional completa de V possui forma conhecida, e e tambem

Wishart Invertida, dada por

(V | Θ−V ) ∼ WIn∗0 [S∗0], (3.10)

onde n∗0 = pT + 1 e S∗0 = 1n∗0

[∑Tt=1(θt −Gtθt−1)(θt −Gtθt−1)′ + n0S0

].


3.5 Simulacao

Este exemplo simulado, considera o seguinte modelo de regressao dinamica com n = 40

e T = 400,

Y t = F ′tθt + ωt, ωt ∼ N40[0,Ω]; t = 1, . . . , 400 (3.11a)

θt = θt−1 + νt, νt ∼ N3[0,V ] (3.11b)

onde,

θt = (θt,1, θt,2, θt,3)′,

F ′t =

1 x11 x12

1 x21 x22

.... . .

...

1 xn1 xn2

, ∀t,

V =

σ2

1v 0 0

0 σ22v 0

0 0 σ23v

.

Deste modo, F ′t corresponde a uma matriz de planejamento constituıda por duas

covariaveis supostamente observadas e uma covariavel para o intercepto. As covariaveis

observadas foram geradas de duas distintas distribuicoes normais. Alem disso,Gt = I3.

No processo de simulacao foi considerada a mistura de M = 4 diferentes proces-

sos estacionarios. Os parametros considerados foram σ2 = (1.0, 2.0, 5.0, 1.5), φ =

(0.04, 0.06, 0.08, 0.04). Para obter as componentes de mistura γm(.) foram considera-

dos parametros de suavizacao h = (12, 15, 10, 20). Para simular os estados latentes, a

matriz de covariancia V = diag(0.07, 0.04, 0.09) foi considerada. As observacoes foram

geradas em 40 localizacoes aleatoriamente distribuıdos em D = [0, 100]× [0, 100].


Prioris pouco informativas foram assumidas para todos os parametros. As dis-

tribuicoes a priori para cada um dos parametros da matriz de covariancia nao-estacionaria

foram semi-Cauchy com c = 30. Para θ0 foi considerado θ0 ∼ N3[0, 100I3].

Foi assumido que V = diag(σ21v, σ

22v, σ

23v). Para o conjunto de variancias σ2

lv,

l = 1, 2, 3, duas prioris foram testadas. Inicialmente distribuicoes a priori Gama In-

vertidas independentes de media 2 e variancia infinita, σ2lv ∼ GI(2, 2) que sao condi-

cionalmente conjugadas, e distribuicoes a priori semi-Cauchy (com c = 30) para a raiz

quadrada destes parametros, σlv, l = 1, 2, 3. No caso em que foram consideradas as

prioris Gamma Invertidas, as distribuicoes condicionais completas de σ2lv, l = 1, 2, 3,

sao independentes e dadas por,

σ2lv | Θ−σ2

lv∼ GI

(α +

T

2, β +

1

2

[T∑t=1

(θt,l − θt−1,l)2

]). (3.12)

Na Tabela 3.1 sao apresentadas as estimativas e intervalos de credibilidade de 95%

dos parametros que nao evoluem no tempo, considerando as prioris Gama Invertida

e semi-Cauchy para os elementos de V . Podemos observar que ao utilizar a priori

Gama Invertida, as estimativas σ1v e σ2v sao sobre-estimadas. Esta situacao parece

fazer referencia a situacao na qual uma priori Gama Invertida aparentemente pouco

informativa influencia as estimativas a posteriori de parametros de variancia, quando

estas variancias sao de fato pequenas (Gelman, 2006).

Nas Figuras 3.1, 3.2 e 3.3 sao apresentadas as densidades marginais a posteriori

dos conjuntos de parametros que constituem a matriz de covariancia espacial nao-

estacionaria, quando sao consideradas as prioris semi-Cauchy para os parametros σlv,

l = 1, 2, 3. As figuras tambem mostram os verdadeiros valores utilizados no processo

de simulacao. Essas figuras e a tabela 3.1 mostram que todos os parametros sao

relativamente bem estimados. Todos os intervalos de credibilidade de 95% contiveram

os parametros utilizados no processo de simulacao.

A Figura 3.4 descreve o comportamento a posteriori dos coeficientes de regressao

ao longo do tempo. Os valores estimados desses parametros acompanham a evolucao

dos verdadeiros coeficientes.


Tabela 3.1: Valores estimados a posteriori dos parametros e seus respectivos intervalos de cre-

dibilidade de 95% utilizando prioris Gama Invertida (Modelo I) para os parametros σ2lv, l = 1, 2, 3 e

utilizando prioris Semi-Cauchy para os parametros σlv (Modelo II).

Parametro Verdadeiro Modelo I Modelo II

σ21 1.00 0.99 (0.91, 1.09) 0.95 (0.87, 1.06)

σ22 2.00 1.98 (1.83, 2.13) 2.04(1.89, 2.21)

σ23 5.00 5.07 (4.67, 5.52) 4.72 (4.35, 5.08)

σ24 1.50 1.55 (1.43, 1.69) 1.45 (1.35, 1.55)

φ1 0.040 0.040 (0.035, 0.044) 0.044 (0.038, 0.050)

φ2 0.060 0.060 (0.054, 0.067) 0.056 (0.050, 0.062)

φ3 0.080 0.076 (0.066, 0.086) 0.089 (0.080, 0.098)

φ4 0.040 0.040 (0.036, 0.045) 0.042 (0.038, 0.045)

h1 12.0 12.4 (11.8, 12.9) 11.8 (11.2, 12.3)

h2 15.0 16.2 (14.7, 17.8) 14.5 (13.5, 15.8)

h3 5.0 5.1 (4.9, 5.3) 5.0 (4.8, 5.2)

h4 20.0 20.8 (19.5, 22.0) 19.9 (18.2, 21.5)

σ21ν 0.070 0.169 (0.118, 0.252) 0.109 (0.052, 0.194)

σ22ν 0.040 0.048 (0.042, 0.055) 0.044 (0.039, 0.051)

σ23ν 0.090 0.096 (0.084, 0.111) 0.086 (0.074, 0.099)


1 2 3 4 5

02

46

8

Pos

terio

ri

σσ12 σσ2

2 σσ32σσ4

2

Figura 3.1: Densidades marginais a posteriori dos parametros σ21 , σ2

2 , σ22 e σ2

4 e verdadeiros valores

utilizados na simulacao representados pelas linhas verticais.

0.04 0.06 0.08 0.10

050

100

150

200

Pos

terio

ri

φφ1 == φφ4 φφ2 φφ3

Figura 3.2: Densidades marginais a posteriori dos parametros φ1, φ2, φ3 e φ4 e verdadeiros valores



5 10 15 20

01

23

4

Pos

terio

ri

h1 h2h3 h4

Figura 3.3: Densidades marginais a posteriori dos parametros h1, h2, h3 e h4 e verdadeiros valores



tempo

θθ 1

0 100 200 300 400

1618

2022

2426

tempo

θθ 2

0 100 200 300 400

01

23

45

6

tempo

θθ 3

0 100 200 300 400

−6

−4

−2

0

Figura 3.4: Evolucao do coeficiente de regressao θt,1, t = 1, . . . , 400, (grafico de cima) ao longo

do tempo (linha cheia escura), estimativas a posteriori (linha azul cheia) e intervalos de credibilidade

de 95% (linhas azuis tracejadas). Os outros dois graficos apresentam os mesmos resultados para θt,2

(grafico do meio) e para θt,3 (grafico de baixo). Os intervalos de credibilidade para θt,2 e θt,3 nao

podem ser vistos nas figuras pois sao demasiadamente pequenos e proximos dos verdadeiros valores.

Capıtulo 4

Modelo dinamico com estrutura de

covariancia espacial dinamica

nao-estacionaria

4.1 Introducao

O objetivo deste capıtulo e permitir que a estrutura de covariancia espacial nao-

estacionaria do modelo dinamico espaco-temporal (3.1) e (3.2) possa mudar ao longo

do tempo. Do ponto de vista das aplicacoes, este modelo se justifica diante a possi-

bilidade de que fenomenos fısicos que ocorram paralelamente ao processo estudado no

tempo possam provocar mudancas na estrutura de covariancia espacial dos processos

estudados.

Como motivacao para o que se pretende, considere o seguinte exemplo. Para M = 4

componentes de misturas e t = 1, . . . , 8 intervalos discretos de tempo, seja Zt(si) =∑4m=1 γm,t(si)Wm(si), onde cada um dos Wm(.) representa um processo estacionario de

media 0 e diferente funcao de covariancia exponencial e cada um dos γm,t(.) representa

uma componente de mistura que evolui de forma diferente no tempo. Observacoes

sao simuladas no tempo t = 1 a partir dos valores gerados Wm(si) e γm,1(si), m =

1, . . . , 4. No tempo, t = 2, os valores observados de Wm(si), sao mantidos fixos, mas

66


as componentes de mistura evoluem para γm,2(si) segundo alguma funcao. O processo

evolui em t = 3, . . . , 8, de forma analoga a evolucao no tempo 2. A Figura 4.2 mostra

como os processos poderiam entao evoluir na medida em que γm,t(.) evoluısse no tempo

atraves de hm,t (Figura 4.1).

Figura 4.1: Evolucao dos parametros de suavizacao hm,t, m = 1, . . . , 4, nos tempos t = 1, . . . , 8.

4.2 Modelo Geral

Seja Y t = (Yt(s1), . . . , Yt(sn))′ o processo observado nas localizacoes s1, . . . , sn no

tempo t, t = 1, . . . , T . O modelo dinamico espaco-temporal de estrutura espacial nao-

estacionaria evoluindo no tempo e dado pelo conjunto de equacoes (4.1), (4.2), (4.3) e

(4.4) descritas a seguir.


Figura 4.2: Evolucao dos processos Zt(si) na medida em que os parametros de suavizacao hm,t

variam em t = 1, . . . , 8. Os valores de hm,t estao representados como raios das circunferencias indi-

cadas.

Y t = F ′tθt + ωt, ωt ∼ Nn[0,Ωt] (4.1)

θt = Gtθt−1 + νt, νt ∼ Np[0,V ] (4.2)

e distribuicao a priori θ0 ∼ Np(m0, c0),V ∼ WIn0(S0).

Como definido no capıtulo 3, θt = (θt,1, . . . , θt,p)′ e um vetor de estados, F t e Gt sao

matrizes constituıdas de elementos conhecidos e WI representa a distribuicao Wishart

Invertida. E suposto que ωt e νt sejam independentes e mutuamente independentes.

Cada elemento de Ωt e dado por,

Ωi,j,t = Cov[Yt(si), Yt(sj)] (4.3)

=M∑m=1

γm,t(si)γm,t(sj)Cov[Wm(si),Wm(sj)],

ondeWm(.) , m = 1, . . . ,M , sao processos estacionarios tais que Cov[Wm(si),Wm′(sj)] =

0 para m 6= m′.


Cada uma das componentes de mistura γm,t(s) e definida a partir de pesos relativos

de nucleos contınuos km,t(s) que mudam ao longo do tempo atraves da relacao

γm,t(s) =km,t(s)∑M

m′=1 km′,t(s),

e portanto, satisfazem a restricaoM∑m=1

γm,t(s) = 1. Alem disso, os nucleos km,t(s)

satisfazem a km,t(s) = k(s − um;hm,t), isto e, a variacao desses nucleos ao longo do

tempo e causada apenas pela mudanca de hm ao longo do tempo.

A estrutura dinamica em hm,t e construıda da seguinte forma. Seja ht = (h1,t, . . . , hM,t),

os M parametros de suavizacao considerados no tempo t. Assumindo que hm,t e hm′,t

sejam independentes para todo m 6= m′, ht segue uma distribuicao log-normal multi-

variada:

(ht | ht−1) ∼ logN(τ t−1,ψt−1) (4.4)

onde τ t−1 = (τ1,t−1, . . . , τM,t−1) e ψt−1 = diag(ψ1,t−1, . . . , ψM,t−1). Cada um dos

parametros ψm,t−1 e τm,t−1 sao definidos de modo que cada um dos parametros de

suavizacao hm,t tenha distribuicao log-normal com E [hm,t|hm,t−1] = hm,t−1 e variancia

constante no tempo para cada componente m, V ar [hm,t|hm,t−1] = λ2h,m. Para isso,

τm,t−1 e ψm,t−1 sao parametrizados da seguinte forma:

ψm,t−1 = log

(1 +

λ2h,m

h2m,t−1

), (4.5)

τm,t−1 = log(hm,t−1)− 1

2ψm,t−1. (4.6)

Modelo de parametros de suavizacao constantes no tempo

O modelo dinamico descrito no Capıtulo 3, segue diretamente, fazendo λ2h,m = 0 para

todo m.

Neste caso, segue que V ar [hm,t|hm,t−1] = 0 e portanto, ψm,t−1 = 0 e τm,t−1 =

log(hm,t−1) de modo que a distribuicao log-normal se degenera a hm,t = hm,t−1 = hm,


para todo t.

Modelo com parametros de suavizacao variando em intervalos de tempos

maiores

Embora o modelo geral descrito na secao 4.2 seja teoricamente plausıvel, nossas con-

clusoes obtidas a partir de exercıcios de simulacao mostram que ao considerar os

parametros de suavizacao evoluindo concomitantemente com os dados, conduz a pro-

blemas de identificabilidade, que interfere na estrutura da posteriori dos parametros,

implicando em problemas de convergencia.

Deste modo, vamos assumir que os parametros de suavizacao evoluam em intervalos

de tempos discretos maiores do que evoluam as observacoes Y t, t = 1, . . . , T . Para fa-

cilitar as expressoes que virao a seguir, suponha que o intervalo de tempo t = 1, . . . , T

possa ser particionado em B sub-intervalos de mesmo tamanho l. De fato, o inter-

valo t = 1, . . . , T tambem poderia ser particionado em B sub-intervalos de tamanhos

quaisquer bastando fazer pequenas modificacoes nas equacoes que vem a seguir.

Vamos supor que para cada tempo b′, b = 1, . . . , B, hm,l(b−1)+1 = . . . = hm,lb e

denotar este conjunto de parametros de suavizacao que sao iguais no tempo b′, por

hm,b′ . Por exemplo, hm,1′ representa os primeiros l parametros de suavizacao que sao

iguais, hm,1 = . . . = hm,l = hm,1′ , na unidade de medida do tempo t. Feitas estas

consideracoes, segue que Ωl(b−1)+1 = . . . = Ωlb = Ωb′ .

A estrutura dinamica em hm,b′ e construıda analogamente aquela descrita na secao

4.2. Seja hb′ = (h1,b′ , . . . , hM,b′), as M componentes de mistura no tempo b′. Assumindo

que hm,b′ e hm′,b′ sao independentes para todo m, hb′ segue uma distribuicao log-normal

multivariada:

(hb′ | h(b−1)′) ∼ logN(τ (b−1)′ ,ψ(b−1)′) (4.7)

onde τ (b−1)′ = (τ1,(b−1)′ , . . . , τM,(b−1)′) e ψ(b−1)′ = diag(ψ1,(b−1)′ , . . . , ψM,(b−1)′). Os

parametros ψm,(b−1)′ e τm,(b−1)′ sao definidos de modo que cada um dos parametros

de suavizacao hm,b′ tenha distribuicao log-normal com E[hm,b′ |hm,(b−1)′

]= hm,(b−1)′ e


variancia constante no tempo dentro de cada sub-regiao m, V ar[hm,b′ |hm,(b−1)′

]= λ2

h,m.

Para isso, τm,(b−1)′ e ψm,(b−1)′ sao parametrizados como:

ψm,(b−1)′ = log

(1 +

λ2h,m

h2m,(b−1)′

), (4.8)

τm,(b−1)′ = log(hm,(b−1)′)−1

2ψm,(b−1)′ . (4.9)

O procedimento de inferencia e aspectos computacionais descritos a seguir levam

em consideracao o modelo descrito acima que considera os parametros de suavizacao

variando em intervalos de tempos maiores.

As seguintes suposicoes a respeito do nucleo e funcao de covariancia sao assumidas.

Vamos considerar o seguinte nucleo Gaussiano que pode evoluir no tempo a partir de

hm,t, dado por

km,t(si) =1

2π|Λm,t|−

12 exp

−1

2(si − um)′Λ−1

m,t(si − um)

,

onde Λm,t =

h2m,t 0

0 h2m,t

.

Vamos supor que os processos estacionarios de mistura sejam isotropicos e Gaus-

sianos com variancia σ2m e funcao de correlacao espacial de amplitude 1/φm, para

m = 1, . . . ,M ,

Wm(.) ∼ PG(0, σ2

m exp(−φm|si − sj|)), (4.10)

de modo que a funcao de covariancia do modelo proposto (4.3) e dada por,

Cov[Z(si), Z(sj)] =∑M

m=1 γm,t(si)γm,t(sj)σ2m exp(−φm|si − sj|).

Neste caso, podemos utilizar a seguinte notacao matricial para representar a matriz

Ωt,

Ωt =∑M

m=1σ2mAm,tR(φm)Am,t, (4.11)

onde Am,t = diag(γm,t(s1−um), . . . , γm,t(sn−um)) e [R(φm)]i,j = exp(−φm|si− sj|).


4.3 Inferencia


Sejam θ = (θ1, . . . ,θT ), σ = (σ1, . . . , σM), φ = (φ1, . . . , φM), h = (h2′ , . . . ,hB′),

λh = (λh,1, . . . , λh,M) e defina Θ = (θ,θ0,σ,φ,h,h1′ ,λh,V ) o conjunto de todos os

parametros do modelo. Dado um conjunto de observacoes realizadas Y = (Y 1, . . . ,Y T )′,

a funcao de verossimilhanca e dada por,

p(Y | Θ) ∝T∏t=1

| Ωt |−12 exp

−1

2(Y t − F ′tθt)′Ω−1

t (Y t − F ′tθt). (4.12)

Desde que Ωl(b−1)+1 = . . . = Ωlb = Ωb′ para b = 1, . . . , B, esta verossimilhanca pode

ser reescrita como,

p(Y | Θ) ∝B∏b=1

| Ωb′ |−l2 exp

−1

2

T∑t=1

(Y t − F ′tθt)′Ω−1t (Y t − F ′tθt)

. (4.13)

4.3.2 Distribuicao a Priori

A distribuicao a priori dos parametros e dada por,

p(Θ) =

[T∏t=1

p(θt|θt−1,V )

][M∏m=1

p(σm)p(φm)p(hm)

][B∏b=2

p(hb′|h(b−1)′)

][

M∏m=1

p(λh,m)

]p(θ0)p(V )p(h1′). (4.14)

Distribuicoes a priori semi-Cauchy (Gelman, 2006) pouco informativas sao uti-

lizadas para cada um dos parametros σm, φm, hm e λh,m(ver sub-secao 2.3). As prioris

para θt seguem uma estrutura hierarquica descrita pela equacao (4.2). Sao estabele-

cidas tambem, prioris pouco informativas para θ0 ∼ Np[m0, c0], V ∼ Wn0(S0) e

h1′ ∼ logN [τ 0′ ,ψ0′ ].


4.3.3 Distribuicao a Posteriori

A distribuicao a posteriori e obtida combinando a funcao de verossimilhanca e in-

formacao a priori,

p(Θ | Y ) ∝B∏b=1

| Ωb′ |−l2 exp

−1

2

T∑t=1

(Y t − F ′tθt)′Ω−1t (yt − F ′tθt)

[

T∏t=1

p(θt|θt−1,V )

][M∏m=1

p(σm)p(φm)p(hm)

][B∏b=2

p(hb′|h(b−1)′)

][M∏m=1

p(λh,m)

]p(θ0)p(V )p(h1′). (4.15)

Como a expressao acima nao possui forma fechada o algoritmo MCMC e utilizado

para amostrar todos os parametros.


O esquema de amostragem a partir do algoritmo MCMC se da a partir das condicionas

completas de todos os parametros (ou blocos de parametros). Utilizamos o amostrador

de Gibbs caso o parametro apresente a condicional completa com forma conhecida ou

caso contrario, utilizamos passos de Metropolis.

Amostrando θ e θ0

As componente de θ sao amostrados utilizando o mesmo algoritmo FFBS descrito

na sub-secao 3.4, substituindo Ω por Ωt. O mesmo vale para θ0.

Amostrando os elementos hb′, b = 2, . . . , B − 1

A condicional completa do vetor de parametros de suavizacao hb′ = (h1,b′ , . . . , hM,b′),

para b = 2, . . . , B − 1 e dada por,


p(hb′ | Θ−hb′ ) ∝ | Ωb′ |−l2

l∏q=1

exp

−1

2(ybl−q+1 − F ′tθbl−q+1)′Ω−1

j′ (ybl−q+1 − F ′tθbl−q+1)

p(h(b+1)′ |hb′)p(hb′ |h(b−1)′). (4.16)

Desde que (hb′ | h(b−1)′) segue a distribuicao log-normal descrita em (4.4),

p(h(b+1)′|hb′) ∝M∏m=1

ψ−0.5m,b′ exp

−1

2(g(b+1)′ − τ b′)′ψ−1

b′ (g(b+1)′ − τ b′),

p(hb′ |h(b−1)′) ∝M∏m=1

h−1m,b′ exp

−1

2(gb′ − τ (b−1)′)

′ψ−1(b−1)′(gb′ − τ (b−1)′)

,

onde gb′ = log(hb′).

Como a condicional completa nao tem forma conhecida, utilizamos passos de Metropolis-

Hastings para amostrar cada bloco hb′ . Alguns cuidados devem ser tomados na escolha

da distribuicao proposta. O algoritmo e descrito a seguir.

Algoritmo: Metropolis-Hastings para hb′ .

A cada passo da iteracao i+ 1 do algoritmo de Metropolis-Hastings,

• Gere cada elemento h(p)m,b′ do vetor h

(p)b′ proposto, a partir de uma proposta de

transicao log-normal com media h(i)m,b′ e variancia ∆h constante,

h(p)m,b′ ∼ logN

[log(h

(i)m,b′)−

1

2log

(1 +

∆h

(h(i)m,b′)

2

), log

(1 +

∆h

(h(i)m,b′)

2

)].

• Aceite h(p)b′ com probabilidade,

B(h(i)b′ ,h

(p)b′ ) = min

1,p(h

(p)b′ | Θ

(i)−hb′

)q(h(p)b′ → h

(i)b′ )

p(h(i)b′ | Θ

(i)−hb′

)q(h(i)b′ → hpb′)

onde q(.) e a densidade da proposta de transicao e p(.) e fornecido pela equacao

(4.16).


Amostrando h1′

A condicional completa de h1′ = (h1,1′ , . . . , hM,1′) e dada por,

p(h1′ | Θ−h1′) ∝ | Ω1′ |−

l2

l∏q=1

exp

−1

2(yq − F ′tθq)′Ω−1

1′ (yq − F ′tθq)

p(h2′|h1′)p(h1′). (4.17)

onde,

p(h2′ |h1′) ∝M∏m=1

ψ−0.5m,1′ exp

−1

2(g2′ − τ 1′)

′ψ−11′ (g2′ − τ 1′)

.

p(h1′) ∝M∏m=1

h−1m,1′ exp

−1

2(g1′ − τ 0′)ψ

−10 (g1′ − τ 0′)

.

O esquema de amostragem e feito atraves por Metropolis-Hastings de modo similar

ao algoritmo utilizado em hb′ , b = 2, . . . , B − 1.

Amostrando hB′

A condicional completa de hB′ = (h1,B′ , . . . , hM,B′) e dada por,

p(hB′ | Θ−hB′ ) ∝ | ΩB′ |−l2

l∏q=1

exp

−1

2(yBl−q+1 − F ′tθBl−q+1)′Ω−1

B′ (yBl−q+1 − F ′tθBl−q+1)

p(hB′ |h(B−1)′). (4.18)

onde,

p(hB′ | h(B−1)′) ∝M∏m=1

h−1m,B′ exp

−1

2(gB′ − τ (B−1)′)

′ψ−1(B−1)′(gB′ − τ (B−1)′)

.

O esquema de amostragem e feito por Metropolis-Hastings de modo similar ao

algoritmo utilizado em hb′ , b = 2, . . . , B − 1.


Amostrando os elementos σm, φm, m = 1, . . . ,M

Por aproveitarem da mesma estrutura da funcao de verossimilhanca e terem dis-

tribuicoes a priori semi-Cauchy, os passos de amostragem de cada um dos parametros

σme φm no algoritmo MCMC sao parecidos aos do Capıtulo 2.

As condicionais completas para σm e φm sao proporcionais a,

p(σm | Θσm) ∝B∏b=1

| Ωb′ |−l2

T∏t=1

exp

−1

2(yt − F ′tθt)′Ω−1

t (yt − F ′tθt)p(σm). (4.19)

p(φm | Θσm) ∝B∏b=1

| Ωb′ |−l2

T∏t=1

exp

−1

2(yt − F ′tθt)′Ω−1

t (yt − F ′tθt)p(φm). (4.20)

Amostrando V

A condicional completa de V possui forma conhecida , e e Wishart Invertida,

(V | Θ−V ) ∼ WIn∗0 [S∗0], (4.21)

onde n∗0 = pT + 1 e

S∗0 =1

n∗0

[T∑t=1

(θt −Gtθt−1)(θt −Gtθt−1)′ + n0S0

].

Amostrando λh,m, m = 1, . . . ,M

A condicional completa de cada λh,m, o desvio-padrao de (hm,1′ , . . . , hm,B′), m =

1, . . . ,M e proporcional a,

p(λh,m | Θ−λh,m) ∝

[B∏b=2

p(hm,b′|hm,(b−1)′)

]p(λh,m).

Por nao ter distribuicao condicional completa conhecida, λh,m e amostrado por

Metropolis-Hastings.


4.5 Simulacao

Neste exercıcio de simulacao, um modelo dinamico de primeira ordem e considerado

fixando n = 100, T = 1200, B = 15 e M = 4,

Y t = θt + ωt, ωt ∼ N100[0,Ωt] (4.22)

θt = θt−1 + νt, νt ∼ N100[0, σ2νIn],

onde Ωi,j,t = Cov[Yt(si), Yt(sj)] =∑4

m=1 γm,t(si)γm′,t(sj)σ2m exp(−φm|si−sj|) e hm,l(b−1)+1 =

. . . = hm,lb = hm,b′ , para b = 1, . . . , 15.

Alem disso,

(hm,b′ | hm,(b−1)′) ∼ logN(τm,(b−1)′ , ψm,(b−1)′). (4.23)

onde ψm,(b−1)′ e τm,(b−1)′ sao como definidos em (4.5) e (4.6).

Note que como T = 1200 e B = 15 estamos assumindo que os parametros de

suavizacao estao evoluindo a cada 80 unidades de tempo dos dados.

Os valores considerados para os parametros foram σ2 = (1.0, 2.0, 5.0, 1.5), φ =

(0.04, 0.06, 0.08, 0.04). Para simplificar, foi considerada uma mesma variancia na evolucao

dos parametros de suavizacao, λ2h,1, . . . , λ

2h,4 = λ2

h igual a 0.5. Os valores de ht foram

gerados a partir de um vetor inicial (12, 15, 10, 20). O valor de σ2ν foi fixado em 0.5.

Fixados os parametros e valores acima, observacoes foram geradas em 100 loca-

lizacoes em D = [0, 100]× [0, 100] ao longo de T = 1200 unidades de tempo.

Foram assumidas prioris pouco informativas, θ0 ∼ Nn[0, 100I1200] ,σ2ν ∼ GI(2, 2)

e h1′ ∼ logN [1, diag(100)]. Como neste caso, V = σ2νIn, a distribuicao a condicional

completa de σ2ν e:

σ2ν | Θ−σ2

ν∼ GI

(α +

nT

2,1

2

[T∑t=1

(θt − θt−1)′(θt − θt−1) + β

]).

Na Tabela 4.1 sao apresentadas as estimativas e intervalos de credibilidade de 95%

dos parametros que nao evoluem no tempo.


Tabela 4.1: Valores estimados a posteriori dos parametros que nao evoluem no tempo e seus

respectivos intervalos de credibilidade de 95% .

Intervalo de credibilidade

Parametro Verdadeiro Estimado 2.5% 97.5%

σ21 1.00 1.03 0.99 1.07

σ22 2.00 2.01 1.93 2.09

σ23 5.00 5.07 4.94 5.23

σ24 1.50 1.54 1.49 1.60

φ1 0.040 0.038 0.036 0.040

φ2 0.060 0.060 0.057 0.063

φ3 0.080 0.080 0.078 0.083

φ4 0.040 0.039 0.037 0.041

λ2h 0.50 0.55 0.32 0.97

σ2ν 0.50 0.50 0.48 0.51

Nas Figuras 4.3 e 4.4 sao apresentadas as densidades marginais a posteriori dos

conjuntos de parametros σ2 = (σ21, . . . , σ

24) e φ = (φ1, . . . , φ4) e os verdadeiros valores

utilizados no processo de simulacao. Essas figuras e a tabela 4.1 mostram que todos os

parametros sao bem estimados. Todos os intervalos de credibilidade de 95% contiveram

o verdadeiros valores dos respectivos parametros utilizados na simulacao.

A Figura 4.5 descreve a evolucao de θt,10, t = 1, . . . 1200 e os intervalos de credi-

bilidade de 95% desses parametros. Os intervalos acompanham a verdadeira evolucao

temporal deste parametro. O mesmo ocorre para outros parametros θt,i, i 6= 10.

A Figura 4.6 apresenta a evolucao de cada um dos parametros de suavizacao hm,b′ ,

b′ = 1, . . . , 15. Cada um dos 4 parametros estimados a posteriori acompanham as

evolucoes temporais de seus respectivos verdadeiros valores. Os verdadeiros valores

hm,b′ utilizados na simulacao dos dados estao contidas nos respectivos intervalos de

credibilidade de 95%.


1 2 3 4 5

05

1015

20

Pos

terio

ri

σσ12 σσ2

2 σσ32σσ4

2

Figura 4.3: Densidades marginais a posteriori dos parametros σ21 , σ2

2 , σ22 e σ2

4 e verdadeiros valores


0.04 0.05 0.06 0.07 0.08

010

020

030

040

0

Pos

terio

ri

φφ1 == φφ4 φφ2 φφ3

Figura 4.4: Densidades marginais a posteriori dos parametros φ1, φ2, φ3 e φ4 e verdadeiros valores



0 200 400 600 800 1000 1200

−50

510

1520

25

tempo

θθ 10

Figura 4.5: Evolucao de θt,10 (pontos), para t = 1, . . . , 1200 e intervalo de credibilidade de 95%

desses parametros a posteriori (linha tracejada azul).

2 4 6 8 10 12 14

1011

1213

14

tempo

h 1

2 4 6 8 10 12 14

1113

15

tempo

h 2

2 4 6 8 10 12 14

910

1214

tempo

h 3

2 4 6 8 10 12 14

1820

2224

tempo

h 4

Figura 4.6: Evolucoes de hm,b′ , m = 1, . . . , 4 e b = 1, . . . , 15 simulados (linha cheia escura), valores

estimados (linha cheia azul) e intervalos de credibilidade de 95% (linha tracejada azul).

Capıtulo 5

Aplicacao

5.1 Introducao

As alteracoes climaticas tem sido tema de grande interesse publico e muitas linhas de

investigacao em todo o mundo, em diferentes areas da Ciencia. O ultimo relatorio

divulgado pela instituicao ganhadora do Nobel International Panel Climate Change

(IPCC, 2007) indica um aumento medio no aquecimento global e uma projecao de

aproximadamente 0.2 graus Celsius por decada. De fato, ha evidencias conflitantes no

que diz respeito deste controverso topico.

Embora o European Climate Support Network (ECSN, 1995) relate um aumento

da temperatura no ultimo seculo na maior parte da Europa, ha evidencias conflitantes

sobre este assunto. Por exemplo, Sahsamanoglou e Makrogiannis (1991) analisam

as diferencas de tendencias de temperaturas ao longo da regiao do Mediterraneo nos

ultimos 40 anos. A diversidade fisiografica e os inumeros efeitos sofridos pela regiao

levam a uma diferenciacao dos efeitos ocorridos no leste e oeste da regiao. Um efeito

resfriamento ou falta de qualquer tendencia e observada no leste, enquanto um efeito

de aquecimento e observado na regiao oeste.

Esta variabilidade e esperada e e principalmente devida a variabilidade climatica

natural, alem da nebulosidade, gases de efeito estufa, aerossois e ciclos de carbono

troposfericos, bem como fatores locais, tais como o crescimento urbano, irrigacao e

81


desertificacao. A presenca de importantes caracterısticas locais e as mudancas ocorridas

tornam relevantes o estudo das alteracoes climaticas ao nıvel regional.

O estado do Rio de Janeiro possui uma area de 43.653 km2 e uma ampla diversidade

fısica, incluindo montanhas, planıcies e alguns declives acentuados nas proximidades do

litoral. Essa diversidade contribui para um clima variando de tropical a mesotermal,

com perıodos que variam de secos a umidos. A dispersao populacional e altamente

irregular, de muito alta densidade em torno da baıa de Guanabara onde reside sua

capital, ate uma baixa densidade em suas regioes centro-sul e nordeste.

O objetivo deste trabalho e analisar as mudancas climaticas que ocorreram no

estado a partir de dados de temperatura coletados por estacoes de monitoramento em

todo o estado durante as ultimas decadas. O modelo estatıstico e baseado em uma

especificacao hierarquica na qual os dados sao modelados atraves de uma estrutura de

regressao para a media e uma representacao espacial nao-estacionaria para o termo de

ruıdo espacialmente estruturado. Ao termo de regressao tambem contem dependencia

espacial atraves de seus coeficientes e suas covariaveis.

5.2 Modelo para variacoes na temperatura

Os dados se constituem das temperaturas mınimas mensais registradas entre Janeiro de

1961 a Dezembro de 2000, totalizando T = 480 instantes de tempo, em n = 37 locais

no estado do Rio de Janeiro, como apresentado na Figura 5.1. Os dados sao reco-

lhidos por diferentes instituicoes governamentais, ligadas a rede do Instituto Nacional

de Meteorologia (INMET) e ao Centro de Previsao de Tempo e Estudos Climaticos

(CPTEC) do Instituto Nacional de Pesquisas Espaciais (INPE). Cada valor observado

de temperatura yt(si) no local si, i = 1, . . . , n e tempo t, t = 1, . . . , T , e explicado como

a soma de uma media dependente do espaco e do tempo e uma componente espacial

nao-estacionaria

yt(si) = µt(si) + εt(si). (5.1)


Formulacoes alternativas podem tambem ser consideradas. Por exemplo, a inclusao

de termos de erros nao estruturados com um erro de medida e/ou efeito pepita ou

correlacao temporal entre os erros.

A estrutura media e dada por uma forma linear com variaveis explicativas que

podem mudar no tempo e espaco

µt(s) =∑p

βp(s, t)fp(s, t), (5.2)

onde fp(s, t) sao variaveis relevantes para temperatura e βp(si, t) sao seus respectivos

coeficientes. Os efeitos dos valores das covariaveis podem mudar no espaco, como em

Gelfand, Kim, Sirmans, e Banerjee (2003) e Paez, Gamerman, e de Oliveira (2005b), ou

no espaco e tempo, como em Huerta, Sanso, e Stroud (2004) e Gelfand et al. (2005a).

As covariaveis consideradas nesta aplicacao foram a altitude, um harmonico de

perıodo 12 para representar o ciclo sazonal da temperatura em um ano e o tempo.

A mudanca climatica no tempo e determinada por caracterısticas globais e locais.

Portanto, poderia ser inadequado assumir ela como fixa no espaco. O efeito do tempo

βt e assumido ser espacialmente estruturado e modelado de acordo com um processo

Gaussiano isotropico para permitir aumentos/decrescimentos diferentes em cada local.

Assim, o vetor de coeficientes do efeito do tempo β = (β1, . . . , βn) e estruturado como

β ∼ N[b1, σ2

β exp (−φβ|si − sj|)]. (5.3)

A media global b representa a tendencia global para toda regiao, mas cada localizacao

pode ter um efeito diferente do tempo em sua temperatura. Os hiperparametros σ2β e

φβ controlam a incerteza sobre a regiao de interesse a respeito da tendencia media b

dos efeitos especıficos dos locais e a similaridade espacial entre eles, respectivamente.

5.2.1 Modelando a estrutura nao-estacionaria

O modelo de misturas apresentado nos capıtulos anteriores e considerado aqui. Como

definido no Capıtulo 2, sao considerados M diferentes centros u1, . . . ,uM localizados

em uma regiao D. Processos estacionarios Wm(.) de media 0, com funcao de covariancia

que dependem de parametros ηm e componentes de mistura γm(.) sao associados a


Figura 5.1: Localizacoes das estacoes de monitoramento e centros das componentes no estado do

Rio de Janeiro.

cada um dos centros um, para m = 1, ...,M . Deste modo, a componente espacial

nao-estacionaria e dada por

ε(s) =M∑m=1

γm(s)Wm(s), (5.4)

para cada localizacao s ∈ D.

Novamente e utilizada a restricaoM∑m=1

γm(si) = 1, de modo que cada componente

de mistura γm(s) representa o peso relativo de uma funcao nucleo contınua km(s), que

depende de um parametro de suavizacao hm de acordo com

γm(s− um) =km(s)∑M

m′=1 km′(s), ∀s ∈ D. (5.5)


A formulacao do modelo (5.4) implica que ε(·) e um processo espacial de media 0 e

funcao de covariancia

Cov (ε(si), ε(sj)) =M∑m=1

γm(si)γm(sj)Cov (Wm(si),Wm(sj)) . (5.6)

Nesta aplicacao os nucleos sao assumidos como Gaussianos e os processos esta-

cionarios de mistura assumidos serem isotropicos e Gaussianos com variancia σ2m e

funcao de correlacao exponencial de amplitude 1/φm, para m = 1, ...,M .

5.2.2 Inferencia dos parametros do modelo

Sejam as observacoes denotadas por Y = (Y ′1, . . . ,Y′T )′, onde Y t = (Yt(s1), . . . , Yt(sn))′,

t = 1, . . . , T . Entao, o modelo dado por (5.1), (5.2), (5.3) e (5.4) pode ser escrito na

forma matricial como

Y ∼ N [Xβ,Ω⊗ IT ]

β ∼ N [b1n,W ]

onde β =(β, β

)′, β = (β0, βx, βc, βs) e o vetor dos coeficientes de regressao fixados,

os elementos de Ω sao dados por Ωi,j =∑M

m=1 γm(si)γm(sj)σ2m exp (−φm|si − sj|) e os

elemento de W por W i,j = σ2β exp (−φβ|si − sj|), 1n e o vetor n-dimensional de 1’s e

In e a matriz identidade n-dimensional.

Denotando por σ = (σ1, . . . , σM), φ = (φ1, . . . , φM) e h = (h1, . . . , hM), o conjunto

de parametros a serem estimados e Θ = (β, b, σβ, φβ,h,σ,φ). A funcao de verossimi-

lhanca e

L(Θ;Y ) ∝T∏t=1

|Ω|−12 exp

−1

2(Y t −X tβ)′Ω−1 (Y t −X tβ)

, (5.7)

onde X t e a sub-matriz de X correspondente ao tempo t. A distribuicao a priori para

os coeficientes de regressao fixados e a media b dos coeficientes do tempo sao assumidas

independentes e vagas com variancias grandes ω e c0 respectivamente. Alem disso, e


assumido que os parametros da media mencionados acima sao independentes dos outros

parametros.

Como discutido no Capıtulo 2, distribuicoes a priori semi-Cauchy foram assumi-

das para os parametros de suavizacao, variancias e amplitude (range) que definem o

processo nao-estacionario.

A combinacao das especificacoes das prioris acima com a verossimilhanca (5.7)

fornece a distribuicao a priori via teorema de Bayes. Dada a complexidade desta ex-

pressao um algoritmo hıbrido MCMC com amostrador de Gibbs, passos de Metropolis-

Hastings e slice sampling foi utilizado.

5.2.3 Aspectos computacionais

A distribuicao condicional completa de β e dada por

[β|Θ−β,Y ] ∼ N[βp,V p

], (5.8)

onde V p =(V −1 +

∑Tt=1X

′

tΩ−1X t

)−1

, βp = V p

(V −1µβ +

∑Tt=1X

′

tΩ−1Y t

), µβ =(

b, b1n

)′e V = diag (ω,W ). A distribuicao condicional completa de b e facil-

mente obtida como [b|Θ−b,Y ] ∼ N [m1, c1], com c1 =(1′W−11 + c−1

0

)e m1 =

c1

(1′W−1β +m0c

−10

).

Cada parametro σβ, φβ, σm e φm, m = 1, . . . ,M e amostrado individualmente com

passos de Metropolis. A amostragem dos parametros de suavizacao e realizada atraves

de slice sampling (Capıtulo 2).

5.2.4 Dados faltantes e interpolacao

Uma caracterıstica comum de conjunto de dados de meio-ambiente e a ausencia de

dados em trechos de tempo. Isto pode ser causado por uma serie de razoes, incluindo

a manutencao da estacao, falha inesperada do equipamento de medicao, dentre outros.

Dados faltantes devem ser manuseados e contabilizados no procedimento de inferencia.

Para algum tempo t dado, o modelo induz uma distribuicao normal conjunta para o

vetor de dados nao-observados Y not e dados observados Y o

t com matriz de covariancia


Ω, que depende dos parametros do modelo Θ. Assim, um exercıcio simples para obter

a distribuicao condicional completa dos dados faltantes e considerar [Y not |Y o

t ,Θ] ∼N[Xno

t β + Ω12Ω−122 (Y o

t −Xotβ) ,Ω11 −Ω12Ω

−122 Ω′12

]. As matrizes Ω11, Ω12 e Ω22 sao

os blocos correspondentes de Ω. Como essa e a distribuicao e condicional completa

de Y not , ela pode ser facilmente amostrada. Assim, o dado faltante e naturalmente

incorporado ao ciclo do MCMC.

Comentarios parecidos podem ser feitos sobre a interpolacao das observacoes para

qualquer outro local da regiao de interesse. Este exercıcio de interpolacao e usualmente

referido como krigagem e fornece uma estrutura natural de passar as informacoes dos

locais monitorados para toda a regiao de interesse. Para isso no entanto, uma operacao

semelhante a krigagem deve ser realizada previamente aos coeficientes de regressao

dependentes espacialmente.

5.3 Exemplo simulado

Esta secao e utilizada para ilustrar a adequacao do modelo, do procedimento de

amostragem e tambem da selecao do numero de componentes. Dados foram simulados

em condicoes semelhantes aos dados originais de temperatura observados no estado do

Rio de Janeiro. Em particular, as localizacoes seguem o padrao das 37 estacoes de

monitoramento dos dados reais e os centros das componentes sao os mesmos utilizados

na analise dos dados reais.

Quatro centros de componentes u2, u3, u4 e u5 foram utilizados na geracao dos

dados. Um centro adicional u1, com h1 ≈ 0 que nao possui efeito no processo e

incluıdo no modelo para ser testado. Como nos exemplos simulados do Capıtulo 2, e

esperado que um procedimento inferencial apropriado possa estimar h1 proximo de 0,

confirmando a irrelevancia desta componente.

Os hiperparametros restantes foram tomados como σ22 = 1.4, σ2

3 = 35.7, σ24 = 2.9,

σ25 = 0.5, φ2 = 8.4 × 10−1, φ3 = 1.7 × 10−4, φ4 = 6.7 × 10−2, φ5 = 8.3 × 10−1 ,

h2 = 90.2, h3 = 110.5, h4 = 77.1 e h5 = 78.1. Foi tambem assumido que b = 9.6× 10−4

e β foi gerado com φβ = 10−3 e σ2β = 3.4 × 10−5. A maior parte destes valores


representam condicoes semelhantes as encontradas na analise dos dados reais. As

correlacoes observadas variam de 0.6 a 0.9. Finalmente, Y = (Y ′1, . . . ,Y′T )′ foi gerado

com T = 480 tempos e utilizando a matriz de planejamento dos dados reais X.

Na pratica, informacoes sobre os hiperparametros nao sao geralmente disponıveis.

Assim, foi inicialmente assumido um suporte grande de 0 a 250 kms para os 5 parametros

de suavizacao hm consideradas, m = 1, . . . , 5. Este e um limite bastante significativo

desde que maxm hm = 110.5. Foram assumidas distribuicoes semi-Cauchy para os

parametros com c = 30. A Figura 5.2 apresenta a verossimilhanca perfilada de cada

parametro de suavizacao hm, onde a concentracao em torno de 0 somente e observada

para a componente nao existente no. 1. A Figura 5.3 mostra os tracos dos 5 parametros

de suavizacao. Ela mostra que elas sao mal estimadas e a distribuicao de h1 se acu-

mulam no limite superior do intervalo. Esses sao claras indicativos de problemas de

convergencia e como discutido no Capıtulo 2 indicam h1 concentrado proximo de 0.

O algoritmo foi reiniciado com um suporte menor [0, 50], somente para h1. A Figura

5.4 mostra que todos os parametros de suavizacao sao agora bem estimadas com os

valores gerados se concentrando em torno de seus verdadeiros valores, incluindo h1 que

se concentra em torno de 0. Isso indica que esta componente nao precisa ser incluıda

no modelo.

A Tabela 5.1 confirma os achados visuais da Figura 5.4 para os outros hiper-

parametros do modelo. Quando os parametros de suavizacao e o numero de compo-

nentes sao corretamente estimados, todos os parametros do modelo podem ser identifi-

cados e sao bem estimados. Assim, o procedimento de estimacao e sua implementacao

computacional parecem estar robustas para serem utilizados em dados reais.

5.4 Resultados da analise

5.4.1 Descricao dos dados e modelo

Os dados consistem da temperatura mınima mensal de Janeiro de 1961 a Dezembro de

2000 (T = 480 meses) observados em n = 37 estacoes de monitoramento, localizados


0 100 200 300 400 500 600

−56

00−

5200

−48

00−

4400

h1

log(

Lile

lihoo

d)

0 100 200 300 400 500 600

−11

000

−90

00−

7000

−50

00

h2

log(

Lile

lihoo

d)

0 200 400 600 800 1000

−18

000

−14

000

−10

000

−60

00

h3

log(

Lile

lihoo

d)

0 200 400 600 800

−65

00−

6000

−55

00−

5000

−45

00

h4

log(

Lile

lihoo

d)

0 200 400 600 800

−60

00−

5500

−50

00−

4500

h5

log(

Lile

lihoo

d)

Figura 5.2: Funcoes das verossimilhancas perfiladas para h1, h2, h3, h4 e h5.

como indicado na Figura 5.1. A temperatura mınima foi escolhida por ser mais sensitiva

a mudancas ambientais que podem ocorrer em uma regiao geografica (Oke, Klysik, e

Bernhofer, 2006). Os valores de temperatura mınima somente foram considerados

quando pelo menos 5 medias diarias foram coletadas para cada mes. Caso contrario, a

medida foi registrada como dado faltante.

Cinco centro u1, u2, u3, u4 e u5 foram originalmente considerados. Eles foram

dispostos em 5 sub-regioes A1, A2, A3, A4 e A5 que se distinguem de acordo com suas

caracterısticas geograficas. A1 possui um terreno acidentado com densidade alta de

florestas; A2 e uma regiao plana densamente povoada, com vegetacao esparsa; A3 e

a regiao mais elevada com cadeias montanhas e floresta densa; A4 possui topografia

plana com lagos e lagoas enquanto que A5 e tambem uma regiao plana.

Tres modelos foram considerados e avaliados em suas performances preditivas. O

modelo geral baseado em uma mistura de processos estacionarios latentes discutidos

aqui e referido como M1. O modelo que considera processos localmente estacionarios

e referido como M2. E finalmente, o modelo mais simples com um unico processo

estacionario e referido como M3. O modelo M2 pode ser obtido como um caso especial


iteration

h 1

0 1000 2000 3000 4000 5000

100

150

200

250

iteration

h 2

0 1000 2000 3000 4000 5000

8090

100

110

120

iteration

h 3

0 1000 2000 3000 4000 5000

110

120

130

140

150

iteration

h 4

0 1000 2000 3000 4000 5000

7075

8085

9095

iteration

h 5

0 1000 2000 3000 4000 5000

7075

8085

9095

100

105

Figura 5.3: Tracos gerado pelo slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha

horizontal) considerando o suporte [0, 250] para h1.

do modelo M1 assumindo as amplitudes dos parametros de suavizacao pequenos e

conhecidos. O modelo M3 pode ser obtido como caso especial do modelo M2 assumindo

uma unica componente com parametro de suavizacao grande.

As componentes escalares de β e b tem distribuicao normal a priori com variancia

grande. Cada um dos σm, φm, hm, σβ e φβ possuem distribuicoes a priori semi-Cauchy

vagas.

5.4.2 Comparacao dos modelos

O desempenho dos modelos foi avaliado atraves de sua capacidade preditiva. Uma

possibilidade e um procedimento de validacao cruzada mas este poderia ser um processo

demorado para os modelos estruturados considerados. Outra possibilidade e a de deixar

apenas um pequeno numero de estacoes fora da analise e avaliar a performance conjunta

de previsao para eles. Este procedimento sempre suscita duvidas sobre a influencia que

a escolha de uma estacao particular possa ter sobre o resultado. Optamos por uma

estrategia mais geral realizando a previsao de 500 observacoes (cerca de 3% do total


iteration

h 1

0 1000 2000 3000 4000 5000

05

1015

2025

30

iteration

h 2

0 1000 2000 3000 4000 5000

8590

9510

010

511

0

iteration

h 3

0 1000 2000 3000 4000 5000

110

115

120

125

130

iteration

h 4

0 1000 2000 3000 4000 5000

7580

8590

95

iteration

h 5

0 1000 2000 3000 4000 5000

7580

8590

9510

010

5

Figura 5.4: Tracos gerado pelo slice sampling para hm, m = 1, . . . , 5 e verdadeiros valores (linha

horizontal) considerando o suporte[0, 50] para h1.

do conjunto de dados) amostradas aleatoriamente nos locais e tempo, que ficaram de

fora da analise.

A Tabela 5.2 apresenta o sumario de alguns criterios para comparacao dos modelos.

O modelo nao-estacionario M2 se saiu melhor que o modelo estacionario M3, como

esperado para uma regiao diversificada. Isto mostra que a hipotese de estacionaridade

espacial nao e apropriada para este estudo. Os erros preditos para o modelo M1 sao

substancialmente menores do que os erros encontrados para os outros modelos. Os

intervalos de credibilidade do modelo M1 tambem sao menores do que os dos outros

dois modelos. Isto e obtido sem sacrificar a sua cobertura, que e igualmente boa. A

Figura 5.5 mostra um sumario grafico das predicoes. Os intervalos de M1 parecem

exibir um padrao mais estavel enquanto que aqueles de M2 mostram um padrao mais

disperso de largura. Isso e devido provavelmente ao pressuposto de estacionaridade

local, que favorece a uma maior variacao nas informacoes associadas a cada local, em

contraste com a completa nao-estacionaridade.

A media a posteriori dos hiperparametros para M2 foram σ21 = 1.61, σ2

2 = 1.60,


Tabela 5.1: Sumario da estimacao para as variancias e amplitudes dos processos (media a posteriori

e intervalos de credibilidade de 95% para os dados simulados).

Parametros Valores verdadeiros Medias a posteriori I.C. 95 %

σ22 1.4 1.4 (1.3, 1.5)

σ23 35.7 35.1 (31.2, 39.7)

σ24 2.9 2.9 (2.8, 3.2)

σ25 0.5 0.5 (0.5, 0.6)

φ2 0.84 0.88 (0.70, 1.10)

φ3 1.7× 10−4 1.1× 10−4 (1.6× 10−4, 2.0× 10−4)

φ4 0.067 0.063 (0.056, 0.072)

φ5 0.83 0.062 (0.22, 1.32)

Tabela 5.2: Comparacao entre os modelos M1, M2 e M3 na predicao de 500 observacoes sele-

cionadas dentre as localizacoes e tempos. SDA e a soma dos desvios absolutos entre valores observados

e preditos. PROP e a proporcao de vezes que os intervalos preditivos de 95% contem o verdadeiro

valor e TMI e o tamanho medio dos 500 intervalos.

Model SDA PROP TMI

M1 290.5 94.2 2.9

M2 371.7 94.2 3.3

M3 389.3 93.8 3.6


σ23 = 1.42, σ2

4 = 1.04 , σ25 = 1.03, φ1 = 0.012, φ2 = 0.044, φ3 = 0.009, φ4 = 0.011 e

φ5 = 0.012. Os resultados para as componentes 1 e 2 e resultados para as componentes

4 e 5 foram parecidos. Assim, um modelo M22 revisado contendo apenas 3 componentes

foi considerado com os centros formados por u3 e as localizacoes medias entre os pares

de centros vizinhos. Este modelo revisado teve uma performance preditiva pior que

M2 e nao foi considerado.

A mensagem geral destes resultados e de que os dados parecem exibir uma de-

pendencia espacial complexa. Considerar estacionaridade local e superior a estacionari-

dade global mas a descricao nao-estacionaria que e proposta e superior a uma descricao

localmente estacionaria. Assim, somente os resultados para M1 sao descritos abaixo.

5.4.3 Estimacao do termo espacial

As densidades marginais a posteriori e medias a posteriori dos hiperparametros do

termo espacial sao mostrados na Figura 5.6. O problema de convergencia descrito

na estimacao do parametro de suavizacao h1 foi observado. O algoritmo foi entao

reiniciado sem esta componente e os problemas nao mais estiveram presentes. Deste

modo, o modelo com somente 4 componentes e considerado para M1.

Mais uma vez, nenhuma tentativa de interpretacao para as componentes do modelo

e feita aqui e nao ha sentido em se fazer comparacoes entre os parametros estimados

a partir dos modelos M1 e M2, apesar da sua aparente similaridade. A Figura 5.7

apresenta os nıveis de contorno das funcoes componentes estimadas γh(.), para h =

2, 3, 4, 5, geradas pelo modelo M1. Elas indicam correlacoes afetadas para pontos em

regioes bem alem de suas imediacoes, o que reforca a necessidade de um modelo que

va alem da estacionaridade local.

5.4.4 Estimacao da estrutura media

A Tabela 5.3 mostra a estimacao dos coeficientes fixados da estrutura media. Os efeitos

sazonais e da elevacao sao significantes como esperado. Este ultimo e compatıvel com

o decrescimo de 0.65oC usualmente citado na literatura (Ahrens, 2000).


10 15 20 25

10

15

20

25

Predicted values

Ob

se

rve

d v

alu

es

10 15 20 25

10

15

20

25

Predicted values

Ob

se

rve

d v

alu

es

Figura 5.5: Sumario do exercıcio de predicao para 500 observacoes selecionadas aleatoriamente:

as medianas sao representadas pela linha cheia e os intervalos de predicao de 95% pela linha de pontos.

Acima: M1; abaixo: M2.

A tendencia no tempo apresenta resultados conflitantes. A tendencia linear global

b parece ser irrelevante, mas um numero substancial de locais indicam um aumento

significante na temperatura e um numero substancial de locais tambem indicam um

significante decrescimo na temperatura. Resultados similares foram obtidos para os


−1 0 1 2 3 4

02

46

810

Poste

rior

σσ22 σσ3

2σσ42σσ5

2

−10 −8 −6 −4 −2 0 2 4

01

23

4

Poste

rior

φφ2 == φφ5φφ3 φφ4

70 80 90 100 110

0.00

0.05

0.10

0.15

0.20

0.25

Poste

rior

h2 h3h4 h5

Figura 5.6: Densidades marginais estimadas a posteriori dos hiperparametros e suas medias a

posteriori (linha vertical). Os parametros σ2m e φm, m = 2, . . . , 5 aparecem na escala logarıtmica para

melhor inspecao visual. As unidades de medida para os parametros de suavizacao e amplitudes e

kilometros.


0 100 200 300 4000

5015

025

0

γγ2

x

y

u2

u3

u4

u5

0.2

0.4

0.6

0 100 200 300 400

050

150

250

γγ3

x

y

u2

u3

u4

u5

0.1

0.2

0.3

0.4

0.5

0 100 200 300 400

050

150

250

γγ4

x

y

u2

u3

u4

u5

0.2

0.2

0.3

0.4

0.5

0.6

0.7

0 100 200 300 4000

5015

025

0

γγ5

x

y

u2

u3

u4

u5

0.2

0.4

0.6

0.8

Figura 5.7: Nıveis de contorno das componentes estimadas γh(.), for h = 2, 3, 4, 5, geradas pelo

modelo M1.

modelo M2 e M3. A Figura 5.8 apresenta um sumario grafico da estimacao destes

coeficientes. Suas medias estimadas variaram desde um aumento de 1.26oC por decada

ate um declınio de 1.68oC por decada. Estes resultados confirmam a importancia

tendencias no tempo especıficas para cada local.

As estimativas pontuais para φβ e σ2β foram respectivamente 1.01 e 2.5 × 10−5. A

Figura 5.9 apresenta os valores interpolados dos coeficientes de tendencia no tempo

para todo o estado do Rio de Janeiro. Embora a amplitude tenha sido pequena,

duas regioes distintas emergem na figura. Os resultados mostram um significativo

aumento na temperatura em torno do centro u2, correspondente a area metropolitana

da cidade do Rio de Janeiro. Esta e a area do estado onde os valores de temperaturas

mais elevadas sao observadas. Isto parece especialmente relevante na parte oeste da

baıa da Guanabara onde a cidade do Rio de Janeiro se encontra, com urbanizacao e

industrializacao intensas. Os resultados mostram tambem um declınio na regiao em

torno do centro u3. Esta e uma regiao caracterizada por uma densa cobertura de


Tabela 5.3: Sumario dos coeficientes de regressao estimados.

coeficiente media I.C. 95%

intercepto 20.07 (19.88, 20.27)

cos 1.43 (1.25, 1.58)

sin 1.80 (1.69, 1.91)

altitude -0.0062 (-0.0063, -0.0061)

tendencia global -0.0010 (-0.0033, 0.0010)

0 5 10 15 20 25 30 35

−0.

015

−0.

010

−0.

005

0.00

00.

005

0.01

0

ββi

Est

imat

ed v

alue

s

_ _ _

_

_

_

_

_

_ _

_ _

_

_

_

_ _

_

_

_

_

_

_

_

_

_

_

_

_

_

_

_

__

_

_

_

Figura 5.8: Estimativas pontuais e intervalos para as tendencias no tempo especıfica para os locais

βi,i = 1, . . . , 37 para o modelo M1.

floresta e cadeias de montanhas. Esta e a area do estado onde valores de temperatura

mais baixas sao observadas. O mapa tambem mostra um aumento de temperatura

ao redor das cidades de Cabo Frio e Campos. Ambas cidades tiveram um substancial

aumento na urbanizacao nas ultimas decadas devido ao turismo e atividade petrolıfera,

respectivamente. Pela figura 5.9 parece haver um padrao de que locais de temperaturas

mais elevadas estao se tornando mais quentes e que locais de temperatura menos elevada


estao se tornando mais frios, mas isso precisa ser confirmado por estudos posteriores.

0 100 200 300 400

050

100

150

200

250

x

y

Figura 5.9: Valores interpolados do aumento/declınio na temperatura para o estado do Rio de

Janeiro.

5.5 Conclusao

Esta aplicacao propoem uma nova abordagem para predicao e analise de mudancas

climaticas no estado do Rio de Janeiro. A implementacao utiliza um modelo nao-

estacionario baseado em uma convolucao finita de processos estacionarios latentes. As

vantagens de nossa abordagem sao apresentadas e constrastadas favoravelmente contra

alternativas que assumem estacionaridade global e nao-estacionaridade baseada em

estacionaridade local. Os modelos sao bem estimados e o numero de componentes

podem ser identificados do procedimento de estimacao a partir dos dados.

A analise de temperatura mınima mostrou a necessidade em se considerar modelos

espaco-temporais mais elaborados para regioes de grande diversidade geografica. Os


dados revelaram a ausencia de uma tendencia global no tempo para a temperatura indi-

cada, mas um crescimento disperso da temperatura em areas largamente urbanizadas,

areas de temperaturas elevadas se aquecendo e areas de florestas, de temperaturas mais

baixas se tornando mais frias.

Capıtulo 6

Consideracoes finais

Os trabalhos desenvolvidos ate o presente momento indicam novas possibilidades de

procedimentos para analise de dados espacialmente correlacionados de forma nao-

estacionaria. O modelo espacial proposto no Capıtulo 2 e modelos espaco-temporais

propostos nos Capıtulos 3 e 4 foram testados e analisados para exemplos simulados.

Os exemplos simulados do capıtulo 2 mostraram basicamente que: a) o mecanismo de

amostragem adotado e capaz de identificar o numero correto de misturas, bem como

quais devem ser as componentes a serem consideradas no modelo. Este esquema e

robusto no sentido de recuperar os parametros utilizados no processo de simulacao;

b) quando um processo nao-estacionario subjacente e de fato dado pela mistura de

processos estacionarios latentes o modelo proposto e superior no que diz respeito a ca-

pacidade preditiva quando comparado ao modelo de misturas de processos localmente

estacionario (Fuentes e Smith (2001); Fuentes (2002); Banerjee et al. (2004b)); e c) nao

se perde em nada ao considerar estruturas mais gerais, como as descritas no capıtulo

2 mesmo que a estrutura de dependencia espacial dos dados seja mais simples.

No capıtulo 5 dados de temperatura mınima mensal do estado do Rio de Janeiro

foram analisados atraves de um modelo hierarquico definido como uma media de-

pendente de covariaveis explicativas mais a realizacao de um processo espacial nao-

estacionario como descrito no Capıtulo 2. Alem disso, um conjunto de coeficientes

espacialmente estruturados foram considerados as tendencias lineares em cada uma

das localizacoes. Nossa abordagem se mostrou util ao estabelecer o numero de mis-

100


turas a serem consideradas, a partir dos dados. Nosso modelo tambem apresentou

melhor capacidade preditiva quando comparado ao modelo de misturas de processos

localmente estacionarios e ao modelo que considera um unico processo estacionario

para toda a regiao.

O exemplo simulado apresentado do capıtulo 4, mostrou que o esquema de inferencia

e de amostragem adotados foi robusto no sentido de recuperar os parametros utilizados

no processo de simulacao. Deste modo, o proximo passo e realizar a aplicacao com

dados reais de estrutura espaco-temporal complexa que possa ser modelado de acordo

com a abordagem discutida no Capıtulo 4. Uma opcao seriam considerar os dados

de temperatura em escala maior, por exemplo, a regiao sudeste ou o Brasil. Embora

tenhamos desenvolvido o procedimento de inferencia e amostragem na situacao onde

a matriz de covariancia evolui em intervalos de tempos maiores do que evoluem os

dados, outras possıveis abordagens deverao tambem ser testadas. Por exemplo, a de se

considerar as variancias dos processos estacionarios mudando no tempo, conjuntamente

ou nao com os parametros de suavizacao; ou ainda, considerar que estes parametros

mudam no tempo, condicionadas a um conjunto de possıveis covariaveis associadas.

Outras extensoes, podem ser consideradas, no sentido de imputar efeito sazonal

a evolucao da estrutura espacial nao-estacionaria. Seguindo nesta linha, poderia se

assumir uma matriz de covariancia para cada uma das estacoes do ano (ou meses do

ano). Neste caso, cada uma dessas matrizes levariam em conta todas as informacoes

relativas aquela estacao (ou mes) ao mesmo tempo mas estariam condicionalmente as-

sociadas ao conjunto de informacoes da estacao (ou mes) anterior. Isto e naturalmente

obtido no contexto de modelos dinamicos. Para todo m, m = 1, . . . ,M , seja o vetor

sazonal de fatores no tempo t, dado por, hm,t = (hm,t, hm,t−1, . . . , hm,t−p+1)′ onde p

e o comprimento do ciclo sazonal. A equacao dinamica para o vetor e definida no

logaritmo,

log(ht) = α log(ht−1) + ξt. (6.1)


onde α =

0 1 0 . . . 0

0 0 1 . . . 0...

......

. . ....

0 0 0 . . . 1

1 0 0 . . . 0

e uma matriz permutacao e ξt ∼ N(0,W ).

Outra abordagem para processos espaciais nao-estacionarios sao os modelos de de-

formacao espacial, que tiveram como ponto de partida, o artigo de Sampson e Guttorp

(1992). Partindo do espaco das observacoes originais D, eles consideram uma funcao

de correlacao do processo espacial referente a um espaco latente L, onde se aplicam

isotropia e estacionaridade. Seguindo essa ideia, Schmidt e O’Hagan (2003) propoem

uma abordagem Bayesiana na qual o mapeamento entre os espacos D e L seja repre-

sentado por uma funcao d(.) que segue um processo Gaussiano. Recentemente, Castro

(2009) considera extensoes desse modelo, inserindo-o no contexto de modelos dinamicos

para a media do processo e mantendo a estrutura espacial nao-estacionaria definida

pelo modelo de deformacao constante ao longo de todos os tempos. Em outro modelo

Castro (2009) considera que esta estrutura espacial possa tambem mudar ao longo do

tempo atraves da inclusao de uma equacao dinamica para a funcao deformacao d(.).

Embora modelos baseados na deformacao espacial sejam um dos metodos mais estuda-

dos para processos nao-estacionarios, assim como os modelos baseados em convolucao,

nao existe na literatura trabalhos considerando estudos simulados, ou aplicacoes com

dados reais, com o objetivo de comparar os modelos de deformacao espacial com os

modelos baseados na convolucao de processos. Nos esperamos, em breve, comparar as

capacidades preditivas destes modelos atraves de aplicacoes com dados reais e dados

simulados.

Neste trabalho e considerado que os dados tenham distribuicao gaussiana. Entre-

tanto, muitas vezes nao e possıvel fazer esta suposicao para dados gerais, dado sua

distribuicao assimetrica. Neste caso, utilizar transformacoes como por exemplo, o log-

aritmo dos dados, levam a dificuldades na interpretacao dos parametros do modelo.

Deste modo, a extensao dos modelos descritos para dados da famılia exponencial seria

bastante interessante. Neste caso, poderiam ser ajustados tambem, dados de contagem

que rotineiramente aparecem na area da saude. Outra saıda seria considerar uma


abordagem nao-parametrica mais irrestrita sobre a distribuicao dos dados atraves de

os modelos nao-parametricos de misturas baseados em processos de Dirichlet (Gelfand

et al. (2005b); Kottas et al. (2008)).

Referencias Bibliograficas

Agarwal, D. e Gelfand, A. (2005), “Slice sampling for simulation based fitting of spatial

data models,” Statistics and Computing, 15, 61–69.

Ahrens, A. D. (2000), Essentials of Meteorology: an invitation to the atmosphere,

Brooks Cole.

Banerjee, S., Carlin, B., e Gelfand, A. (2004a), Hierarchical Modeling and Analysis for

Spatial Data, Chapman & Hall.

Banerjee, S., Gelfand, A., Knight, J., e Sirmans, C. (2004b), “Spatial modelling of house

prices using normalized distance-weighted sums of stationary processes,” Journal of

Business and Economic Statistics, 22, 206–213.

Barber, J. (1999), “Modeling and prediction of non-stationary environmental pro-

cesses,” Tese de doutorado, North Carolina State University.

Berger, J. O., Oliveira, V., e Sanso, B. (2001), “Objective bayesian analysis of spatially

correlated data,” Journal of the American Statistical Association, 96, 1361–1374.

Calder, C. A. (2003), “Exploring latent structure in spatial temporal processes using

process convolutions,” Tese de doutorado, ISDS, Duke University.

— (2007), “Dynamic factor process convolution models for multivariate sace-time data

with application to air quality assessment,” Environmental and Ecological Statistics,

14, 229–247.

Carter, C. K. e Kohn, R. (1994), “On Gibbs sampling for state space models,”

Biometrika, 81, 541–553.

104


Castro, F. (2009), Modelos dinamicos para deformacoes espaciais, Exame de Quali-

ficacao de Doutorado, Instituto de Matematica, DME-UFRJ, Rio de Janeiro, Brasil.

Cressie, N. (1993), Statistical for spatial Data, New York: Wiley.

ECSN (1995), Climate of Europe: Recent Variation, Present State and Future

Prospects, KNMI, De Bilt.

Epanechnikov, V. (1969), “Nonparametric estimates of a multivariate probability den-

sity,” Theory of Probability and its Applications, 14, 153–158.

Fruhwirth-Schnater, S. (1994), “Data augmentation and dynamic linear models,” Jour-

nal of Time Series Analysis, 15, 183–202.

Fuentes, M. (2002), “Spectral methods for nonstationary spatial processes,”

Biometrika, 89, 197–210.

Fuentes, M., Chaudhuri, A., e Holland, M. D. (2007), “Bayesian entropy for spatial

sampling design of environmental data,” Environmental and Ecological Statistics, 14,

323–340.

Fuentes, M., Chen, L., Davis, J. M., e Lackmann, G. M. (2005), “Modeling and pre-

dicting complex space-time structures and patterns of coastal wind fields,” Environ-

metrics, 16, 449–464.

Fuentes, M., Henry, J., e Reich, B. J. (2009), “Nonparametric spatial models for ex-

tremes: application to extreme temperature data,” Tech. rep., Statistics Department,

North Carolina State University. Institute of Statistics Mimeo Series 2621.

Fuentes, M. e Reich, B. J. (2009), “Multivariate spatial nonparametric modelling via

kernel processes mixing,” Tech. rep., Statistics Department, North Carolina State

University. Instute of Statistics Mimeo Series 2622.

Fuentes, M. e Smith, R. (2001), A New class of nonstationary spatial models, Relatorio

tecnico, North Carolina State University.


Gelfand, A., Banerjee, S., e Gamerman, D. (2005a), “Spatial process modelling for

univariate and multivariate dynamic spatial data,” Environmetrics, 16, 465–479.

Gelfand, A., Kim, H., Sirmans, C., e Banerjee, S. (2003), “Spatial modeling with spa-

tially varying coefficient processes,” Journal of the American Statistical Association,

98, 387–396.

Gelfand, A., Kottas, A., e MacEachern, S. (2005b), “Bayesian nonparametric spatial

modeling with Dirichlet process mixing,” Journal of the American Statistical Asso-

ciation, 100, 1021–1035.

Gelfand, A., Schmidt, A., Banerjee, S., e Sirmans, C. (2004), “Nonstationary mul-

tivariate process modeling through spatially varying coregionalization,” Test, 13,

263–312.

Gelman, A. (2006), “Prior Distributions for Variance Parameters in Hierarchical Mod-

els,” Bayesian Analysis, 1, 515–533.

Green, P. J. (1985), “Reversible jump Markov chain Monte Carlo computation and

Bayesian model determination,” Biometrika, 82, 711–732.

Hass, T. (1995), “Local prediction of a spatio-temporal process with an application to

wet sulfate deposition,” Journal of the American Statistical Association, 90, 1189–

1199.

Higdon, D. (1998), “A process-convolution approach to modelling temperatures in the

North Atlantic Ocean,” Environmental and Ecological Statistics, 5, 173–190.

— (2002), “Space and space-time modeling using process convolutions,” in Quantitative

Methods for Current Environmental Issues, eds. Anderson, C., Barnett, V., Chatwin,

P. C., e El-Shaarawi, A. H., Springer Verlag, pp. 37–56.

Higdon, D., Swall, J., e Kern, J. (1999), “Non-stationary spatial modelling,”

in Bayesian Statistics 6 - Proceedings of the Sixth Valencia Meeting, eds.

J. M. Bernardo, J. O. Berger, A. P. D. e Smith, A. F. M., Claredon Press - Ox-

ford, pp. 761–768.


Huerta, G., Sanso, B., e Stroud, J. R. (2004), “A spatio-temporal model for Mexico

city ozone levels,” Applied Statistics, 53, 231–248.

IPCC (2007), Climate Change 2007: The Physical Science Basis. Contribution of

Working Group I to the Fourth Assessment Report of the IPCC, Cambridge Uni-

versity Press.

Kern, J. C. (2000), “Bayesian Process-Convolution Approaches to specifying Spatial

Dependence Structure,” Tese de doutorado, Tese de doutorado, Duke University,

Durham, NC 27708.

Kim, H., Mallick, B., e Holmes, C. (2006), “Analyzing spatial data using piecewise

Gaussian processes,” Journal of the American Statistical Association, 100, 653–668.

Kottas, A., Duan, J., e Gelfand, A. (2008), “Modeling disease incidence data with

spatial and spatio temporal Dirichlet process mixtures,” Biometrical Journal, 50,

29–42.

Lopes, H., Salazar, E., e Gamerman, D. (2008), “Spatial dynamic factor analysis,”

Bayesian Analysis, 3(4), 759–792.

Nadaraya, E. A. (1964), “On Estimating Regression,” Theory of Probability and its

Applications, 9, 141–142.

Neal, R. M. (2003), “Slice sampling,” The Annals of Statistics, 3, 705–767.

Nychka, D. e Saltzman, N. (1998), “Design of Air Quality Networks,” in Case studies

in environmental statistics, ed. in Statistics, L. N., Springer verlag, pp. 51–76.

Oke, T., Klysik, K., e Bernhofer, C. (2006), “Editorial: Progress in urban climate,”

Theoretical and Applied Climatology, 84, 1–2.

Paciorek, C. e Schervish, M. (2004), “Nonstationary covariance functions for Gaussian

process regression,” Advances in Neural Information Processing Systems, 16, 273–

280.


— (2006), “Spatial modelling using a new class of nonstationary covariance functions,”

Environmetrics, 17, 483–506.

Paez, M., Gamerman, D., Landim, D., e Salazar, E. (2008), “Spatially-varying dynamic

coefficient models,” Journal of Statistical Planning and Inference, 138, 1038–1058.

Paez, M., Gamerman, D., e Oliveira, V. (2005a), “Interpolation performance of a

spatio-temporal model with spatially varying coefficients: application to PM10 con-

centrations in Rio de Janeiro,” Environmental and Ecological Statistics, 12, 169–193.

Paez, M. S., Gamerman, D., e de Oliveira, V. (2005b), “Interpolation performance

of a spatio-temporal model with spatially varying coefficients: application to PM10

concentrations in Rio de Janeiro,” Environmental and Ecological Statistics, 12, 169–

193.

Sahsamanoglou, H. e Makrogiannis, T. (1991), “Temperature trends over the Mediter-

ranean region,” Theoretical and Applied Climatology, 45, 183–192.

Sampson, P. e Guttorp, P. (1992), “Nonparametric estimation of nonstationary spatial

covariance structure,” Journal of the American Statistical Association, 87, 108–119.

Sanso, B., Schmidt, A., e Nobre, A. (2008), “Interpolation performance of a spatio-

temporal model with spatially varying coefficients: application to PM10 concentra-

tions in Rio de Janeiro,” Canadian Journal of Statistics, 36, 239–258.

Schmidt, A. e Gelfand, A. (2003), “A Bayesian Coregionalization Approach for Multi-

variate Pollutant Data,” Journal of Geophysical Research-Atmospheres, 108, 6783.

Schmidt, A. e O’Hagan, A. (2003), “Bayesian inference for nonstationary spatial co-

variance structures via spatial deformations.” Journal of the Royal Statistical Society

Series B, 65, 743–758.

Schmidt, A. e Sanso, B. (2006), Modelagem Bayesiana da Estrutura de Covariancia

de Processos Espaciais e Espaco-Temporais, Minicurso 17o Simposio Nacional de

Probabilidade e Estatıstica, ABE. Caxambu, MG, Brasil.


Smith, R. (2001), Environmental Statistics, Relatorio Tecnico, University of North

Carolina.

Stroud, J. R., Muller, P., e Sanso, B. (2001), “Dynamic model for spatio-temporal

data,” Journal of the Royal Statistical Society, B,, 63, 673–689.

Swall, J. L. (1999), “Non-Stationary spatial modeling using a process convolution

approach,” Tese de doutorado, Tese de doutorado, Duke University, Durham, NC

27708.

Watson, G. S. (1964), “Smooth regression analysis,” Shankya series A, 26, 359–372.

West, M. e Harrison, J. (1997), Bayesian Forecasting and Dynamic Models, New York:

Springer - Second Edition.

Documents

Universidade Federal do Rio de Janeiro Instituto de Matem ... · Universidade Federal do Rio de Janeiro Instituto de Matem atica Modelos de convolu˘c~ao para dados espa˘co-temporais