Estimação Paramétrica De Modelos Auto-Regressivos via Estatística Bayesiana e Simulação De Monte Carlo

ESTIMAÇÃO PARAMÉTRICA DE MODELOS AUTO-REGRESSIVOS VIA ESTATÍSTICA BAYESIANA E SIMULAÇÃO DE MONTE CARLO

GUILHERME DE A. BARRETO

MARINHO G. DE ANDRADE

Departamento de Engenharia Elétrica Universidade de São Paulo Av. Dr. Carlos Botelho, 1465 CEP. 13650-970, São Carlos, SP

Inst. de Ciências Matemáticas e de Computação Universidade de São Paulo Av. Dr. Carlos Botelho, 1465 CEP. 13650-970, São Carlos, SP

Resumo - Neste trabalho é proposta uma abordagem bayesiana para a estimação paramétrica de modelos estocásticos autoregressivos de ordem p, AR(p). Serão também apresentados procedimentos para selecionar a ordem p do modelo, fazer previsões e avaliar a robustez dos estimadores por meio de técnicas de simulação de Monte Carlo em Cadeias de Markov (MCMC). A abordagem proposta é comparada com o método clássico de Box-Jenkins (máxima verossimilhança) em uma série histórica de vazões médias mensais do reservatório de Furnas. Conclui-se que o uso de estatística bayesiana e técnicas de simulação MCMC tornam o processo de estimação paramétrica mais flexível e poderoso quando comparado com técnicas clássicas. Abstract - In this work we propose a Bayesian approach for the parameter estimation problem of stochastic autoregressive models of order p, AR(p). Procedures for order determination, forecasting and robustness evaluation through Monte Carlo Markov Chain (MCMC) simulation techniques are also presented. The proposed approach is compared with the classical one by Box-Jenkins (maximum likelihood estimation) on a streamflow time series from Furnas reservoir. We conclude that the use of Bayesian statistics and MCMC simulation implies in more flexible and powerful results than those obtained from the classical approach. Keywords - System Identification, autoregressive models, Bayesian statistics, Markov chain, streamflow forecasting, robustness.

1 Introdução Seja um processo estocástico {Zt , t ∈ N+}. Um modelo linear auto-regressivo de ordem p, AR(p), é aquele cujo o valor corrente do processo Zt é expresso como uma combinação linear de seus p valores passados Zt-1, Zt-2, ...., Zt-p e de um ruído branco at. Considera-se que o ruído {at , t ∈ N+} é um processo i.i.d. N(0,σ2). Um modelo AR(p) pode, então, ser escrito da seguinte forma: Zt = φ1Zt-1 + φ2Zt-2 + .... + φpZt-p + at

(1)

Neste trabalho, utilizam-se técnicas de estatística bayesiana e simulação de Monte Carlo em cadeias de Markov (MCMC), para realizar as seguintes etapas: (i) determinação da ordem p que melhor representa a série temporal observada, (ii) estimação dos parâmetros φ1, φ2 ,...., φp e τ=1/σ2, (iii) previsão de valores e (iv) análise de robustez do estimador bayesiano. Estas técnicas são comparadas com a abordagem “clássica” para este tipo de problema usando o estimador de máxima verossimilhança (Box et al., 1994). As técnicas desenvolvidas são aplicadas em um problema de previsão de vazões mensais médias do reservatório de Furnas. Este é um problema de fundamental importância para o planejamento da operação do sistema energético, em que visa-se uma tomada de decisão que resulte na operação ótima do sistema (Pereira et al., 1990; Soares, 1990). Na Seção 2, discute-se diferenças entre a abordagem clássica e bayesiana para o problema de estimação paramétrica. Na Seção 3, conceitos relacionados a simulação estocástica MCMC. Na Seção 4 e 5, são apresentados técnicas para determinação da ordem do modelo e previsão de

valores. Na Seção 6, são apresentados os resultados da simulação com a série de Furnas. O artigo é concluído na Seção 7. 2 Estimação Paramétrica do Modelo AR(p) via Estatística Clássica e Bayesiana A abordagem padrão em estimação paramétrica de um processo estocástico AR(p) consiste na especificação da função de verossimilhança L(Z|θ), onde Z = (Z1, Z2, ..., ZN) é a série temporal e θ = (φ1, φ2 ,...., φp, τ) é o vetor de parâmetros (desconhecidos). Pode-se mostrar que a função de verossimilhança, dado o modelo AR(p) em (1) e conhecendo-se as pprimeiras observações Zp, é dada por: N−p τ 2

τ [(φ- φˆ )X'X(φ- φˆ ) 2 + (Z- Zˆ )'(Z- Zˆ )]}

L(Z|φ,τ) ∝

exp{ −

(2)

onde a estimativa de Máxima Verossimilhança φˆ é dada por (Box et al., 1994; Davis e Vinter, 1985): φˆ = (X ' X )−1 (X ' Z ) σˆ 2 = τ −1 =

(

)(

1 Z − Xφˆ ' Y − Xφˆ N−p

)

Z1   Zp  Z p +1      Z2   Z p +1  Z p+2  onde Z =  .  e X=       Z   Z N −1  Z − N p N     O método de máxima verossimilhança considera o parâmetro como uma constante, enquanto a abordagem bayesiana entende o parâmetro θ como








uma variável aleatória que obedece a uma determinada distribuição de probabilidade. O processo de inferência bayesiana requer também a especificação de uma densidade de probabilidade a priori p(θ) que reflete o conhecimento prévio a respeito da distribuição do parâmetro θ em questão. A verossimilhança L(Z|θ) e a densidade a priori p(θ) são combinadas usando o teorema da Probabilidade de Bayes (Gelman et al., 1995) para fornecer a densidade a posteriori f(θ|Z): f (θ | Z ) =

p(θ)L(Z | θ) p(Z )

(3)

onde p(Z) é um fator de normalização. Em geral, escreve-se apenas, f(θ|Z) ∝ p(θ)L(Z|θ)

(4)

A densidade a posteriori f(θ|Z) diz como a variável aleatória θ está distribuída após os dados terem sido observados (daí o adjetivo a posteriori), e a partir dela várias informações podem ser obtidas. Dois pontos devem ser enfatizados aqui. Primeiro, o método bayesiano é mais flexível do que a técnica de máxima verossimilhança, pois permite testar várias possibilidades para a densidade a priori p(θ) e escolher aquela que resultar em uma densidade a posteriori f(θ|Z) que melhor se adequa ao problema. Segundo, é possível que f(θ|Z) não seja uma densidade de probabilidade padrão. Assim, valores médios e variâncias deverão ser obtidos por integração, o que pode se tornar uma tarefa difícil. A solução encontrada é estimar a distribuição a posteriori f(θ|Z) a partir de uma amostra gerada por simulação de Monte Carlo em Cadeias de Markov. 3 Simulação de Monte Carlo em Cadeias de Markov (MCMC)

A forma da função de verossimilhança dada na Eq. (2) sugere uma densidade a priori Normal para φ e uma Gama para τ, ou seja, p(φ, τ)=p1(φ|τ)p2(τ)

(6)

tal que p 2 p1(φ|τ) ∝ τ

exp{ −

τ [(φ-µ)’P-1(φ-µ)]} ~ N(µ, P-1) 2

p 2 (τ ) ∝ τ α −1 exp{− βτ} ~ Gama(α, β)

onde o vetor de médias µ e a matriz de covariância P1 em p1(φ|τ), α e β em p2(φ|τ) são especificados pelo usuário. Assim, a posteriori conjunta, calculada de acordo com a Eq. (4), é dada por: ( N + 2α) −1 f(φ,τ|Z) ∝ τ 2 exp

{ −τD −

τ (φ-φb)’V(φ-φb)} 2

1 [(Z ' Z + µ' Pµ) − 2 ( X ' Z + Pµ)' ( X ' X + P) −1( X ' Z + Pµ)] . As densidades condicionais são facilmente obtidas devido à conjugação: onde V = (X'X + P) e D = β +

f(φ|τ, Z) ~ N(φb, (τV)-1) (7) N + 2α 1 f(τ|φ, Z) ~ Gama( , D + (φ − φ b )'V (φ − φ b ) 2 2 As densidades em (7) estão na forma apropriada para a introdução do Algoritmo Gibbs-Sampling (Geman e Geman, 1984), descrito a seguir: Passo 1. Atribua valores φ(0) = φ1(0) , , φ (p0) e τ(0).

(



)

iniciais aos parâmetros

Obtenha um novo valor φ(j+1) a partir da densidades condicional: φ(j+1) ~ f(φ|τ(j), Z) Passo 2.

Obtenha um novo valor τ(j+1) a partir da densidade condicional: τ(j+1) ~ f(τ|φ(j+1), Z) Passo 3.

O objetivo da simulação MCMC é gerar uma amostra {θ(j) = (φ(j), τ(j)), j = 1, ..., m} a partir das densidades condicionais f(φ|τ, Z), φ = (φ1, φ2 ,...., φp), e f(τ|φ, Z) obtidas a partir da densidade conjunta a posteriori f(φ,τ|Z). Essa amostra é utilizada para obter inferências a cerca dos parâmetros φ e τ, como por exemplo, as estimativas pontuais φˆ e τˆ , dadas por: 1 φˆ = m

m

∑ φ ( j) j =1

e τˆ =

1 m

Verifique se houve convergência. Caso contrário, faça j=j+1 e vá para o passo 2. Passo 4.

Após um número suficientemente grande de iterações, as observações geradas (φ(j), τ(j)) convergirão para uma amostra aleatória da densidade conjunta f(φ,τ|Z) a posteriori.

m

∑ τ( j)

(5)

j =1

3.1. Densidade a Priori Normal-Gama Uma primeira escolha para densidade a priori em geral leva em consideração a possibilidade de se obter uma densidade a posteriori, ou que pelo menos, produza densidades condicionais a posteriori f(φ|τ, Z) e f(τ|φ, Z) que gerem facilmente um número aleatório.

3.2. Densidade a Priori t-Gama No caso mais geral, a densidade a posteriori não produz densidades condicionais padrões. Isto ocorre se for escolhida para φ uma densidade a priori tStudent, p-dimensional, com k graus de liberdade, vetor de localização µ e matriz de precisão P, e uma Gama para τ, ou seja,

 (φ − µ )' P (φ − µ ) p1 (φ) ∝ 1 +  k  

−

(k + p) 2

(8)

p 2 (τ ) ∝ τ α −1 exp{− βτ} Para este caso, a posteriori é dada por: f(φ,τ|Z) ∝ tal que

N − p +2α −1 exp{-τ[β+ τ 2

( )

4 Seleção da Ordem p via Simulação MCMC

B (φ) ]}p1(φ) 2

( )(

)(

B(φ) = φ − φˆ ' X ' X φ − φˆ + Z − Zˆ ' Z − Zˆ  (φ − µ )' P (φ − µ ) p1(φ) ∝ 1 +  κ  

−

)

( )

(κ+ p ) 2

τ [ (φ − φˆ )' ( X' X)(φ − φˆ ) ]}Π1(φ) 2 = Ψ(φ).p1(φ) 1 N − p + 2α ,β+ B(φ) ) (9) f(τ|φ, Z) ~ Gama( 2 2 f(φ|τ, Z) ∝ exp{-

A densidade condicional f(φ|τ, Z) não pode ser amostrada diretamente. Assim, utiliza-se o algoritmo Metropolis-Hastings (Brooks, 1998) cujo os passos são dados a seguir.

(

)



Dentre os critérios mais utilizados para seleção de modelos estão os critérios AIC (Akaike, 1974) e BIC (Schwarz, 1978) dados por: AIC = ln σˆ 2a +

As densidades condicionais são dadas por:

Passo 1. Atribua valores (0) φ(0) = φ1 , , φ (p0) e τ(0).

A convergência das cadeia será verificada pelo método de Gelman-Rubin (Gelman et al., 1995). (4) Após a convergência das cadeias, as mesmas são concatenadas para formar uma única longa cadeia de tamanho M.n.

iniciais aos parâmetros

Obtenha novos valores a partir da função de transição na Eq. (8): φ(j+1) ~ p1(φ(j)). Passo 2.

Calcule a probabilidade de aceitação do novo valor: α(φ

 Ψ(φ( j +1) )  , φ ) = min1, , î Ψ (φ( j ) )  (j)

Mp

ck

Gere u de Uniforme[0,1] e faça φ u > α(φ(j+1),φ(j)). Passo 4.

Passo 5.

(j+1)

= φ se (j)

Gerar τ(j+1) ~ f(τ|φ, Z) a partir da Eq. (9).

Verifique se houve convergência. Caso contrário, faça j=j+1 e repita passos 2-6. Passo 6.

3.3. Procedimentos de Simulação Os seguintes procedimentos são efetuados com o intuito de gerar uma amostra {θ(j) = (φ(j), τ(j)), j = 1, ..., n} por simulação MCMC. (1) Um número M de cadeias independentes são simuladas ao mesmo tempo por n iterações e para diferentes valores iniciais do parâmetro. Ou seja, tem-se uma matriz { θ ij }, i = 1, ..., M; j = 1, ..., n, cujas linhas indicam a cadeia e as colunas a iteração correspondente dentro daquela cadeia. (2) A metade inicial da cadeia é descartada para evitar estimativas viciadas (biased estimates) devido a influência dos valores iniciais. (3)

( )

BIC = ln σˆ 2a +

M ln(N ) N

 τ = f(ZN+k|ZN+k-1, φ, τ) ∝ ∫φ∫τ τ1/2exp − (Z N + k − î 2

∑i =1 φ i Z N + k −i )  2

p

f(φ, τ|Z)dφdτ

(10)

onde f(φ, τ|Z) é a densidade a posteriori para os parâmetros φ e τ, e Mp é o modelo AR(p) sendo avaliado. A Eq. (10) pode também ser calculada usando as observações φ ( j ) e τ ( j ) obtidas por simulação MCMC. Neste caso, a densidade preditiva é estimada por cˆ k = fˆ ( Z N + k | Z N + k −1 ) : Mp

k = cˆM p

Ψ(φ(j) ) ≠ 0

e

onde M = p+1 é o número de parâmetros e N é o tamanho da série em questão. Bastante comuns também são aqueles baseados na densidade preditiva a posteriori (Gelman et al., 1995). A densidade preditiva ordenada para a observação ZN+k dada a série ZN+k-1 = (Z1, ..., ZN, ..., ZN+k-1)’, φ e τ, é dada por (Andrade e Barreto, 1999):

Passo 3.

(j+1)

2M N

2  ( j)  p   1 m ( j) 1/ 2  τ ( j) exp− Z N + k − ∑ φi Z N + k −i   ∑τ m j =1  2  i =1   î

( )

Calcula-se cˆ M p =

∏k =1 cˆ Mk p K

para cada modelo Mp

e escolhe-se o que produzir o maior valor para cˆ M p . Outra possibilidade é determinar o modelo via fator de Bayes, definido como Bij = p(Z|Mi )/p(Z|Mj), onde p(Z|Mp) é o fator de normalização (ver Eq. (3)) definida como: p(Z| Mp) = ∫φ∫τ L(Z|φ, τ, Mp)p(φ, τ|Mp)dφdτ onde p(φ, τ|Mp) é a densidade a priori e L(Z|φ, τ, Mp) é a função de verossimilhança dada na Eq. (2). Note que p(Z| Mp) também pode ser estimada usando a amostra gerada por simulação MCMC como:

(

)

pˆ Z | M p =

1 m

∑ L(Z | φ ( j ) , τ ( j ) , M p ) m

j =1

Prefere-se o modelo Mj ao modelo Mi se Bˆ ij < 1 .

(11)

5

Previsão via Simulação MCMC

Seja Zf = (Z(N+1), Z(N+2), ..., Z(N+k))’ os k valores da série a serem previstos e Zp = (Z(N-(p-1), Z(N-(p-2),, ..., ZN)’ as p últimas observações da mesma. Pode-se mostrar (Andrade e Barreto, 1999) que Zf é um vetor de variáveis aleatórias que obedece a uma distribuição N[E(Zf|φ, τ, Zp), Var(Zf|φ, τ, Zp)], tal que: E(Zf|φ, τ, Zp) = -A-1(φ)B(φ)Zp

(12)

Var(Zf|φ, τ, Zp) = diag{τ -1(A-1(φ))( A-1(φ))'} (13) onde as matrizes A(φ) e B(φ) são dadas por: 0  1  1  − φ1 − φ1 A(φ)=  − φ2    − φ k −1 − φk − 2 

0 0 1



 −φ p   0 B(φ)=  0    0

−φ p −1 −φp 0









− φ k −3

−φ p − 2 0 −φp





Para a série Zt ajusta-se um modelo AR(p). A Fig. 1 traz a série Furnas antes e depois da transformação. Nota-se que a série original apresenta sazonalidades referentes a períodos de seca e chuvas intensas ao longo dos anos, cuja transformação não linear tem como objetivo eliminar.

−φ1   − φ2  − φ 3    − φ k  k× p 





− φ k −3

0

log y t − µ m ~ N(0,1) σm







log yt ~ N(µm, σm) e Z t =

0  0 0    1  k ×k 



seja, yt ~ log-N(µm, σm) onde m = 1, 2, ..., 12, representa o número de meses, e µm e σm são as médias e desvios-padrão amostrais para cada mês do ano. Desta forma, tem-se que:





As Eqs. (12) e (13) podem ser calculada usando simulação de MCMC através da seguinte expressão: Ê(Zf|Zp) =

( )( )

1 m ∑ − A −1 φ ( j ) B φ ( j ) Z p m j =1

1 Vˆar (Zf|Zp) = diag m  +

1 m

m

1

j =1 τ

( j)

∑

(14)

( ) ( )

' A −1 φ ( j ) A −1 φ ( j )   

∑ (A −1 (φ ( j ) )B(φ ( j ) )Z p ) m

2

(15)

j =1



m



j =1

( )( )



2

+  1 ∑ A−1 φ( j ) B φ( j ) Z p  m  

Assim, as Eqs. (14) e (15) são, respectivamente, as estimativas de Monte Carlo dos valores previstos e da variância destes valores.

Figura 1. Série de vazões para a usina de Furnas A seleção da ordem do modelo AR(p) segundo os critérios AIC, BIC, densidade preditiva ordenada e Fator de Bayes são mostrados na Tabela 1. Para o Fator de Bayes, as razões B12 = 0,0107/0,5581 = 0,02 < 1 e B32 = 0,1636/0,5581 = 0,29 < 1 indicam que deve ser dada preferência ao modelo AR(2). Tabela 1. Seleção do Modelo

6 Estudo de Caso – Série Furnas Neste trabalho é considerada uma série temporal cujos dados são referentes às vazões médias mensais que chegaram ao reservatório de Furnas no período de Janeiro de 1960 a Dezembro de 1990. Todos os procedimentos descritos até a Seção 5 foram implementados na linguagem do software MATLAB® com o auxílio das toolboxes de Identificação e Estatística. As observações reais das medidas de vazão média são denotadas yt. Assume-se que estas observações seguem uma distribuição log-Normal, ou

AIC BIC DPO (k=20) Fator de Bayes

-1,0350 -1,0243 2,4215 0,0107

-1,0467 -1,0255 6,4706 0,5581

-1,0456 -1,0141 1,8854 0,1636

A Fig. 2 apresenta os histogramas das amostras obtidas para o caso de uma priori t-gama para os parâmetros φi, i=1 e 2, e τ da série de vazões de Furnas. O índice b indica que a estimativa é bayesiana. O caso da priori normal-gama produz resultados semelhantes e não são mostrados. A Tabela 2 traz as estimativas dos parâmetros do modelo AR(2) ajustado a série Furnas. EMV são

as estimativas de máxima verossimilhança. PNG são as estimativas obtidas com a priori normal-gama via Amostrador de Gibbs. PTS são as estimativas obtidas com a priori t-gama via Metropolis-Hastings. Nesta tabela também é mostrada a variância da estimativa. Estes valores foram calculados usando a Eq. (5) depois que o critério de Gelman-Rubin acusou convergência das cadeias. Método EMV * PNG

Tabela 2. Estimativas pontuais e variância φ1 φ2 Var(φ1) Var(φ2) Var( τ) τ 0,6478 0,0027 0,1623 0,0027 2,0285 --0,6435 0,0026 0,1644 0,0028 1,9987 0,0217

** 0,6480 0,0025 0,1619 0,0027 1,9488 0,0317 PTG * Parâmetros priori's normal-gama: µ = (0 0)', P = I2, α = 2, β = 2.

ˆ ** Parâmetros priori's t-gama: µ = φ EMV , P = 75I2, α = 2, β = 2.

φˆ1b

onde yto e y tp são os valores de vazão observados e previstos, respectivamente. A Tabela 3 mostra as previsões num horizonte de 12 passos à frente fornecidas pela abordagem clássica e pela bayesiana em que foi utilizada uma densidade a priori t-gama para a série Furnas. Tabela 3. Previsão 12 passos a frente: Real Prev. clássica Prev. bayesiana 1450 1264,7 1270,2 1512 1269,1 1275,4 1977 1246,4 1253,5 971 896,2 899,9 839 677,5 679,7 639 571,2 572,7 624 487,4 488,3 513 411,4 412,0 354 416,3 417,0 421 506,4 506,9 587 715,0 715,7 1447 1224,4 1225,3 73,10 72,27 MSE 18,34% 18,18% MRE

série Furnas IC (95%) 1223,4 1320,8 1224,8 1322,7 1202,1 1300,3 864,4 930,6 655,1 700,1 552,2 589,1 473,9 499,1 400,5 420,2 402,4 426,8 492,9 515,7 695,2 727,6 1191,8 1243,5 -----------

Os valores na Tabela 3 correspondem aos valores no domínio original de medidas de vazões (e não àquele da série Zt) obtidos a partir da seguinte transformação inversa: E(yt|φ, τ) = exp{µm+

σm [2E(Zt|φ, τ) + σmV(Zt|φ, τ)]} 2

onde E(Zt|φ, τ) e V(Zt|φ, τ) são, respectivamente, o valor esperado e a variância gerados pela previsão clássica e bayesiana calculados na Seção 5. Percebese que o desempenho da técnica bayesiana é ligeiramente melhor que o da técnica clássica. A Fig. 3 mostra os valores das previsões fornecidos pela técnica bayesiana mostrados na Tabela 3, juntamente com os valores reais de uma parte da série Furnas.

φˆ 2b

τˆ b

Figura 2: Histograma gerado por simulação MCMC para os parâmetros φ1, φ2 e τ com priori t-gama.

As próximas simulações avaliam a capacidade preditiva das formulações clássicas de Box-Jenkins e a bayesiana. A previsão de vazões pelo método clássico e bayesiano são avaliados pelo erro médio quadrático (MSE) e erro médio relativo (MRE%), definidos como: 1 k p 2 MSE = ∑ y to − y t k t =1

(

100 k MRE = ∑ k t =1

)

y to − y tp y tp

Figura 3: Previsão 12 passos à frente para a série Furnas. Valores previstos (círculo) e reais (asterisco).

Nos testes a seguir, foram calculados diversos valores pontuais das estimativas bayesianas correspondendo a uma determinada faixa de variação do vetor µ presente em ambas densidades a priori. A matriz de precisão P é mantida em uma valor constante e alto para reduzir a variabilidade de modo a facilitar a visualização. A Fig. 4 mostra o resultado. Nesta figura, φˆ i , i=1 e 2, é a estimativa de máxima verossimilhança e φˆ ib , i=1 e 2, é a bayesiana.

•

•

φˆ 1b − φˆ 1

• •

µ1 − φˆ 1

• •

φˆ 2b − φˆ 2

µ2 − φˆ 2

Figura 4: Análise de Robustez usando priori normalgama (contínuo) e t-gama (tracejado) - série Furnas.

Principalmente na região central a inclinação da curva para a densidade a priori t-gama é menor do que aquela para a curva da densidade a priori normalgama. Isto indica que a densidade a priori t-gama produz estimativas bayesianas com menor sensibilidade a variações do vetor µ. É importante ressaltar que a realização deste tipo de teste só é possível graças a simulação MCMC. A Tabela 4 traz os valores numéricos referentes a Fig. 4. Tabela 4. Valores de

φˆ 1b NG

TG

0,5699 0,5811 0,5932 0,6042 0,6149 0,6288 0,6408 0,6521 0,5780 0,5959 0,6012 0,6076 0,6172 0,6227 0,6306 0,6434

φˆ1b e φˆ2b

Int. Cred. 95% 0,5211 0,5320 0,5416 0,5584 0,5663 0,5763 0,5876 0,6017 0,5100 0,5261 0,5373 0,5317 0,5476 0,5538 0,5588 0,5774

0,6205 0,6301 0,6473 0,6517 0,6618 0,6746 0,6916 0,6978 0,6535 0,6713 0,6680 0,6745 0,6816 0,6901 0,7100 0,7076

para a Furnas.

φˆ 2b 0,0822 0,0943 0,1062 0,1157 0,1278 0,1410 0,1527 0,1653 0,0976 0,1028 0,1149 0,1229 0,1289 0,1379 0,1438 0,1495

Int. Cred. 95% 0,0366 0,0467 0,0571 0,0674 0,0805 0,0948 0,1068 0,1167 0,0213 0,0333 0,0443 0,0517 0,0555 0,0671 0,0767 0,0785

0,1290 0,1430 0,1592 0,1666 0,1838 0,1882 0,2016 0,2107 0,1775 0,1747 0,1803 0,1891 0,1967 0,2146 0,2103 0,2247

7 Conclusão Uma das principais aplicações de simulação de Monte Carlo está na análise bayesiana de séries temporais (Tweedie, 1998). Para este fim, os procedimentos propostos neste trabalho foram os seguintes:

Determinação da função de verossimilhança dos dados, especificação da densidade a priori dos parâmetros, e cálculo da densidade a posteriori via teorema de Bayes. Análise inicial da série buscando observar propriedades particulares da mesma como tendências, sazonalidades, etc. e posterior determinação da FAC e FACP. Determinação da ordem do modelo via métodos clássicos como AIC, BIC, densidade preditiva ordenada e fator de Bayes. Estimação dos parâmetros via métodos clássicos como o de máxima verossimilhança, e através de técnicas bayesianas + simulação MCMC. Previsão de valores para a série sendo analisada por meio técnicas clássicas e bayesianas + simulação MCMC. Análise da robustez dos estimadores bayesianos via simulação MCMC.

Estes dois últimos procedimentos são de particular importância e compreendem os principais resultados teóricos e práticos deste trabalho. Para trabalhos futuros, sugere-se a comparação do método bayesiano proposto aqui com outras técnicas de previsão de séries temporais, como, por exemplo, redes neuro-fuzzy (Ballini et al., 1999). Referências Akaike, H. (1974). A new look at the statistical identification model. IEEE Trans. on Automatic Control, 19:716-723. Andrade, M.G. e Barreto, G. A. (1999). Estimação e análise de robustez na abordagem bayesiana de modelos AR(p): Estudo de caso em previsão de vazão. Relatório Técnico. ICMC-USP. Ballini, R., Soares, S. e Andrade, M.G. (1999). Seasonal streamflow forecasting via a neural fuzzy system. Proc. 14th IFAC World Congress, Beijing, China, pp. 81-86. Box, G., Jenkins, G.M. e Reinsel, G. (1994). Time series analysis: Forecasting and control. 3rd ed., Prentice-Hall. Brooks, S.P. (1998). Markov chain Monte Carlo method and its application. The Statistician, 47:69-100. Davis, M.H.A. e Vinter, R. (1985) Stochastic Modelling and Control. Monographs on Statistics and Applied Probability. Chapman & Hall. Gelman, R., Carlin, J.B., Stern, H.S. e Rubin, D.B. (1995). Bayesian Data Analysis. Chapman & Hall. Geman, S. e Geman, D. (1984). Stochastic relaxation, Gibbs distributions and the bayesian restoration of images. IEEE Trans. on Pattern Analysis and Machine Intelligence.6:721-741. Pereira, M.V.F., Cunha, S.H.F., Terry, L.A. e Mossé, A. (1990). Modelos computacionais para planejamento e operação de sistemas hidrotérmicos de grande porte. Revista SBA: Controle e Automação, 1(1):31-41. Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2):461-464. Soares, S. (1990). Planejamento da operação de sistemas hidrotérmicos. Revista SBA: Controle e Automação, 1(2):122-127. Tweedie, R.L. (1998). Markov chains: Structure and applications. Technical Report, Department of Statistics, Colorado State University.