ads:
Universidade Federal de Campina Grande
Centro de Ciências e Teconologia
Programa de Pós-Graduação em Matemática
Curso de Mestrado em Matemática
Estudo de Modelos ARIMA com
Variáveis Angulares para Utilização
na Perfuração de Poços Petrolíferos
por
Areli Mesquita da Silva
sob orientação do
Prof. Dr. Francisco Antônio Morais de Souza
Dissertação apresentada ao Corpo Docente do Programa
de Pós-Graduação em Matemática - CCT - UFCG, como
requisito parcial para obtenção do título de Mestre em
Matemática.
Campina Grande - PB
Julho/2007
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
Estudo de Modelos ARIMA com
Variáveis Angulares para Utilização
na Perfuração de Poços Petrolíferos
por
Areli Mesquita da Silva
Dissertação apresentada ao Corpo Docente do Programa de Pós-Graduação em
Matemática - CCT - UFCG, como requisito parcial para obtenção do título de Mestre
em Matemática.
Área de Concentração: Probabilidade e Estatística
Aprovada por:
Prof. Dr. André Gustavo Campos Pereira
Prof. Dr. Antonio José da Silva
Prof. Dr. Francisco Antônio Morais de Souza
Orientador
Universidade Federal de Campina Grande
Centro de Ciências e Tecnologia
Programa de Pós-Graduação em Matemática
Curso de Mestrado em Matemática
Julho/2007
ads:
Agradecimentos
A Deus p o r mais essa dádiva em minha vida!
A meus pais, Eri e Manoel, pelo investimento e incentivo dados em todos os
momentos.
A Fúlvio (Vinho) pelo apoio, companheirismo e por sempre procurar deixar meu
ego nas alturas!
Ao meu orientador, professor Francisco Antônio Morais de Souza, por todos os
ensinamentos, pacientemente, compartilhados, sem os quais, teria sido inviável desen-
volver este trabalho.
À ANP (Agência Nacional do Petróleo, Gás e Biocombustíveis) e aos demais
orgãos financiadores pela concessão da bolsa.
Aos professores André Gustavo Campos Pereira e Antonio José da Silva por terem
aceito participar da banca.
Ao professor Brandão por suas brilhantes sugestões (veja como o Apêndice B
ficou lindo!).
A todos os professores de graduação e pós-graduação da UAME/UFCG que es-
tiveram sempre na torcida!
A todos os funcionários da UAME/UFCG que nunca economizam esforços na
hora de ajudar!
A Joelma (Joca) por acreditar que, um dia, daria tudo certo...(nunca esquecerei
da pergunta: E a integral?).
A Cris, Grayci (minhas irmãs acadêmicas), a Rosângela (Rosinha), Tatiana (Chaty),
Juliana, Jacqueline, Hallyson, Jesualdo (Nash), Josea ne, Leomaques,..., pelo carinho e
convivência.
A todos que, com simples gestos, contribuíram para que este trabalho foss e con-
cluído.
Dedicatória
A minha família.
“A grandeza de um ser humano não está no quanto ele
sabe, mas no quanto ele tem consciência que não sabe. O
destino não é freqüentemente inevitável, mas uma ques-
tão de escolha. Quem faz escolha, escreve sua própria
história, constrói seus próprios caminhos.”
Augusto Cury
Resumo
Séries temporais envolvendo dados angulares aparecem nas mais diversas áreas
do conhecimento. Por exemplo, na perfuração de um poço petrolífero direcional, o
deslocamento da broca de perfuração, ao longo da trajetória do poço, pode ser consi-
derado uma realização de uma série temporal de dados angulares. Um dos interesses,
neste contexto, consiste em realizar previsões de posicionamentos futuros da broca de
perfuração, as quais darão mais apoio ao engenheiro de petróleo na tomada de deci-
são de quando e como interferir na trajetória de um poço, de modo que este siga o
curso planejado. Neste trabalho, estudamos algumas classes de modelos que podem
ser utilizados para a modelagem desse tipo de série.
Abstract
Time series involving angular data appear in many diverse areas of scientific
knowledge. For example, in the drilling of a directional oil well, the displacement of
the drill, along the path of the well, can b e considered a s an a ngular data time series.
One of the objectives, in this context, consists in carrying out forecasts of the future
positions of the drill, which will give more support to the petroleum engineer in the
decision-making of when and how interfere in the path of a well, so that this follows
the planned course. In this work, we study some classes of models that can be utilized
for the modeling of that kind of series.
Sumário
Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1 Quantificação de Incertezas de Subsuperfície 4
1.1 Incertezas de Subsup erfície . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Impacto de Incertezas Dinâmicas sobre um Prog rama de Perfuração . . 6
1.3 Justificativa da Aquisição de Dados Complementares . . . . . . . . . . 7
2 Séries Temporais 9
2.1 Modelos Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.1 Modelos Auto-Regressivos . . . . . . . . . . . . . . . . . . . . . 16
2.1.2 Modelos de Médias Móveis . . . . . . . . . . . . . . . . . . . . . 18
2.1.3 Modelos Auto-Regressivos e de Médias Móveis . . . . . . . . . . 19
2.1.4 Modelos Auto-Regressivos Integrados e de Médias Móveis . . . . 20
2.2 A Função de Autocorrelação Parcial . . . . . . . . . . . . . . . . . . . . 22
2.3 Alguns Casos Particulares de Modelos Lineares . . . . . . . . . . . . . 24
2.3.1 Modelo Auto-Regressivo de Ordem 1 - AR(1) . . . . . . . . . . 24
2.3.2 Modelo Auto-Regressivo de Ordem 2 - AR(2) . . . . . . . . . . 24
2.3.3 Modelo de Médias Móveis de Ordem 1 - MA(1) . . . . . . . . . 26
2.3.4 Modelo de Médias Móveis de Ordem 2 - MA(2) . . . . . . . . . 26
2.3.5 Modelo Auto-Regressivo e de Médias Móveis de Ordem (1,1) -
ARMA(1,1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4 Identificação de Modelos ARIMA . . . . . . . . . . . . . . . . . . . . . 29
2.4.1 Procedimentos de Identificação . . . . . . . . . . . . . . . . . . 29
2.4.2 Estimativas Preliminares . . . . . . . . . . . . . . . . . . . . . . 34
ii
2.5 Estimação de Modelos ARIMA . . . . . . . . . . . . . . . . . . . . . . 35
2.5.1 Método dos Momentos . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.2 Método de Máxima Verossimilhança . . . . . . . . . . . . . . . 37
2.5.3 Variância dos Estimadores . . . . . . . . . . . . . . . . . . . . . 40
2.6 Diagnóstico de Modelos ARIMA . . . . . . . . . . . . . . . . . . . . . . 41
2.6.1 Teste de Autocorrelação Residual . . . . . . . . . . . . . . . . . 41
2.6.2 Teste de Box-Pierce . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.6.3 Teste da Autocorrelação Cruzada . . . . . . . . . . . . . . . . . 42
2.7 Previsão com Modelos ARIMA . . . . . . . . . . . . . . . . . . . . . . 43
2.7.1 Previsão de Erro Quadrático Médio (EQM) mínimo . . . . . . . 44
2.7.2 Formas Básicas de Previsão . . . . . . . . . . . . . . . . . . . . 45
2.7.3 Equação de Previsão . . . . . . . . . . . . . . . . . . . . . . . . 46
2.7.4 Atualização das Previsões . . . . . . . . . . . . . . . . . . . . . 47
2.7.5 Intervalos de Confiança . . . . . . . . . . . . . . . . . . . . . . . 47
3 Séries Temporais Envolvendo Dados Angulares 49
3.1 Modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.1 Processo Gaussiano Transformado . . . . . . . . . . . . . . . . . 50
3.1.2 Processo Arqueado . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.1.3 Processos Baseados em Funções de Ligação . . . . . . . . . . . . 52
3.2 Seleção do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.3 Identificação do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.4 Ajuste do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.4.1 Modelo Gaussiano Arqueado . . . . . . . . . . . . . . . . . . . . 55
3.4.2 Modelo Gaussiano Transformado . . . . . . . . . . . . . . . . . 56
3.4.3 Modelos com Ligação Direta e Inversa . . . . . . . . . . . . . . 57
Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Apêndices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
A Demonstração da Desiguadade (3.1) 59
A.1 Resultados Utilizados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
A.2 Demostração da Desigualdade (3.1) . . . . . . . . . . . . . . . . . . . . 59
iii
B Demonstração do Teorema (3.1) 61
B.1 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
B.2 Propriedades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
B.3 Demonstração do Teorema (3.1) . . . . . . . . . . . . . . . . . . . . . . 62
Bibliografia 71
Introdução
A perfuração de poços direcionais na indústria de petróleo, é uma técnica utilizada
de forma cada vez mais freqüente, tanto para atingir formaçõ es produtoras situadas
abaixo de locaçõ es verticalmente inacessíveis, como também para perfurar vários poços
a partir de um mesmo ponto [Thomas, 2001]. Sua utilização se dá, em particular, em
poços offshore.
A primeira etapa no projeto de um poço direcional é determinar o tipo de tra-
jetória a ser seguida para se atingir o alvo desejado, que pode ser uma formação com
acúmulo de hidrocarbonetos. Nessa etapa são levados em consideração os seguintes
elementos:
- A profundidade do(s) ponto(s) de mudança de trajetória;
- O afastamento horizontal;
- A direção/locação do objetivo;
- A profundidade vertical final do poço;
- As inclinações dos diversos trechos.
A mudança de orientação da trajetória do poço é uma operação dispendiosa que
envolve a retirada da coluna de perfuração e a introdução de uma ferramenta especial
contendo um motor de fundo [Lima], que tem a finalidade de iniciar a deflexão do poço e
orientá-lo para a direção desejada. Feita a deflexão, a ferramenta com o motor de fundo
é retirada e retorna-se com a coluna normal de perfuração, continuando até um próximo
desvio ou até atingir o alvo desejado (formação com acúmulo de hidrocarbonetos).
2
Sob o ponto de vista operacional, em cada mudança da direção do poço, a sua
nova orientação é feita a partir de informações obtidas em superfície, sobre a inclinação
e direção do poço [Thomas, 2001]. Essas informações podem ser enviadas pelo fluido
de perfuração ou através de um cabo elétrico e são registradas de forma contínua e
instantânea (no caso do cabo elétrico). É com base nessas informações que o engenheiro
de petróleo toma a decisão sobre interferências na trajetória do poço.
O deslocamento da broca de perfuração, ao longo da trajetória do poço, pode
ser visto como uma realização de uma série temp oral, onde a componente aleatória
corresponde à posição real da broca em cada momento. Por mais controle que se tenha
do processo, essa posição não é determinística, isto é, po de ser vista como uma variável
aleatória seguindo uma determinada distribuição de probabilidade.
Uma série temporal consiste de um co njunto de observações ordenadas no tempo
[Morettin e Toloi, 2004]. São exemplos de séries temporais:
- Cotações diárias do barril de petróleo;
- Índice de poluição de uma região produtora de petróleo;
- Registros de marés em um porto marítimo;
- Preços diários das ações de uma empresa de petróleo, por exemplo, a Petrobras.
Em geral, na análise de uma série temporal, estamos interessados em:
- Investigar o mecanismo gerador dessa série;
- Fazer previsões para va lores futuros da série;
- Procurar periodicidade relevantes nos dados.
Considerando essa abordagem do poço direcional como uma série temporal, o
nosso interesse consiste em estudar modelos adequados para fazer previsões de posi-
cionamentos futuros da broca, ou seja, previsões sobre a inclinação e direção do poço.
Essas previsões podem ser feitas a partir dos registros obtidos contínua e instantanea-
mente, além da litologia da rocha atravessada pelo próprio poço ou das litologias das
rochas atravessadas por outros poços do campo em desenvolvimento.
3
Como proposta de modelo para fazer as previsões, temos o modelo ARIMA (Auto-
Regressivo Integrado e de Médias Móveis), dada a sua ampla divulgação e utilização
[Morettin e Toloi, 2004].
Estatisticamente falando, estamos tratando de uma modelagem através de um
modelo ARIMA, aplicada à inclinação e direção de um poço direcional. O desafio
e a contribuição desse trabalho co nsiste na utilização de uma variável a ngular como
variável resposta, uma vez que na literatura são utilizadas variáveis lineares.
No Capítulo 1, discutimos sobre a quantificação de incertezas de subsuperfície,
bem como o impacto dessas incertezas sobre um programa de perfuração de poços
petrolíferos.
No Capítulo 2, apresentamos um estudo sobre séries temporais, enfatizando a
classe de modelos Auto-Regressivos Integrados e de Médias Móveis (ARIMA), utiliza-
dos para descrição, interpretação e previsão de séries temporais.
Dedicamos o Capítulo 3 a análise de séries temporais envolvendo dados circulares,
já que a variável de interesse é a posição da broca ao longo da trajetória de um poço
petrolífero.
Nos Apêndices, recordamos algumas definições e enunciamos os principais resul-
tados utilizados nas demonstrações.
Capítulo 1
Quantificação de Incertezas de
Subsuperfície
A modelagem de reservatórios é uma tarefa bastante árdua, devido à complexida-
de física envolvida na predição do escoamento e à dificuldade de se obter dados para a
modelagem [Charles et al., 2001]. Com isso, as atividades relacionadas à predição dos
parâmetros de interesse econômico, tais como: o volume total da rocha, a localização
do óleo, perfis de produção e estimativas de reservas, são de difícil realização.
Os responsáveis por tomadas de decisão devem realizar uma quantificação siste-
mática dos riscos técnicos associados a qualquer desenvolvimento recente. Além disso,
uma quantificação do impacto das diversas incertezas de subsuperfície (estruturais,
geológicas e dinâmicas) sobre os parâmetros de interesse econômico, pode auxiliar a
justificar a aquisição e processamento de mais dados e, com isso, reduzir a incerteza
inerente ao processo de tomada de decisão.
A seguir, discutiremos o valor da quantificação de incertezas de subsuperfície
no processo de tomada de decisões de investimentos, apresentaremos uma experiência
relativa ao impacto das incertezas dinâmicas sobre o programa de perfuração e, por
fim, veremos como utilizar a quantificação de incertezas estruturais para justificar a
aquisição de dados complementares.
5
1.1 Incertezas de Subsuperfície
Atualmente, existem ferramentas para a co nstrução de modelos geológico s 3D
e para a quantificação de incertezas sobre os parâmetros associados a esses modelos
[Charles et al., 2001]. Tais ferramentas têm facilitado a co mpreensão do impacto de
cada uma da s incertezas de subsuperfície sobre o campo de produção. Por exemplo,
existe uma cadeia de ferramentas desenvolvidas para lidar com a quantificação de
incertezas de subsuperfície cujo suporte é constituído po r três softwares principais, a
saber, ALEA, JACT A
T M
e EST.
A partir de mapas de incerteza produzidos por intérpretes sísmicos, o ALEA si-
mula diversos modelos estruturais do reservatório, calcula os correspondentes volumes
totais de rochas, além de exportar estas superfícies para o JACT A
T M
. Com isso,
é possível quantificar o impacto das incertezas oriundas, por exemplo, da conversão
tempo-profundidade sobre as incertezas associadas ao volume total da ro cha.
Após identificar os parâmetros que afetam as propriedades do reservatório, pode-
se simular diversos modelos geológicos de reservatórios, bem como calcular seus res-
pectivos volumes.
O software JACT A
T M
realiza uma combinação entre ambientes deposicionais,
tipos de rochas e simulações de parâmetros petrofísicos e permite que as incertezas que
afetam os parâmetros geo e statísticos sejam incorporadas. As realizações resultantes
podem ser visualizadas em 3D, analisadas e exp ortadas para um simulador de fluxo.
O software EST possibilita a simulaçã o de fluxo para cada realização geoestatís-
tica proveniente do JACT A
T M
, no entanto, estas simulações podem ser muito dispen-
diosas e incompletas, já que o ALEA e o JACT A
T M
consideram apenas as incertezas
estáticas (estruturais e geológicas), não levando em conta as incertezas dinâmicas, tais
como: permeabilidade relativa, transmissividade defeituosa, ou qualquer parâmetro
de fluxo. Assim, para minimizar a quantidade de operações de simulação de fluxo,
é necessário utilizar ferramentas e métodos para incorporar incertezas dinâmicas na
quantificação de incertezas associadas a perfis de produção ou a estimativas de reser-
vas.
Após identificar as principais incertezas de subsuperfície, pode-se utilizar os três
softwares ALEA, JACT A
T M
e EST para transformar essas incertezas de subsuperfície
6
em incertezas dos parâmetros de interesse, durante o proc esso de tomada de decisão.
Lembrando que a qualidade do resultado dependerá da confiabilidade das incertezas
de subsuperfície.
1.2 Impacto d e Incertezas Dinâmicas sobre um Pro-
grama de Perfuração
A seguir descreveremos uma experiência, relatada em [Charles et al., 2001], que
foi realizada em dois campos petrolíferos, denotados por campo X e campo Y (por
motivos relacionados a sigilo).
Um determinado campo X, desenvolvido recentemente, possui um campo satélite
Y, ambos de alta pressão e temperatura e separados por uma falha principal e por
uma falha tectônica contendo um fluido desconhecido. A produção do campo X pode
reduzir o volume de fluidos do campo Y. Tal redução precisa ser quantificada, já que
é impossível perfurar qualquer poço após uma depleção de 100 barris. Surge, então,
um questionamento: após o início da produção do campo X, por quanto tempo o
desenvolvimento do campo Y pode ser adiado?
Há duas alternativas para o desenvolvimentto do campo Y: ou a realização de
uma perfuração vertical a partir de uma nova plataforma ou uma perfuração direcional
a partir de uma plataforma já existente no campo X.
Como o interesse é modelar apenas a depleção, basta considerar somente as incer-
tezas dinâmicas, ou seja , não é necessário fazer uma representação das heterogeneidades
ou estruturas do campo Y.
Com a utilização de um modelo de simulação de fluxo, construído pelo operador,
foi possível realizar as seguintes atividades:
1 - Definição da variável resposta: depleção média em todas as camadas do campo
Y;
2 - Identificação dos principais parâmetros de incerteza: valor da permeabilidade
absoluta da falha tectônica e do aqüífero; permeabilidade relativa; tamanho do
aqüífero a oeste; variação de porosidade dentro do aqüífero a norte; principais
falhas de transmissividade; faixas de permeabilidade no reservatório; anisotropia
7
vertical e fluido acumulado na falha tectônica;
3 - Identificação da ordem de incerteza dos parâmetros selecionados;
4 - Estimação da função densidade de probabilidade (fdp) associada a cada parâme-
tro (tratado como variável aleatória);
5 - Utilização da metodologia experimental para identificar os parâmetros de im-
pacto mais significativo na depleção. Os únicos parâmetros que interferiram na
depleção foram: a permeabilidade relativa, o fluido acumulado na falha tectô-
nica e a transmissividade defeituosa. A partir da simulação da variável resposta,
construiu-se uma superfície resposta como função dos valores assumidos por esses
três parâmetros;
6 - Realização de uma simulação Monte Carlo utilizando, tanto as fdp’s associadas
aos parâmetros, como o modelo analítico da superfície de resposta, fornecendo
perfis de prováveis depleções.
Para cada uma das hipóteses (distribuição triangular da transmissividade defei-
tuosa e falha de escoa mento) foi construído um perfil de depleção. No primeiro caso,
ocorreu uma depleção de 100 barris, após 1,7 ano de produção no campo X, enquanto
que, no segundo caso, o tempo de depleção foi de 1,2 ano.
Com base no cenário mais pessimista, decidiu-se desenvolver o campo Y, a partir
de uma plataforma já existente no campo X. Associado a esta decisão, admitiu-se um
risco de 5% (após 1,3 ano de produção do campo X).
1.3 Justificativa da Aquisição de Dados Complemen-
tares
A partir da interpretação de dados sísmicos 2D migrados no tempo e de 50 poços
disponíveis em um campo maduro, construíram-se mapas de profundidade do nível do
reservatório [Charles et al., 2001]. Esses mapas facilitaram a estimativa do possível
volume total da rocha e a compreensão dos parâmetros de maior incerteza, tornando
possível justificar uma aquisição de dados sísmicos e definir um processamento mais
ajustado.
8
Os parâmetros de incerteza mais sig nificativos foram registrados e ordenados da
seguinte maneira:
- Campo de velocidade utilizado para a migração sísmica no tempo;
- Interpolação de dados sísmicos 2D e valores de poços;
- Conversão tempo-profundidade;
- OWC (ponto de contato água-óleo).
A combinação de todas estas incertezas gerou um intervalo de confiança em torno
da profundidade do topo do reservatório. Além disso, a simulação de 200 mapas de
possíveis profundidades em torno do mapa base, e dentro do interva lo de confiança,
resultou numa série de possíveis valores para o volume total da rocha. Porém, para
melhor estimar o potencial deste campo, as duas maiores incertezas deveriam ser re-
duzidas.
Portanto, mesmo quando se trata de um campo maduro com muitos poços, as
incertezas geométricas podem influenciar no volume total da rocha. Assim, os parâ-
metros a ser melhorados precisam ser identificados, a fim de reduzir essas incertezas.
Recomenda-se, também, a aquisição de dados sísmicos 3D e que esses tais dados sejam
processados utilizando-se uma migração de profundidade melhor do que uma migração
de tempo clássica.
Capítulo 2
Séries Temporais
Uma série temporal pode ser vista como um conjunto de observações Z
t
, geradas
sequencialmente no tempo [Box e Jenkins, 1976]. Fazemos referência ao parâmetro t
como sendo o tempo, mas a série Z
t
poderá ser função de algum outro parâmetro fí-
sico, como espaço , volume, profundidade, etc. Se o conjunto de instantes de tempo for
discreto (enumerável) ou não-enumerável, a série será discreta ou contínua, respectiva-
mente. De um modo mais formal, uma série temporal é uma realização ou trajetória
de um pro ce sso estocástico.
Um processo estocástico é uma família de variáveis aleatórias {Z
t
; t ∈ T } definidas
num mesmo espaço de probabilidades. Ou seja, para cada t ∈ T , Z
t
é uma variável
aleatória definida sobre o espaço amostral Ω. Portanto, Z
t
é uma função de dois
argumentos, Z(t, w), onde t ∈ T e w ∈ Ω.
Na Figura 2.1, podemo s observar que, para cada t ∈ T , Z(t, w) é uma variável
aleatória com uma distribuição de probabilidade. Por outro lado, para cada w ∈ Ω
fixado, obtemos uma função do tempo, ou seja, uma realização do processo.
São exemplos de séries temporais:
1- Valores diários de poluição numa região produtora de petróleo;
2- Preços diários das ações de uma empresa de petróleo;
3- Cotações diária s do barril de petróleo;
4- Rendimento anual per capita;
10
Figura 2.1: Processo estocástico como uma família de variáveis aleatórias.
5- Inflação mensal de uma determinada cidade;
6- Intensidade da corrente elétrica num dado ponto;
7- Intensidade do som num determinado local;
8- Registro de marés em um porto marítimo;
As séries 1 a 5 são discretas, enquanto que, as séries 6 a 8 são contínuas.
Os principais objetivos da análise de uma série temporal são
- Investigar o mecanismo gerador dessa série;
- Descrever o comportamento da série;
- Procurar periodicidades relevantes nos dados;
- Realizar previsões de valores futuros da série.
Para atingir esses objetivos, lançamos mão de modelos estocásticos (ou proba-
bilísticos). Uma classe importante de modelos estocásticos para descrição de séries
temporais é a dos modelos estacionários, que são baseados na hipótese de que o pro-
cesso permanece em equilíbrio em torno de um nível médio consta nte. Em outras
11
palavras, o processo evolui no tempo de modo que a escolha de uma origem dos tem-
pos não é importante, ou seja, as características de Z
t+k
, para todo k, são as mesmas
de Z
t
[Morettin e Toloi, 2004]. Desta forma, a média µ(t) e a variância V (t) de Z
t
são
constantes para todo t ∈ T , ou seja,
µ(t) = E[Z
t
] = µ e V (t) = Var[Z
t
] = E[(Z
t
− µ)
2
] = σ
2
.
A covariância entre Z
t
e Z
t+k
, é denominada função de autocovariância (facv), e
é definida por
γ
k
= Cov[Z
t
, Z
t+k
] = E[(Z
t
− µ)(Z
t+k
− µ)].
Pela própria definição de γ
k
, temos que γ
0
= Var[Z
t
] = σ
2
e, sendo o processo
estacionário, |γ
k
| → 0 quando k → ∞. Este comportamento pode ser observado na
Figura 2.2.
Figura 2.2: Representação da função de auto covariância.
Como a facv pode ser sensíve l às unidades em que são medidas as observações, é
comum utilizarmos a função de autocorrelação (fac), dada por
ρ
k
=
Cov[Z
t
, Z
t+k
]
Var[Z
t
]Var[Z
t+k
])
, k ∈ T.
Se o processo for estacionário, então a variância σ
2
= γ
0
é a mesma, tanto no
tempo t + k co mo em t. Assim,
ρ
k
=
γ
k
γ
0
=
γ
k
σ
2
, k ∈ T.
12
Observe que a fac é simétrica em torno do zero e ρ
k
= ρ
−k
, para todo k. A Figura
2.3 mostra a fac como um gráfico dos valores localizados nas diagonais da matriz de
autocorrelação.
Figura 2.3: Uma matriz de autocorrelação e a fac corresp ond ente.
A fac ρ
k
pode ser estimada através da expressão:
r
k
=
c
k
c
0
,
onde
c
k
=
1
N
N−k
t=1
(Z
t
− Z)(Z
t+k
− Z) , k = 0, 1, ..., N − 1
é a estimativa da função de autocovariância γ
k
e Z é a média amostral da série temporal.
Para que seja viável descrever uma série temp o ral através de modelos estacioná-
rios, devemos supor que tal série é estacionária. No entanto, na prática, a maioria das
séries que encontramos apresentam algum tipo de não-estacionariedade, por exemplo,
existem séries não-estacionárias quanto ao nível e outras quanto ao nível e à inclinação,
como mostram as Figuras 2.4 e 2.5. Outro tipo de não-estacionariedade é a explosiva,
13
que surge em séries que representam o crescimento de uma colônia de bactérias, por
exemplo.
Figura 2.4: Representação de uma série não-estacionária quanto ao nível.
Figura 2.5: Representação de uma série não-estacionária quanto ao nível e à inclinação.
Mais adiante abordaremos a classe de modelos ARIMA, que será útil para des-
crever de maneira satisfatória séries estacionárias e séries não-estacionárias que não
apresentam comportamento explosivo.
A fim de facilitar a manipulação dos modelos abo rdados mais adiante, utilizare-
mos o operador translação para o passa do, denotado por B e definido por
BZ
t
= Z
t−1
B
m
Z
t
= Z
t−m
, 2 ≤ m < t.
Mesmo quando uma série é não-estacionária, podemos transformar os dado s ori-
ginais, a fim de tentar obter uma série estacionária. O procedimento mais utilizado
consiste em diferenciar sucessivamente a série original, até se obter uma série estacio-
nária. Diferenciar, aqui, significa considerar diferenças sucessivas da série original.
A primeira diferença de Z
t
é definida por
∆Z
t
= Z
t
− Z
t−1
= Z
t
− BZ
t
= (1 − B)Z
t
,
14
onde B é o operador translação para o passado.
A segunda diferença é
∆
2
Z
t
= ∆[∆Z
t
] = ∆[Z
t
− Z
t−1
] = Z
t
− 2Z
t−1
+ Z
t−2
= (1 − 2B + B
2
)Z
t
= (1 − B)
2
Z
t
.
A n-ésima diferença de Z
t
é definida por
∆
n
Z
t
= ∆[∆
n−1
Z
t
] .
Em geral, pode-se considerar vários modelos diferentes para descrever o com-
portamento de uma série. No entanto, devemos utilizar critérios de comparaçã o en-
tre eles, a fim de escolher o modelo mais parcimonioso, ou seja, aquele com uma
quantidade mínima de parâmetros e que forneça previsões bastante precisas. A es-
colha do modelo a dequado baseia-se num ciclo iterativo do método de Box e Jenkins
[Morettin e Toloi, 2004], cujas etapas consistem em:
1- Fazer uma descrição da série, através do cálculo de estatísticas resumo e da
representação gráfica dos dados e, a partir daí, escolher uma classe de modelos
para a análise;
2- Identificar um modelo através da análise de autocorrelações, dentre o utros crité-
rios;
3- Estimar os parâmetros do modelo identificado;
4- Realizar uma análise de resíduos, a fim de verificar se o modelo ajustado é ade-
quado para fazer previsões de valores futuros da série.
Se o modelo identificado não for adequado, o ciclo deve ser repetido a partir da
etapa 2.
2.1 Modelos Lineares
Os modelos abordados a seguir são casos particulares de um modelo de filtro
linear. A principal suposição deste modelo é que a série temporal tenha sido gerada a
15
partir de um filtro linear, ilustrado na Figura 2.6, cuja entrada é um ruído branco a
t
,
ou seja, para cada t ∈ T , a
t
é uma variável aleatória com
E[a
t
] = 0, ∀t,
Var[a
t
] = σ
2
a
, ∀t,
E[a
t
a
s
] = 0, s = t.
Assim, a série pode ser expressa da seguinte maneira
Z
t
= µ + a
t
+ ψ
1
a
t−1
+ ψ
2
a
t−2
+ ···
= µ + ψ(B)a
t
, (2.1)
onde µ, em geral, é o pa râmetro que determina o nível da série e
ψ(B) = 1 + ψ
1
B + ψ
2
B
2
+ ···
é o operador linear, cuja finalidade é tranformar a
t
em Z
t
, denominado função de
transferência do filtro.
Figura 2.6: Série temp oral gerada por um filtro linear.
Quando a série de pesos ψ
1
, ψ
2
, . . . for finita ou infinita convergente, então Z
t
é
estacionária com média µ. Caso contrário, Z
t
é não-estacionária e µ não tem significado
específico [Morettin e Toloi, 2004].
Lembrando que a
t
é um ruído branco e supondo que
∞
k=0
ψ
2
k
< ∞, temos que a
facv de Z
t
pode ser escrita da seguinte maneira
γ
k
= σ
2
a
∞
k=0
ψ
i
ψ
i+k
,
com ψ
0
= 1. Assim, para k = 0, obtemos a variância de Z
t
,
γ
0
= Var[Z
t
] = σ
2
a
∞
k=0
ψ
2
k
.
16
A série
˜
Z
t
= Z
t
− µ, pode ser escrita como uma soma de valores passados mais
um ruído a
t
, ou seja,
˜
Z
t
= π
1
˜
Z
t−1
+ π
2
˜
Z
t−2
+ ··· + a
t
,
ou ainda,
˜
Z
t
− π
1
˜
Z
t−1
− π
2
˜
Z
t−2
− ··· = a
t
,
donde segue que
π(B)
˜
Z
t
= a
t
, (2.2)
onde π(B) = 1 − π
1
B − π
2
B
2
− ··· .
Comparando as expressões (2.1) e (2.2), temos que
π(B)ψ(B)a
t
= a
t
,
daí,
π(B) = ψ
−1
(B) , (2.3)
mostrando que os pesos π
k
podem ser obtidos a partir dos pesos ψ
k
e vice-versa.
Quanto às condições de estacionariedade e invertibilidade, um processo linear será
estacionário se a série ψ(B) convergir para |B| ≤ 1 e será invertível se π(B) convergir
para |B| ≤ 1 [Morettin e Toloi, 2004].
2.1.1 Modelos Auto-Regressivos
Considerando o caso especial de (2.2), em que π
k
= 0, k > p, e renomeando os pesos
de π
k
para φ
k
, obtemos o modelo auto-regressivo de ordem p, denotado por AR(p)
˜
Z
t
= φ
1
˜
Z
t−1
+ φ
2
˜
Z
t−2
+ ··· + φ
p
˜
Z
t−p
+ a
t
, (2.4)
ou equivalentemente,
φ(B)
˜
Z
t
= a
t
, (2.5)
onde
φ(B) = 1 − φ
1
B − φ
2
B
2
− ··· − φ
p
B
p
é chamado operador auto-regressivo de ordem p.
De (2.5) temos que
˜
Z
t
=
1
φ(B)
a
t
= φ
−1
(B)a
t
,
17
ou seja, o modelo AR(p) pode ser visto como a saída
˜
Z
t
de um filtro linear, com função
de transferência φ
−1
(B), desde que a entrada a
t
seja um ruído branco.
Para que o processo Z
t
seja estacionário, a série ψ(B) = φ
−1
(B) deve convergir
para |B| 1, ou seja, as raízes de φ(B) = 0 devem cair fora do círculo unitário.
Por outro lado, como a série π(B) = φ(B) = 1 −φ
1
B −φ
2
B
2
−···−φ
p
B
p
é finita,
conseqüentemente, π(B) é convergente para |B| 1, então não há restrições sobre os
parâmetros de um processo auto-regressivo para garantir a invertibilidade de Z
t
.
Para encontrar a fac de um processo AR(p), devemos, primeiramente, multiplicar
ambos os membros de (2.4) por Z
t−k
e, em seguida, calcular o valor esperado
E[
˜
Z
t
˜
Z
t−k
] = φ
1
E[
˜
Z
t−1
˜
Z
t−k
] + φ
2
E[
˜
Z
t−2
˜
Z
t−k
] + ···+ φ
p
E[
˜
Z
t−p
˜
Z
t−k
] + E[a
t
˜
Z
t−k
] (2.6)
Mas, para k > 0, temos E[a
t
˜
Z
t−k
] = 0, pois
˜
Z
t−k
envolve ruídos apenas até a
t−k
,
não-correlacionados. Com isso,
γ
k
= φ
1
γ
k−1
+ φ
2
γ
k−2
+ ··· + φ
p
γ
k−p
, k > 0 .
Assim, dividindo ambos os membros dessa expressão por γ
0
= Var[Z
t
], obtemos a fac
ρ
k
= φ
1
ρ
k−1
+ φ
2
ρ
k−2
+ ··· + φ
p
ρ
k−p
, k > 0 . (2.7)
Segundo Box e Jenkins, a fac de um processo AR(p), consiste de uma mistura de
exponenciais e senóides amortecidas.
A va riância do processo pode ser o btida fazendo k = 0 na expressão (2.6), obtendo
Var(
˜
Z
t
) = Var(Z
t
) = γ
0
= φ
1
γ
1
+ ··· + φ
p
γ
p
+ σ
2
a
.
Dividindo ambos os membros por γ
0
, obtemos
1 = φ
1
ρ
1
+ ··· + φ
p
ρ
p
+
σ
2
a
ρ
0
,
donde segue que
γ
0
= σ
2
a
/(1 − φ
1
ρ
1
− ··· − φ
p
ρ
p
) . (2.8)
Os parâmetros auto-regressivos φ
1
, . . . , φ
p
podem ser escritos em termos de ρ
1
,
ρ
2
, . . . , ρ
p
. Para tanto, basta substituir k = 1, 2, . . . , p em (2.7), obtendo um sistema
18
com p equações lineares, chamadas equações de Yule-Walker,
ρ
1
= φ
1
+ φ
2
ρ
1
+ ··· + φ
p
ρ
p−1
ρ
2
= φ
1
ρ
1
+ φ
2
+ ··· + φ
p
ρ
p−2
.
.
.
.
.
.
.
.
.
.
.
.
ρ
p
= φ
1
ρ
p−1
+ φ
2
ρ
p−2
+ ··· + φ
p
cuja representação matricial é
P
p
φ = ρ
p
,
onde
P
p
=
1 ρ
1
··· ρ
p−1
ρ
1
1 ··· ρ
p−2
.
.
.
.
.
.
.
.
.
.
.
.
ρ
p−1
ρ
p−2
··· 1
, φ =
φ
1
φ
2
.
.
.
φ
p
e ρ
p
=
ρ
1
ρ
2
.
.
.
ρ
p
.
Donde segue que,
φ = P
−1
p
ρ
p
. (2.9)
Utilizando a expressão (2.9) podemos estimar os coeficientes φ
1
, . . . , φ
p
, substi-
tuindo as fac teóricas ρ
k
por suas estimativas r
k
.
2.1.2 Modelos de Médias Móveis
Se ψ
k
= 0, k > q, na expressão (2.1), obtemos o modelo de médias móveis de ordem q,
denotado por MA(q). Renomeando os pesos de ψ
k
para −θ
k
, temos
Z
t
= µ + a
t
− θ
1
a
t−1
− θ
2
a
t−2
− ··· − θ
q
a
t−q
,
ou ainda,
˜
Z
t
= Z
t
− µ = (1 − θ
1
B − θ
2
B
2
− ··· − θ
q
B
q
)a
t
= θ(B)a
t
,
onde θ(B) = 1 − θ
1
B − θ
2
B
2
− ··· − θ
q
B
q
é chamado operador de médias móveis de
ordem q.
Como a série ψ(B) = θ(B) = 1 − θ
1
B − θ
2
B
2
−··· − θ
q
B
q
é finita, então não há
restrições sobre os parâmetros de um processo MA(q) para garantir a estacionariedade
de Z
t
.
19
A condição de invertibilidade para um processo MA(q) é que π(B) = θ
−1
(B)
convirja para |B| ≤ 1, isto é, as raízes de θ(B) = 0 devem cair fora do círculo unitário.
A facv de um modelo MA(q) é
γ
k
= E[(a
t
− θ
1
a
t−1
− ··· − θ
q
a
t−q
)(a
t−k
− θ
1
a
t−k−1
− ··· − θ
q
a
t−k−q
)]
= E
a
t
−
q
i=1
θ
i
a
t−i
a
t−k
−
q
j=1
θ
j
a
t−k−j
= E[a
t
a
t−k
] −
q
i=1
θ
i
E[a
t−k
a
t−i
] −
q
j=1
θ
j
E[a
t
a
t−k−j
] +
q
i=1
q
j=1
θ
i
θ
j
E[a
t−j
a
t−k−j
] .
Sabendo que
E[a
t
a
t−k
] =
σ
2
a
, k = 0
0, k = 0 ,
obtemos
γ
0
= Var[Z
t
] = (1 + θ
2
1
+ θ
2
2
+ ··· + θ
2
q
)σ
2
a
(2.10)
e
γ
k
=
(−θ
k
+ θ
1
θ
k+1
+ ··· + θ
q−k
θ
q
)σ
2
a
, k = 1, 2, ··· , q
0, k > q .
Donde segue que a fac de Z
t
é
ρ
k
=
−θ
k
+θ
1
θ
k+1
+···+θ
q−k
θ
q
1+θ
2
1
+θ
2
2
+···+θ
2
q
, k = 1, 2, ··· , q
0, k > q .
(2.11)
Ao contrário do que ocorre com um modelo AR(p), a fac de um modelo MA(q)
se anula para lags maiores do que q.
2.1.3 Modelos Auto-Regressivos e de Médias Móveis
Uma das maneiras de tornar um modelo mais parcimonioso, consiste em considerar,
simultaneamente, termos auto-regressivos e termos de médias móveis. Com isso, surge
uma classe de modelos mistos, denominados modelos auto-regressivos e de médias mó-
veis de ordem (p, q), denotados por ARMA(p,q)
˜
Z
t
= φ
1
˜
Z
t−1
+ ··· + φ
p
˜
Z
t−p
+ a
t
− θ
1
a
t−1
− ··· − θ
q
a
t−q
,
isto é,
φ(B)
˜
Z
t
= θ(B)a
t
,
20
onde φ(B) e θ(B) são os operadores auto-regressivos e de médias móveis, respectiva-
mente.
As condições de estacio nariedade e invertibilidade para um processo ARMA(p,
q) é que as raízes de φ(B) = 0 e de θ(B) = 0 caiam fora do círculo unitário.
A facv de um modelo ARMA(p, q) é
γ
k
= E{(φ
1
˜
Z
t−1
+ ··· + φ
p
˜
Z
t−p
+ a
t
− θ
1
a
t−1
− ··· − θ
q
a
t−q
)
˜
Z
t−k
}.
Lembrando que
˜
Z
t−k
depende apenas de choques a
t−k
, ocorridos até o tempo t−k,
temos que a covariância cruzada entre
˜
Z
t
e a
t
, definida p o r
γ
za
(k) = E[a
t
˜
Z
t−k
] ,
se anula para valores de k > 0 e é diferente de zero para k ≤ 0. Daí, a facv fica na
forma
γ
k
= φ
1
γ
k−1
+ ··· + φ
p
γ
k−p
+ γ
za
(k) − θ
1
γ
za
(k − 1) − ··· − θ
q
γ
za
(k − q) . (2.12)
Para k > q, obtemos
γ
k
= φ
1
γ
k−1
+ φ
2
γ
k−2
+ ··· + φ
p
γ
k−p
, k > q .
Portanto, a fa c do modelo é
ρ
k
= φ
1
ρ
k−1
+ φ
2
ρ
k−2
+ ··· + φ
p
ρ
k−p
, k > q ,
mostrando que as autocorrelações, para k > q, se comportam como nos modelos auto-
regressivos.
2.1.4 Modelos Auto-Regressivos Integrados e de Médias Móveis
A seguir, abordaremos uma classe de modelos apropriados para descrever séries tem-
porais não-estacionárias homogêneas, ou seja, séries que, apesar de não evoluirem
em torno de uma média constante ao longo do tempo, quando diferenciadas d ve-
zes, tornam-se estacionárias. Por exemplo, se a série for não-estacionária quanto ao
nível, então d = 1. Isto significa que basta ca lcular sua primeira diferença para torná-la
estacionária. Já séries não -esta cionárias quanto à inclinação, devem ser diferenciadas
duas vezes (d = 2), para obter a estacionariedade [Morettin e Toloi, 2004].
21
Se W
t
= ∆
d
Z
t
for estacionária, podemos representá-la através de um modelo
ARMA(p, q)
φ(B)W
t
= θ(B)a
t
.
Neste caso, dizemos que Z
t
é uma integral de W
t
, já que diferenciando Z
t
(no sentido
de diferenção sucessivas) obtemos W
t
. Dizemos, ainda, que Z
t
segue um modelo auto-
regressivo integrado de médias móveis de ordem (p, d, q), denotado p o r ARIMA(p, d,
q)
φ(B)∆
d
Z
t
= θ(B)a
t
. (2.13)
Sendo W
t
estacionária, então todas as raízes de φ(B) = 0 caem fora do círculo
unitário.
Uma forma alternativa de escrever a expressão (2.13) é
ϕ(B)Z
t
= θ(B)a
t
, (2.14)
em que
ϕ(B) = φ(B)∆
d
= φ(B)(1 − B)
d
é um operador auto-regressivo não-estacionário de ordem p + d, com d raízes sobre o
círculo unitário e as p restantes, fora do círculo unitário.
Com essa notação, o modelo ARIMA pode ser representado pela seguinte expres-
são
Z
t
= ϕ
1
Z
t−1
+ ··· + ϕ
p+d
Z
t−p−d
+ a
t
− θ
1
a
t−1
− ··· − θ
q
a
t−q
, (2.15)
que é denominada equação de diferenças, bastante útil para o cálculo de previsões.
Quando o interesse é calcular a variância dos erros de previsão, é conveniente
expressar o modelo ARIMA na forma de choques aleatórios, ou seja, em termos do
valor atual e prévios de a
t
, ou seja,
Z
t
= a
t
+ ψ
1
a
t−1
+ ψ
2
a
t−2
+ ··· = ψ(B)a
t
. (2.16)
Outra maneira de representar o modelo ARIMA é a forma invertida, que consiste
em expressar Z
t
em termos de seus valores prévios e do valor atual de a
t
, isto é,
Z
t
= π
1
Z
t−1
+ π
2
Z
t−2
+ ··· + a
t
. (2.17)
22
Às vezes, é útil considerar uma extensão do modelo ARIMA, acrescentando um termo
constante θ
0
na expressão (2.13), obtendo
ϕ(B)Z
t
= φ(B)∆
d
Z
t
= θ
0
+ θ(B)a
t
. (2.18)
Se θ
0
= 0, o modelo (2.18) pode ser usado para representar séries com tendências
estocásticas, ou seja , séries que apresentam mudanças aleatórias no nível e/ou na in-
clinação. Se θ
0
= 0, então o modelo (2.18 ) é capaz de representar séries com tendência
polinomial determinística de grau d. Além disso,
E(W
t
) = µ
w
= θ
0
/(1 − φ
1
− φ
2
− ··· − φ
p
) .
O modelo ARIMA é uma generalização dos modelos vistos anteriormente, já que
ARIMA(p, 0, 0) = AR(p) ,
ARIMA(0, 0, q) = MA(q) e
ARIMA(p, 0, q) = ARMA(p, q) .
2.2 A Função de Autocorrelação Parcial
A função de autocorrelação parcial (facp) é um instrumento bastante útil durante
a etapa de identificação do modelo a ser ajustado aos dados observados. Vejamos, a
seguir, como essa função é construída.
Denotando por φ
kj
o j-ésimo coeficiente de um modelo AR(k), temos que φ
kk
é
o último coeficiente. Utilizando essa notação, as equações de Yule-Walker po dem ser
escritas da seguinte maneira:
1 ρ
1
··· ρ
k−1
ρ
1
1 ··· ρ
k−2
.
.
.
.
.
.
.
.
.
.
.
.
ρ
k−1
ρ
k−2
··· 1
=
φ
k1
φ
k2
.
.
.
φ
kk
=
ρ
1
ρ
2
.
.
.
ρ
k
.
23
Resolvendo, sucessivamente, estas equações para k = 1, 2, . . . , obtemos
φ
11
= ρ
1
, φ
22
=
1 ρ
1
ρ
1
ρ
2
1 ρ
1
ρ
1
1
=
ρ
2
− ρ
2
1
1 − ρ
2
1
, φ
33
=
1 ρ
1
ρ
1
ρ
1
1 ρ
2
ρ
2
ρ
1
ρ
3
1 ρ
1
ρ
2
ρ
1
1 ρ
1
ρ
2
ρ
1
1
, ···
De modo geral, para ρ
kk
, a matriz no numerador é a mesma que a matriz no
denominador, exceto pela última coluna, que é substituída pelo vetor de autocorrelação
ρ
k
= (ρ
1
, . . . , ρ
k
)
t
.
A função de autocorrelação parcial é definida como sendo a quantidade φ
kk
, en-
carada como função de k.
Para um processo AR(p), a facp se anula para todos as defasagens maiores do
que p, isto é, o seu gráfico apresenta um "corte" após o defasagem p. Portanto, o
gráfico dessa função permite identificar o grau do polinômio auto-regressivo. Já para
o processo MA(q), a facp é dominada por uma mistura de exponenciais e/ou senóides
amortecidas. Tal comportamento é semelhante ao da fac de um processo AR(p). Por
fim, a facp de um processo ARMA(p,q), comporta-se de maneira similar à facp de um
processo MA puro [Morettin e Toloi, 2004].
Durante o estágio de identificação do modelo precisaremos calcular estimativas
das facp, a fim de compará-las com as respectivas facp teóricas. Por exemplo, no caso
dos modelos AR, tais estimativas podem ser feitas, ajustando-se, sucessivamente, pro-
cessos auto-regressivos de ordem p = 1, 2, 3, . . . por mínimos quadrados e considerando
as estimativas
ˆ
φ
11
,
ˆ
φ
22
,
ˆ
φ
33
, . . . do último coeficiente de cada ordem. A facp estimada
pode ser obtida, de modo alternativo, substituindo-se, nas equações de Yule-Walker,
as fac ρ
j
por suas estimativas r
j
, isto é,
r
j
=
ˆ
φ
k1
r
j−1
+
ˆ
φ
k2
r
j−2
+ ··· +
ˆ
φ
kk
r
j−k
, j = 1, . . . , k
e resolvendo-se essas e quações para k = 1, 2, . . . .
24
2.3 Alguns Casos Particulares de Modelos Lineares
2.3.1 Modelo Auto-Regressivo de Ordem 1 - AR(1)
O modelo AR(1) é dado por
˜
Z
t
= φ
1
˜
Z
t−1
+ a
t
,
ou equivalentemente,
φ(B)
˜
Z
t
= a
t
,
onde φ(B) = 1 − φ
1
B.
Para que o processo seja estacionário é necessário que −1 < φ
1
< 1.
Por (2.7), a fac de um processo AR(1) é da forma
ρ
k
= φ
1
ρ
k−1
, k > 0
cuja solução é
ρ
k
= φ
k
1
, k ≥ 0 .
Donde seg ue que, se φ
1
> 0, a fac decai exponencialmente e, caso φ
1
< 0, ela também
decai exponencialmente, alternando valores positivos e negativos. A Figura 2.7 ilustra
esse comportamento para φ
1
= 0, 8 e φ
1
= −0, 8.
Por (2.8) a variância de um processo AR(1) é
γ
0
=
σ
2
a
1 − ρ
1
φ
1
⇒ γ
0
=
σ
2
a
1 − φ
2
1
.
2.3.2 Modelo Auto-Regressivo de Ordem 2 - AR(2)
O modelo AR(2) é dado por
˜
Z
t
= φ
1
˜
Z
t−1
+ φ
2
˜
Z
t−2
+ a
t
,
ou ainda,
φ(B)
˜
Z
t
= a
t
,
onde φ(B) = 1 − φ
1
B − φ
2
B
2
.
Para que o processo seja estacionário é preciso que
φ
1
+ φ
2
< 1 , φ
2
− φ
1
< 1 , −1 < φ
2
< 1 .
25
Figura 2.7: Processos AR(1) e suas correspondentes funções de autocorrelação.
A fac de um processo AR(2) é
ρ
k
= φ
1
ρ
k−1
+ φ
2
ρ
k−2
, k > 0.
Substituindo p = 2 nas equações de Yule-Walker, obtemos
ρ
1
= φ
1
+ φ
2
ρ
1
ρ
2
= φ
1
ρ
1
+ φ
2
donde segue que
φ
1
= ρ
1
(1 − ρ
2
)/(1 − ρ
2
1
) e φ
2
= (ρ
2
− ρ
2
1
)/(1 − ρ
2
1
) .
Utilizando as equações de Yule-Walker, podemos também expressar ρ
1
e ρ
2
em
função de φ
1
e φ
2
, da seguinte maneira
ρ
1
= φ
1
/(1 − φ
2
) e ρ
2
= φ
2
+ φ
2
1
/(1 − φ
2
) .
A Figura 2.8 ilustra a fac de um processo AR(2) para φ
1
= 1, φ
2
= 0, 89 e
φ
1
= −1, φ
2
= −0, 89.
26
Figura 2.8: Funções de auto correlação para um processo AR(2).
Para obter a variância de um processo AR(2), basta substituir p = 2 em (2.8),
obtendo
γ
0
=
σ
2
a
1 − φ
1
ρ
1
− φ
2
ρ
2
.
2.3.3 Modelo de Médias Móveis de Ordem 1 - MA(1)
O modelo MA(1) é representado por
˜
Z
t
= a
t
− θ
1
a
t−1
= θ(B)a
t
,
onde θ(B) = 1 − θ
1
B
O processo é invertível se −1 < θ
1
< 1.
Substituindo q = 1 na expressão (2 .10), obtemos a variância do processo
γ
0
= (1 + θ
2
1
)σ
2
a
.
Utilizando (2.11), encontramos a função de autocorrelação
ρ
k
=
−θ
1
1+θ
2
1
, k = 1
0, k > 1 .
A Figura 2.9 apresenta a fac de um processo MA(1) para θ
1
= 0, 8.
2.3.4 Modelo de Médias Móveis de Ordem 2 - MA(2)
O modelo MA(2) é dado por
˜
Z
t
= a
t
− θ
1
a
t−1
− θ
2
a
t−2
= θ(B)a
t
,
onde θ(B) = 1 − θ
1
B − θ
2
B
2
.
27
Figura 2.9: Função de autocorrelação de um processo MA(1).
Para que o processo seja invertível é necessário que as raízes de θ(B) = 0 caiam
fora do círculo unitário, ou seja, devemos ter
θ
1
+ θ
2
< 1 , θ
2
− θ
1
< 1 , −1 < θ
2
< 1 .
Observe que essas condições são equivalentes às condições de estacionariedade
para um processo AR(2).
A partir de (2.10) e (2.11) o btemos
γ
0
= (1 + θ
2
1
+ θ
2
2
)σ
2
a
,
ρ
1
=
−θ
1
(1−θ
2
)
1+θ
2
1
+θ
2
2
,
ρ
2
=
−θ
2
1+θ
2
1
+θ
2
2
,
ρ
k
= 0 , k > 2 .
A Figura 2.10 apresenta a fac de um processo MA(2) para θ
1
= 0, 5, θ
2
= −0, 3.
Figura 2.10: Função de autocorrelação de um processo MA(2).
28
2.3.5 Modelo Auto-Regressivo e de Médias Móveis de Ordem
(1,1) - ARMA(1,1)
O modelo ARMA(1,1) é dado por
˜
Z
t
= φ
1
˜
Z
t−1
+ a
t
− θ
1
a
t−1
,
Ou equivalentemente,
φ(B)
˜
Z
t
= θ(B)a
t
,
onde φ(B) = 1 − φ
1
B e θ(B) = 1 − θ
1
B.
O processo é estacionário se −1 < φ
1
< 1 e invertível se −1 < θ
1
< 1.
A partir de (2.12) podemos obter
γ
1
= φ
1
γ
0
+ γ
za
(1) − θ
1
γ
za
(0) ,
γ
0
= φ
1
γ
1
+ γ
za
(0) − θ
1
γ
za
(−1) .
Mas,
γ
za
(1) = 0 ,
γ
za
(0) = E[a
t
˜
Z
t
] = E[a
t
(φ
1
˜
Z
t−1
+ a
t
− θ
1
a
t−1
)] = E[a
2
t
] = σ
2
a
γ
za
(−1) = E[a
t
˜
Z
t+1
] = E[a
t
(φ
1
˜
Z
t
+ a
t+1
− θ
1
a
t
)] = φ
1
E[a
t
˜
Z
t
] + E[a
t
a
t+1
] − θ
1
E[a
2
t
]
= φ
1
E[a
2
t
] − θ
1
E[a
2
t
] = (φ
1
− θ
1
)σ
2
a
.
Portanto,
γ
1
= φ
1
γ
0
− θ
1
σ
2
a
e γ
0
= θ
1
γ
1
+ σ
2
a
− θ
1
(φ
1
− θ
1
)σ
2
a
,
donde segue que
γ
0
=
(1 + θ
2
1
− 2φ
1
θ
1
)σ
2
a
1 − φ
2
1
e γ
1
=
(1 − φ
1
θ
1
)(φ
1
− θ
1
)σ
2
a
1 − φ
2
1
.
Para valores de k > 1, a fac do processo é
ρ
k
= φ
1
ρ
k−1
.
A Figura 2.11 ilustra a fac de um processo ARMA(1,1), com φ
1
= 0, 8 e θ
1
= 0, 3.
29
Figura 2.11: Função de autocorrelação de um processo ARMA(1,1).
2.4 Identificação de Modelos ARIMA
A identificação do particular modelo ARIMA a ser ajustado aos dados é uma
das etapas mais críticas do ciclo iterativo do método de Box e Jenkins, pois, vários
pesquisadores, usando a mesma série, podem identificar modelos diferentes.
O principal objetivo da identificação é encontrar os valores p, d e q do mo-
delo ARIMA(p,d,q), bem como determinar estimativas preliminares dos parâmetros,
as quais serão úteis durante o estágio de estimação.
2.4.1 Procedim entos de Identificação
A primeira etapa do processo de identificação consiste em verificar se é necessário
transformar a série original, a fim de estabilizar sua variância. Neste sentido, a trans-
formação de Box-Cox é bastante útil
Z
(λ)
t
=
Z
λ
t
− c
λ
λ = 0
logZ
t
λ = 0 ,
onde λ e c são parâmetros a serem estimados.
Para se ter uma noção do tipo de transformação a ser utilizada, pode-se construir
um gráfico que traz no eixo das abscissas, médias de subconjuntos de observações da
série original e no eixo das ordenadas, a amplitude de cada um desses conjuntos, isto
é, se Z
1
, Z
2
, . . . , Z
k
for um tal subconjunto, o gráfico será constituído por pontos da
forma (
Z, w), onde
Z =
1
k
k
i=1
Z
i
e w = max(Z
i
) − min(Z
i
) .
30
Se w independer de Z, os pontos desse gráfico ficarão espalhados em torno de
uma reta paralela ao eixo das abscissas; neste caso, não é necessário aplicar nenhuma
transformação à série original. Caso w seja diretamente proporcional a Z, a transfor-
mação logarítmica é adequada. A Figura 2 .12 apresenta alguns gráficos que podem
ocorrer na prática e os respectivos valores de λ.
Figura 2.12: Gráficos amplitude × média, ilustrando alguns valores possíveis de λ.
A segunda etapa do processo de identificação consiste em diferenciar a série, ob-
tida na primeira etapa, até conseguir sua estacionariedade, ou seja, até que o processo
W
t
= ∆
d
Z
t
se reduza a um ARMA(p,q). Uma maneira de saber a quantidade de dife-
renças, d, necessárias para tornar o processo estacionário consiste em observar quando
a fac amostral de W
t
decresce rapidamente para zero. Na prática, d = 0, 1 ou 2
[Morettin e Toloi, 2004].
A terceira etapa do processo de identificação consiste em analisar o comporta-
mento das autocorrelações e autocorrelações parciais estimadas, as quais devem re-
presentar adequadamente as respectivas quantidades teóricas desconhecidas. Através
dessa análise, devemos identificar o processo ARMA(p,q). A Tabela 2.1 apresenta um
resumo das principais características dos modelos mais usuais.
Na literatura, podemos encontrar outras propostas de identificação de modelos
ARMA(p,q). Existem, por exemplo, os métodos baseados em uma função penalizadora,
cuja idéia é escolher as ordens k e l que minimizem a seguinte quantidade
P (k, l) = lnˆσ
2
k,l
+ (k + l)
C(N)
N
,
31
Tabela 2.1 : Características das fac e facp de um processo ARIMA(p,d,q).
Ordem (1, d, 0) (0, d, 1)
comportamento de ρ
k
decai exponencialmente somente ρ
1
= 0
comportamento de φ
kk
somente φ
11
= 0 decaimento exponencial
dominante
estimativas iniciais φ
1
= ρ
1
ρ
1
= −θ
1
/(1 + θ
2
1
)
região de admissibilidade −1 < φ
1
< 1 −1 < θ
1
< 1
Ordem (2, d, 0) (0, d, 2)
comportamento de ρ
k
mistura de exponenciais
ou ondas senóides amor-
tecidas
somente ρ
1
= 0 e ρ
2
= 0
comportamento de φ
kk
somente φ
11
= 0 e φ
22
= 0 dominada por mistura de
exponenciais ou senóides
amortecidas
estimativas iniciais
φ
1
=
ρ
1
(1−ρ
2
)
1−ρ
2
1
,
φ
2
=
ρ
2
−ρ
2
1
1−ρ
2
1
ρ
1
= −
θ
1
(1−θ
2
)
1+θ
2
1
+θ
2
2
,
ρ
2
= −
θ
2
1+θ
2
1
+θ
2
2
região de admissibilidade
−1 < φ
2
< 1
φ
2
− φ
1
< 1
φ
2
+ φ
1
< 1
−1 < θ
2
< 1
θ
2
− θ
1
< 1
θ
2
+ θ
1
< 1
Ordem (1, d, 1)
comportamento de ρ
k
decai exponencialmente após o lag 1
comportamento de φ
kk
dominada por decaimento exponencial ap ós o lag 1
estimativas iniciais ρ
1
= (1 − φ
1
θ
1
)(φ
1
− θ
1
)/(1 + θ
2
1
− 2φ
1
θ
1
) , ρ
2
= ρ
1
φ
1
região de admissibilidade −1 < φ
1
< 1, −1 < θ
1
< 1
onde ˆσ
2
k,l
é uma estimativa da variância residual obtida ajustando um modelo ARMA(k,l)
às N observações da série e C(N) é uma função do tamanho da série.
Quando o número de parâmetros aumenta, o termo penalizador (k + l)
C(N )
N
au-
menta e a variância diminui. Portanto, minimizar P (k, l) é equivalente a identificar as
ordens k e l que equilibrem tal comportamento [Morettin e Toloi, 2004].
32
A seguir, citaremos alguns procedimentos de identificação baseados em funções
penalizadoras particulares.
- Critério de Informação de Akaike
Akaike (1973,19 74 ) propôs que as ordens k e l do modelo deveriam ser escolhidas
de modo a minimizar o seguinte critério
AIC(k, d, l) = Nlnˆσ
2
k,l
+
N
N − d
2(k + l + 1 + δ
d0
) + Nln2π + N , (2.19)
onde
δ
d0
=
1, d = 0
0, d = 0 ,
e ˆσ
2
k,l
é o estimador de máxima verossimilhança de σ
2
a
.
Se o interesse for comparar vários modelos, com N fixado, então os dois últimos
termos de (2.19) podem ser desconsiderados. Nestes casos, supondo d = 0, o critério
para determinação das ordens p e q, se reduz a
AIC(k, l) = N
lnˆσ
2
k,l
+
2
N
(k + l + 2)
, (2.20)
que ainda pode ser reescrito da seguinte maneira
AIC(k, l) = lnˆσ
2
k,l
+
2
N
(k + l) , (2.21)
já que os valores k e l que minimizam (2.21) são os mesmos que minimizam (2.20), pois
lnˆσ
2
k,l
+
2
N
(k + l)
< N
lnˆσ
2
k,l
+
2
N
(k + l + 2)
.
Para os modelos AR(p), o critério AIC se reduz a
AIC(k) = Nlnˆσ
2
k
+ 2k .
Com o intuito de diminuir a probabilidade de selecionar uma ordem maior do que
a verdadeira, Hurvich e Tsa i (1989) sugeriram uma correção para o AIC, dada por
AIC
c
(k) = AIC(k) +
2(k + 1)(k + 2)
N − k + 2
.
- Critério de Informação Bayesiano
Akaike (1977), Rissanem (1978) e Schwarz (1978), sugerem escolher o modelo
cujas ordens k e l minimizam o Critério de Informação Bayesiano, dado por
BIC(k, l) = lnˆσ
2
k,l
+ (k + l)
lnN
N
,
33
onde ˆσ
2
k,l
é a estimativa de máxima vero ssimilhança da variância residual do modelo
ARMA(k,l).
Para os modelos AR(p), o critério se reduz a
BIC(k) = lnˆσ
2
k
+
k
N
lnN .
- Critério de Hannan e Quinn
A proposta de Hannan e Quinn (1979) é minimizar a seguinte quantidade
HQC(k, l) = lnˆσ
2
k,l
+ 2(k + l)c
lnlnN
N
, c > 1 .
Para modelos AR(p), o critério pode ser escrito da seguinte forma
HQC(k) = lnˆσ
2
k
+ 2ck
lnlnN
N
, c > 1 .
- Critério FPE (Final Predictor Error)
Supondo que a série é representada por um modelo AR(p), Akaike (1969) propôs
minimizar a seguinte quantidade
FPE(k) =
1 +
2k
N
ˆσ
2
k
, µ conhecido
1 +
2k+1
N
ˆσ
2
k
, µ desconhecido ,
onde ˆσ
2
k
= c
0
−
k
j=1
ˆ
φ
j
c
j
.
Pode-se mostrar que o FPE é um estimador assintoticamente não-viciado e con-
sistente para o erro quadrático médio da previsão de Z
N+1
[Morettin e Toloi, 2004].
- Critério CAT (Criterion Autoregressive Transfer Function) - Método de Parzen
Este critério é baseado numa filosofia diferente das anteriores. Primeiramente,
deve-se assumir que o verdadeiro modelo é um AR(∞)
π(B)Z
t
= a
t
.
O próximo passo consiste em estimar a função de transferência π(B). Daí, a or-
dem selecionada ˆp é vista como uma aproximação finita ótima para o processo AR(∞).
34
A seleção de uma função de transferência ótima é feita a partir do valor de k que
minimiza a expressão
CAT(k) =
−
1 +
1
N
, k = 0
1
N
k
j=1
ˆσ
−2
j
− ˆσ
−1
k
, k = 1, 2, . . . ,
onde ˆσ
2
j
é a variância residual estimada para o modelo ajustado de ordem j.
2.4.2 Estimativas Prel im inare s
A seguir, veremos como obter, a partir das autocorrelações amostrais da série W
t
=
∆
d
Z
t
, estimativas preliminares dos parâmetros do modelo identificado, a s quais se-
rão utilizadas como valores iniciais para o processo iterativo de estimação de máxima
verossimilhança.
Para processos AR(p) devemos resolver as equações de Yule-Walker, substituindo
as autocorrelações teóricas ρ
k
por suas estimativas r
j
, com isso, obteremos
ˆ
φ
1
,
ˆ
φ
2
, . . . ,
ˆ
φ
p
.
Uma estimativa inicial da variância residual de um processo AR(p), pode ser
obtida substituindo-se, na expressão (2.8), γ
0
por c
0
, os φ
j
por
ˆ
φ
j
e o s ρ
j
por r
j
,
obtendo
ˆσ
2
a
= c
0
(1 −
ˆ
φ
1
r
1
−
ˆ
φ
2
r
2
− ··· −
ˆ
φ
p
r
p
) .
Para processos MA(q), estimativas iniciais para θ
1
, θ
2
, . . . , θ
q
, podem ser obtidas,
substituindo-se ρ
1
, . . . , ρ
q
por r
1
, . . . , r
q
na expressão (2.11) e resolvendo as q equações
não-lineares resultantes.
A variância residual pode ser estimada, inicialmente, através da expressão (2.10),
substituindo-se γ
0
por c
0
e os θ
j
por suas estimativas iniciais, obtendo
ˆσ
2
a
= c
0
/(1 +
ˆ
θ
2
1
+
ˆ
θ
2
2
+ ··· +
ˆ
θ
2
q
) .
Para os processos ARMA(p,q), resolvemos as p equações
ρ
k
= φ
1
ρ
k−1
+ ··· + φ
p
ρ
k−p
, k = q + 1, . . . , q + p ,
substituindo ρ
k
por r
k
, a fim de obter estimativas preliminares para φ
1
, . . . , φ
p
. Em
seguida, através da expressão (2.12), obtemos
ˆ
θ
1
,
ˆ
θ
2
, . . . ,
ˆ
θ
q
e ˆσ
2
a
.
35
Quando utilizamos o modelo ARIMA, com µ
w
= 0, isto é,
φ(B)W
t
= θ
0
+ θ(B)a
t
,
com µ
w
= θ
0
/(1−φ
1
−···−φ
p
), podemos obter uma estimativa inicial de θ
0
, substituindo
µ
w
por W e os φ
j
por
ˆ
φ
j
, obtendo
ˆ
θ
0
= W (1 −
ˆ
φ
1
− ··· −
ˆ
φ
p
) .
2.5 Estimação de Modelos ARIMA
Após identificar um modelo provisório a ser ajustado à série temporal, devemos
obter estimativas eficientes para os seus parâmetros.
Vamos denotar por ξ = (φ, θ, σ
2
a
) o vetor com os p + q + 1 parâmetros de um
modelo ARIMA(p,d,q), onde φ = (φ
1
, . . . , φ
p
) e θ = (θ
1
, . . . , θ
q
). A seguinte notação
também será útil: η = (φ, θ).
Suponha que a série o riginal Z = (Z
1
, Z
2
, . . . , Z
N
) tenha sido gerada por um
processo ARIMA(p,d,q). A partir daí, considerando d diferenças, podemos gerar uma
série W
t
estacionária: W = (W
1
, W
2
, . . . , W
n
), o nde W
t
= ∆
d
Z
t
e n = N − d. Com
isso, o problema de estimar os parâmetros do modelo ARIMA é equivalente a estimar
os parâmetros do modelo modelo ARMA(p,q) estacionário e invertível, representado
por
a
t
=
˜
W
t
− φ
1
˜
W
t−1
− φ
2
˜
W
t−2
− ··· − φ
p
˜
W
t−p
+ θ
1
a
t−1
+ θ
2
a
t−2
+ ··· + θ
q
a
t−q
, (2.22)
em que W
t
= ∆
d
Z
t
,
˜
W
t
= W
t
− µ
w
e µ
w
= E[W
t
] .
Quando d > 0, é conveniente considerar µ
w
= 0. Caso contrário, µ
w
será mais
um parâmetro a ser estimado.
A seguir, vamos descrever alguns métodos que possibilitam a obtenção de esti-
madores para os parâmetros do modelo identificado.
2.5.1 Método dos Momentos
O método dos momentos é um dos métodos de estimação mais simples e antigo. Este
método consiste em substituir, nas equações que relacionam as autocorrelações e os
36
parâmetros do modelo, os momentos teóricos (média, variância e autocorrelação) pelos
respectivos momentos amostrais e, em seguida, resolver as equações resultantes.
As estimativas preliminares descritas em 5.4.2 são obtidas através do método dos
momentos.
Para o modelo AR(p), o estimador de φ, pelo método dos momentos é dado por
ˆ
φ
MM
= (
ˆ
φ
1,MM
, . . . ,
ˆ
φ
p,MM
)
t
= R
−1
p
r
p
,
onde
R
p
=
1 r
1
r
2
. . . r
p−1
r
1
1 r
1
. . . r
p−2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
r
p−1
r
p−2
r
p−3
. . . 1
e r
p
= (r
1
, r
2
, . . . , r
p
)
t
.
Utilizando
ˆ
φ
MM
, podemos também estimar σ
2
a
, através do método dos momentos,
obtendo
ˆσ
2
MM
= c
0
(1 −
ˆ
φ
1,MM
r
1
− ··· −
ˆ
φ
p,MM
r
p
)
= c
0
(1 − r
t
p
φ
MM
) = c
0
(1 − r
t
p
R
−1
p
r
p
) .
Em particular, para p = 1, temos que
ˆ
φ
MM
=
ˆ
φ
1,MM
= r
1
e ˆσ
2
MM
= c
0
(1 − r
2
1
) .
Para o modelo MA(q), o estimador de
ˆ
θ, utilizando o método dos momentos, é
obtido resolvendo as equações
r
k
=
−
ˆ
θ
k,M M
+
ˆ
θ
1,MM
ˆ
θ
k+1,M M
+ ··· +
ˆ
θ
q−k,MM
ˆ
θ
q,MM
1 +
ˆ
θ
2
1,MM
+
ˆ
θ
2
2,MM
+ ··· +
ˆ
θ
2
q,MM
, k = 1, 2, . . . , q .
A variância residual estimada através do método dos momentos é
ˆσ
2
MM
= c
0
/(1 +
ˆ
θ
2
1,MM
+
ˆ
θ
2
2,MM
+ ··· +
ˆ
θ
2
q,MM
) .
Em particular, para q = 1, temos que
r
1
=
−
ˆ
θ
1,MM
1 +
ˆ
θ
2
1,MM
e ˆσ
2
MM
=
c
0
1 +
ˆ
θ
2
1,MM
.
Para o modelo ARMA(p,q), os parâmetros
ˆ
φ e
ˆ
θ, são estimados, através do
método dos momentos, em duas etapas:
37
(1) estimação de φ, através da solução
ˆ
φ
MM
= (φ
1,MM
, . . . , φ
p,MM
) da seguinte equa-
ção
r
k
=
ˆ
φ
1,MM
r
k−1
+ ··· +
ˆ
φ
p,MM
r
k−p
, k = q + 1, . . . , q + p ;
(2) estimação de θ, através da solução
ˆ
θ
MM
= (θ
1,MM
, . . . , θ
q,MM
) da equação (2.12),
utilizando as autocovariâncias amostrais c
k
e os estimadores
ˆ
φ
MM
obtidos na
etapa anterior.
Em particular, para p = q = 1, obtemos
r
2
=
ˆ
φ
1,MM
r
1
r
1
= c
1
/c
0
= (1 −
ˆ
φ
1,MM
ˆ
θ
1,MM
)(
ˆ
φ
1,MM
−
ˆ
θ
1,MM
)/(1 +
ˆ
θ
2
1,MM
− 2
ˆ
φ
1,MM
ˆ
θ
1,MM
) .
2.5.2 Método de Máxima Verossimilhança
Vamos denotar por f(z|ξ) a função densidade (ou de probabilidade) conjunta de Z =
(Z
1
, Z
2
, . . . , Z
N
). Fixado ξ, a função f(z|ξ) associa um determinado valor a cada
conjunto de observações z observado. Ago ra, quando fixamos z e variamos ξ, obtemos a
função de verossimilhança, denotada por L(ξ|z). Essa função é de grande importância
na teoria de estimação, devido ao "princípio da verossimilhança", que diz o seguinte:
dado que o modelo adotado é correto, toda a informação sobre ξ presente na amostra
está contida na função de verossimilhança; o s outros aspectos dos dados são irrelevantes
[Box e Jenkins, 1976]. Em geral, é conveniente trabalharmos com o logaritmo natural
de L(ξ|z), denota do por l(ξ|z) e denominado função de log-verossimilhança.
Os valores dos parâmetros que maximizam a função de verossimilhança (ou equi-
valentemente, a função de log-verossimilhança) são chamados estimadores de máxima
verossimilhança (EMV).
Observe que só é possível calcular os a
t
em (2.22) se tivermos valores iniciais para
os
˜
W ’s e para os a’s. Tais valores podem ser obtidos através de dois procedimentos:
um condicional e o outro incondicional.
- Procedimento Condicional
O procedimento condicional consiste em substituir os valores iniciais desconhe-
cidos por valores supostamente razoáveis, ou seja, supomos que são dados p valores
W
t
e q valores a
t
, que serão denotados por w
∗
e a
∗
, respectivamente. A partir daí, os
38
valores a
1
, a
2
, . . . , a
n
, condicionais à escolha dos valores iniciais w
∗
e a
∗
, poderão ser
calculados através da expressão (2.22).
Supondo que os a
t
’s são normalmente distribuídos, a função densidade conjunta
de a
1
, a
2
, . . . , a
n
é dada por
f(a
1
, a
2
, . . . , a
n
) =
n
t=1
f(a
t
) = (2π)
−n/2
(σ
a
)
−n
exp
−
n
t=1
a
2
t
/2σ
2
a
. (2.23)
Dada uma amostra particular w, a função de verossimilhança associada ao vetor
de parâmetros ξ e condicional à escolha de w
∗
e a
∗
, pode ser obtida a partir das
expressões (2.22) e (2.23)
L(ξ|w, w
∗
, a
∗
) = (2π)
−n/2
(σ
a
)
−n
exp
−
1
2σ
2
a
n
t=1
(
˜
W
t
− φ
1
˜
W
t−1
− ···−
φ
p
˜
W
t−p
+ θ
1
a
t−1
+ ··· + θ
q
a
t−q
)
2
.
Considerando o logaritmo de L, obtemos
l(ξ |w, w
∗
, a
∗
) = −
n
2
log(2π) − nlog(σ
a
) −
1
2σ
2
a
n
t=1
(
˜
W
t
− φ
1
˜
W
t−1
− ···−
φ
p
˜
W
t−p
+ θ
1
a
t−1
+ ··· + θ
q
a
t−q
)
2
.
Isto é,
l(ξ |w, w
∗
, a
∗
) ∝ −nlog(σ
a
) −
1
2σ
2
a
S(η|w, w
∗
, a
∗
) , (2.24)
onde
S(η|w, w
∗
, a
∗
) =
n
t=1
(
˜
W
t
− φ
1
˜
W
t−1
− ··· − φ
p
˜
W
t−p
+ θ
1
a
t−1
+ ··· + θ
q
a
t−q
)
2
=
n
t=1
a
2
t
(η|w, w
∗
, a
∗
) . (2.25)
é a soma de quadrados condicional.
Utilizando um asterisco para denotar l e S condicionais a w, w
∗
, a
∗
, podemos
escrever (2.24) e (2.25) da seguinte maneira
l
∗
(ξ) −nlog(σ
a
) −
1
2σ
2
a
S
∗
(η) ,
S
∗
(η) =
n
t=1
a
2
t
(η|w, w
∗
, a
∗
) .
Nosso interesse é maximizar l
∗
(ξ), que é equivalente a minimizar S
∗
(η). Portanto,
estimadores de máxima verossimilhança serão estimadores de mínimos quadrados e o
estudo de l
∗
(ξ) é equivalente ao de S
∗
(η).
Os valores iniciais w
∗
e a
∗
podem ser escolhidos de duas formas:
39
(1) um procedimento consiste em substituir os elementos de w
∗
e a
∗
por suas es-
peranças. Temos que E(a
t
) = 0 e, se o modelo não tiver parte determinística,
E(W
t
) = 0. Caso o modelo tenha parte determinística, s ubstituímos cada ele-
mento de w
∗
por w ;
(2) se o processo estiver próximo da não-estaciona riedade, ou seja, se alguma raiz
de φ(B) estiver próxima do círculo unitário, um pro ce dimento adequado consiste
em utilizar a expressão (2.22) para calcular a
p+1
, a
p+2
, . . . , colocando os valores
anteriores de a
t
iguais a zero.
Com isso, teríamos
a
p+1
=
˜
W
p+1
− φ
1
˜
W
p
− ··· − φ
p
˜
W
1
+ θ
1
a
p
+ ··· + θ
q
a
p−q+1
e assim por diante.
- Procedimento Não-Condicional
O procedimento não-condicional consiste em estimar os valores iniciais para os
˜
W ’s e para os a’s através de um método chamado backforecasting ("previsão para o
passado"), a fim de gerar valores antes do início da série.
Segundo [Morettin e Toloi, 2004], a função de log-verossimilhança não-condicional
pode ser aproximada por
l(ξ ) −nlogσ
a
−
S(η)
2σ
2
a
,
onde
S(η) = S(φ, θ) =
n
t=−∞
[a
t
(η, W)]
2
(2.26)
é a soma de quadrados não-condicional e
[a
t
(η, W)] = E(a
t
|η, W) . (2.27)
Pode-se obter boas aproximações para os estimadores de máxima verossimilhança
através dos estimadores de mínimos quadrados, obtidos minimizando-se a expressão
(2.26). Dado η, o cálculo da soma de quadrados (2.26) é feito através do cálculo das
esperanças condicionais (2.27) e através da expressão (2.22). Os valores [W
−j
] e [a
−j
],
j = 0, 1, 2, . . . são calculados utilizando-se o procedimento backforecasting.
40
Supondo que os W
t
’s tenham sido gerados por um processo ARIMA usual
φ(B)W
t
= θ(B)a
t
, (2.28)
então eles poderiam ter sido, igualmente, gerados pelo pro cesso
φ(F )W
t
= θ(F )e
t
, (2.29)
onde F é o operador translação para o futuro e e
t
é um ruído branco com a mesma va-
riância que a
t
[Box e Jenkins, 1976]. A representação (2.28) é chamada forma forward
do processo e a representação (2.29) é denominada forma backward. Assim, fazer
previsões antes que a série se inicie é equivalente a prever a série reversa.
2.5.3 Variância dos Estimadores
A precisão dos estimadores encontradas deve ser avaliada através da co nstrução de
intervalos de confiança para os parâmetros. Considerando o vetor de parâmetros η =
(φ, θ), cuja ordem é p + q.
Supondo n suficientemente grande, os estimadores de máxima verossimilhança
têm uma distribuição assintótica normal, isto é,
ˆ
η
D
−→ N
p+q
(η, V),
V = 2σ
2
a
∂
2
S(η)
∂η
2
1
···
∂
2
S(η)
∂η
1
∂η
k
.
.
.
.
.
.
.
.
.
∂
2
S(η)
∂η
k
∂η
1
···
∂
2
S(η)
∂η
2
k
. (2.30)
Além disso, o estimador de máxima verossimilhança de σ
2
a
é dado por
ˆσ
2
a
=
S(
ˆ
η)
n
e, para n suficientemente grande, ˆσ
2
a
e
ˆ
η são não-correlacionados [Morettin e Toloi, 2004].
As estimativas das variâncias dos estimadores e covariâncias entre os estimadores
são obtidas substituindo-se σ
2
a
em (2.30) por ˆσ
2
a
e calculando-se as derivadas
∂
2
S(η)
∂η
i
∂η
j
,
numericamente. Utilizando as estimativas das variâncias, podemos obter intervalos de
confiança para os parâmetros η
i
, i = 1, 2, . . . , p + q.
41
2.6 Diagnóstico de Modelos ARIMA
Após identificar o modelo e estimar seus parâmetros, devemos verificar se ele
representa, satisfatoriamente, os dados observados. Esta verificação pode ser feita
através de uma técnica chamada superajustamento, a qual consiste em estimar um
modelo com parâmetros extras e examinar, primeiramente, se eles são significativos e,
em seg uida, se a inclusão dos mesmos diminue significativamente a variância residual.
Para tanto, precisamos ana lisar os resíduos do modelo ajustado. Se o modelo ajustado
φ(B)W
t
= θ(B)a
t
,
com W
t
= ∆
d
Z
t
, for verdadeiro, então os "erros verdadeiros" a
t
= θ
−1
(B)φ(B)W
t
serão um ruído branco [Morettin e Toloi, 2004].
A seguir, descreveremos alguns testes de diagnósticos de um mo delo ajustado a
uma série temp o ral, baseados nas autocorrelações estimadas dos resíduos.
2.6.1 Teste de Autocorrelação Residual
Após estimar φ e θ, calculamos os resíduos estimados (ou simplesmente resíduos)
através da seguinte expressão
ˆa
t
=
ˆ
θ
−1
(B)
ˆ
φ(B)W
t
.
Se o modelo ajustado for adequado, os resíduos estimados ˆa
t
deverão estar pró-
ximos dos resíduos verdadeiros a
t
, conseqüentemente, deverão ser aproximadamente
não-correlacionados. Ou seja, denotando por ˆr
k
as autocorrelações dos resíduos ˆa
t
,
deveríamos ter ˆr
k
0. Em particular, supondo que o modelo ajustado é adequado,
deveríamos ter, aproximadamente,
ˆr
k
∼ N(0, 1/n) .
O cálculo das autocorrelações ˆr
k
é feito através da expressão
ˆr
k
=
n
t=k+1
ˆa
t
ˆa
t−k
n
t=1
ˆa
2
t
.
Para valores "grandes" de k, podemos obter uma indicação de uma possível
quebra de comportamento de ruído branco em a
t
, comparando ˆr
k
com os limites ±2/
√
n
[Morettin e Toloi, 2004].
42
2.6.2 Teste de Box-Pierce
Box e Pierce (1970) propuseram um teste bastante útil para indicar se os valores das
autocorrelações dos resíduos estimados são muito altos. Se o modelo for apropriado, a
estatística
Q(K) = n(n + 2)
K
j=1
ˆr
2
k
n − j
terá, aproximadamente, uma distribuição χ
2
com K − p − q graus de liberdade. Para
valores grandes de Q(K) rejeitamos a hipótese de ruído branco para os resíduos.
2.6.3 Teste da Autocorrelação Cruzada
Novos termos de médias móveis podem ser incluídos no modelo, a partir da verificação
das autocorrelaçõ es ˆr
k
. Por exemplo, se |ˆr
5
| > 2/
√
n, então um termo θ
5
a
t−5
deve
ser inserido no modelo. Uma maneira alternativa, consiste em investigar a função de
correlação cruzada (fcc), baseada na correlação cruzada entre valores passados da série
e o valor presente do ruído, e definida por
s
k
=
a
t
(Z
t−k
−
Z)
a
2
t
(Z
t
− Z)
2
, k = 1, 2, 3, . . .
Como os verdadeiros a
t
são desconhecidos, utilizamos os resíduos estimados ˆa
t
e subs-
tituímos s
k
por
ˆs
k
=
ˆa
t
(Z
t−k
− Z)
ˆa
2
t
(Z
t
− Z)
2
, k = 1, 2, 3, . . .
Se o modelo for apropriado, então a
t
e Z
t−k
devem ser não-correlacionados, para
k ≥ 1, ou seja,
Cov[a
t
, Z
t−k
] = γ
az
(k) = 0, k ≥ 1 .
Daí, se para um certo k
0
, s
k
0
assumir um valor "grande", o modelo deverá ser
considerado inadequado.
Se |s
k
| > 2/
√
n, então γ
az
(k) é significativamente diferente de zero. É razoável,
portanto, para k suficientemente grande, julgar s
k
significante quando |ˆs
k
| > 2/
√
n
[Morettin e Toloi, 2004].
Os resíduos podem ser utilizados para modificar o modelo da seguinte maneira:
se os resíduos b
t
do modelo ajustado
φ
0
(B)∆
d
0
Z
t
= θ
0
(B)b
t
(2.31)
43
não forem aleatórios, podemos utilizar o método de identificação visto na seção 5.4,
para descrevê-los através do modelo
φ(B)∆
d
b
t
= θ(B)a
t
. (2.32)
Daí, substituindo (2.32) em (2.31), obtemos um novo modelo que deverá ser ajustado
aos dados
φ
0
(B)φ(B)∆
d
0
∆
d
Z
t
= θ
0
(B)θ(B)a
t
,
cujos resíduos são aleatórios. Este ciclo de identificação, estimação e verificação deve
ser repetido até que um modelo adequado seja encontrado.
2.7 Previsão com Modelos ARIMA
Nas seções 5.4, 5.5 e 5.6 seg uimos as etapas do ciclo iterativo de identificação,
estimação e diagnóstico, com o objetivo de construir um modelo ARIMA(p, d, q) que
representasse adequadamente os dados observados. Agora, vamos utilizar esse modelo
para fazer previsões.
Supondo que temos as observações . . . , Z
t−2
, Z
t−1
, Z
t
, até o instante t, nosso in-
teresse é prever um valor Z
t+h
, h ≥ 1. Dizemos que t é a origem das previsões e h o
horizonte e denotamos por
ˆ
Z
t
(h) a previsão de Z
t+h
(ver Figura 2.13).
Figura 2.13: Observações de uma série temporal com previsões de origem t e horizonte h.
Primeiramente, vamos assumir que W
t
= (1 − B)
d
Z
t
é estacionário e invertível e
que os parâmetros do modelo são conhecidos.
44
Substituindo t por t + h nas expressões (2.1 5), (2.16) e (2.17), obtemos o modelo
ARIMA(p, d, q) nas três formas básicas:
(i) forma de equação de diferenças
Z
t+h
= ϕ
1
Z
t+h−1
+ ···+ ϕ
p+d
Z
t+h−p−d
−θ
1
a
t+h−1
−···−θ
q
a
t+h−q
+ a
t+h
; (2.33)
(ii) forma de choques aleatórios
Z
t+h
= a
t+h
+ψ
1
a
t+h−1
+ψ
2
a
t+h−2
+··· =
∞
j=0
ψ
j
a
t+h−j
=
t+h
j=−∞
ψ
t+h−j
a
j
; (2.34)
(iii) forma invertida
Z
t+h
= π
1
Z
t+h−1
+ π
2
Z
t+h−2
+ ··· + a
t+h
=
∞
j=1
π
j
Z
t+h−j
+ a
t+h
. (2.35)
2.7.1 Previsão de Erro Quadrático Médio (EQM) mínimo
Supondo que
ˆ
Z
t
(h) seja uma função linear das observações até o instante t, então, por
(2.34), também será uma função de a
t
, a
t−1
, . . . .
Indicando a melhor previsão por
ˆ
Z
t
(h) = ψ
∗
h
a
t
+ ψ
∗
h+1
a
t−1
+ ψ
∗
h+2
a
t−2
+ ··· =
∞
j=0
ψ
∗
h+j
a
t−j
,
nosso objetivo é encontrar os pesos ψ
∗
j
que minimizem o EQM de previsão, dado por
E[Z
t+h
−
ˆ
Z
t
(h)]
2
= E
∞
j=0
ψ
j
a
t+h−j
−
∞
j=0
ψ
∗
h+j
a
t−j
2
. (2.36)
Observando que
∞
j=0
ψ
j
a
t+h−j
=
∞
j=−h
ψ
h+j
a
t−j
, temos que o erro de previsão é
e
t
(h) = Z
t+h
−
ˆ
Z
t
(h) = ψ
0
a
t+h
+ ψ
1
a
t+h−1
+ ··· + ψ
h−1
a
t+1
−
∞
j=0
(ψ
h+j
− ψ
∗
h+j
)a
t−j
.
Substituindo essa última expressão em (2.36) e usando o fato de que os a
t
são não-
correlacionados, podemos reescrever o EQM de previsão da seguinte forma
E[e
t
(h)]
2
= (1 + ψ
2
1
+ ··· + ψ
2
h−1
)σ
2
a
+
∞
j=0
(ψ
h+j
− ψ
∗
h+j
)
2
σ
2
a
,
45
que é minimizado se ψ
∗
h+j
= ψ
h+j
, j = 0, 1, 2, . . . , h fixo. Assim, a previsão de EQM
mínimo é dada por
ˆ
Z
t
(h) = ψ
h
a
t
+ ψ
h+1
a
t−1
+ ψ
h+2
a
t−2
+ ··· =
∞
j=0
ψ
h+j
a
t−j
.
Conseqüentemente, o erro de previsão é
e
t
(h) = a
t+h
+ ψ
1
a
t+h−1
+ ··· + ψ
h−1
a
t+1
.
Logo,
Z
t+h
= e
t
(h) +
ˆ
Z
t
(h) , h ≥ 1 .
Utilizando a notação
[Z
t+h
] = E[Z
t+h
|Z
t
, Z
t−1
, . . . ] ,
temos que:
(a) a previsão de EQM mínimo é a esperança condicional de Z
t+h
, da das as observa-
ções passadas da série, ou seja,
ˆ
Z
t
(h) = [Z
t+h
];
(b) [e
t
(h)] = 0 e a va riância do erro de previsão é dada por
V (h) = (1 + ψ
2
1
+ ψ
2
2
+ ··· + ψ
2
h−1
)σ
2
a
; (2.37)
(c) os erros de previsão a um passo são não-correlacionados, pois
e
t
(1) = Z
t+1
−
ˆ
Z
t
(1) = a
t+1
;
(d) os erros de previsão para intervalos de tempo maiores que um são correlacionados,
bem como os erros de previsão para o mesmo horizonte h, de diferentes origens t
e t − j [Morettin e Toloi, 2004].
2.7.2 Formas Básicas de Previsão
A previsão
ˆ
Z
t
(h) pode ser calculada de três formas, utilizando as diversas representa-
ções do modelo ARIMA.
46
(i) Previsão utilizando a equação de diferenças
Considerando a esperança condicional em (2.33), temos que
ˆ
Z
t
(h) = ϕ
1
[Z
t+h−1
] + ··· + ϕ
p+d
[Z
t+h−p−d
]
−θ
1
[a
t+h−1
] − ··· − θ
q
[a
t+h−q
] + [a
t+h
] , h ≥ 1 ,
onde devemos usar os seguintes fatos:
[Z
t+k
] =
ˆ
Z
t
(k), k > 0,
[Z
t+k
] = Z
t+k
, k ≤ 0,
[a
t+k
] = 0, k > 0,
[a
t+k
] = a
t+k
, k ≤ 0,
(ii) Previsão utilizando a forma de choques aleatórios
Tomando a esperança condicional em (2.34), obtemos
ˆ
Z
t
(h) = ψ
1
[a
t+h−1
] + ψ
2
[a
t+h−2
] + ··· + ψ
h−1
[a
t+h
] + ψ
h
[a
t
] + ··· + [a
t+h
] .
(iii) Previsão utilizando a forma invertida
Considerando a esperança condicional em (2.35), temos que
ˆ
Z
t
(h) = π
1
[Z
t+h−1
] + π
2
[Z
t+h−2
] + ··· + [a
t+h
] .
2.7.3 Equação de Previsão
Por 2.7.2(i), a equação de previsão, vista como uma função de h, com origem t fixa,
satisfaz a equação de diferenças
ˆ
Z
t
(h) =
p+d
i=1
ϕ
i
ˆ
Z
t
(h − 1) , h > q ,
ou ainda,
ϕ(B)
ˆ
Z
t
(h) = (1 − B)
d
φ(B)
ˆ
Z
t
(h) = 0 , h > q ,
com ϕ(B) operando sobre h.
A função
ˆ
Z
t
(h), para h > q − p − d, consiste de uma mistura de polinômios,
exponenciais e senóides amortecidas [Morettin e Toloi, 2004].
47
2.7.4 Atualização das Previsões
Calculando as previsões de Z
t+h+1
a pa rtir de duas origens t + 1 e t, obtemos, respec-
tivamente,
ˆ
Z
t+1
(h) = ψ
h
a
t+1
+ ψ
h+1
a
t
+ ψ
h+2
a
t−1
+ ··· (2.38)
e
ˆ
Z
t
(h + 1) = ψ
h+1
a
t
+ ψ
h+2
a
t−1
+ ··· (2.39)
Subtraindo (2.39) de (2.38), temos que
ˆ
Z
t+1
(h) =
ˆ
Z
t
(h + 1) + ψ
h
a
t+1
.
Portanto, quando um novo dado for observado, podemos atualizar a previsão de
Z
t+h+1
, feita no instante t. Essa atualização consiste em prever o valor de Z
t+h+1
, na
origem t+1, adicionando à
ˆ
Z
t
(h+ 1) um múltiplo do erro de previsão a
t+1
= Z
t+1
−
ˆ
Z
t
(1).
2.7.5 Intervalos de Confiança
Para obtermos um intervalo de confiança para Z
t+h
, vamos supor que os erros satisfazem
as seguintes condições:
E[a
t
] = 0, ∀t,
Var[a
t
] = σ
2
a
, ∀t,
E[a
t
a
s
] = 0, s = t,
a
t
∼ N(0, σ
2
a
), ∀t.
Dado que conhecemos os valores passados e presente da série, Z
t
, Z
t−1
, Z
t−2
, . . . ,
a distribuição condicional de Z
t+h
será N(
ˆ
Z
t
(h), V (h)), onde V (h) é a variâ ncia do erro
de previsão, calculada através da expressão (2.37).
Assim, temos que
U =
Z
t+h
−
ˆ
Z
t
(h)
V (h)
∼ N(0, 1)
Portanto, fixado o coeficiente de confiança γ, é possível encontrar um valor u
γ
na
distribuição se U, tal que P(−u
γ
< U < u
γ
) = γ. Em outras palavras,
ˆ
Z
t
(h) − u
γ
V (h) ;
ˆ
Z
t
(h) + u
γ
V (h)
(2.40)
48
é um intervalo (aleatório) que contém Z
t+h
com probabilidade γ.
O cálculo de V (h) é feito, substituindo-se σ
2
a
por sua estimativa ˆσ
2
a
(obtida na
etapa de estimação dos parâmetros do modelo), ou seja,
V (h) = (1 + ψ
2
1
+ ψ
2
2
+ ··· + ψ
2
h−1
)ˆσ
2
a
= ˆσ
2
a
1 +
h−1
j=1
ψ
2
j
. (2.41)
Substituindo (2.41) em (2.40), obtemos
ˆ
Z
t
(h) − u
γ
ˆσ
a
1 +
h−1
j=1
ψ
2
j
≤ Z
t+h
≤
ˆ
Z
t
(h) + u
γ
ˆσ
a
1 +
h−1
j=1
ψ
2
j
.
Capítulo 3
Séries Temporais Envolvendo Dados
Angulares
Em diversas áreas de conhecimento aparecem dados da forma (θ
1
, t
1
), . . . , (θ
n
, t
n
),
onde θ
1
, . . . , θ
n
consistem de direções em tempos t
1
, t
2
, . . . , t
n
[Mardia e Jupp, 2000].
Em outras palavras, esses dados constituem uma série temporal de dados angulares
(circulares ou direcionais).
São exemplos de séries temporais de dados circulares:
1- Direção de ventos e correntes marinhas;
2- Direção de migrações de animais;
3- Posição da broca durante a perfuração de um poço petrolífero.
3.1 Modelos
Existem diversos modelos para descrição e análise de séries temporais de dados
angula res, muitos deles construídos a partir de modelos para séries temp o rais lineares.
A escolha do modelo mais adequado é feita em várias etapas: escolha de uma classe geral
de modelos; identificação; estimação dos parâmetros do modelo identificado; ajuste e,
por fim, segue a etapa de previsão.
A seguir faremos uma descrição de quatro classes de modelos para séries temporais
de dados angulares, propostos por [Fisher e Lee, 1994].
50
3.1.1 Processo Gaussiano Transformado
Seja {(X
t
, Y
t
); t ∈ T } um processo no plano, onde T é um conjunto de índices, então a
projeção radial sobre o círculo unitário gera um processo correspondente Θ
t
sobre esse
círculo, definido p o r
X
t
= R
t
cosΘ
t
, Y
t
= R
t
senΘ
t
.
Quando {(X
t
, Y
t
); t ∈ T } é um processo Gaussiano bivariado estacionário en-
tão Θ
t
tem uma distribuição Gaussiana transformada. Além disso, se {X
t
; t ∈ T }
e {Y
t
; t ∈ T } sã o duas realizações independentes de um processo Gaussiano estacio-
nário de média zero e variância unitária então Θ
t
tem distribuição uniforme circular
[Fisher e Lee, 1994].
O ajuste de tais modelos apresenta um problema de falta de dados, já que a
parte radial {R
t
; t ∈ T } de um processo Gaussia no transformado não é observada. No
entanto, esse problema pode ser contornado através da utilização do algoritmo EM, o
qual será abordado na seção 3.2.
A estrutura de correlação de um processo {Θ
t
; t ∈ T } pode ser quantificada atra-
vés de uma medida de correlação entre duas variáveis circulares Θ
t
e Φ
t
, denominada
coeficiente de correlação, introduzido por [Fisher e Lee, 1983] e definido por
ρ
τ
=
E[sen(Θ
1
− Θ
2
)sen(Φ
1
− Φ
2
)]
E[sen
2
(Θ
1
− Θ
2
)]E[sen
2
(Φ
1
− Φ
2
)]
,
onde (Θ
1
, Φ
1
) e (Θ
2
, Φ
2
) são realizações independentes de (Θ, Φ).
De modo análogo ao caso linear, pode-se mostrar que (vide Apêndice A)
−1 ≤ ρ
τ
≤ 1. (3.1)
Além disso, se Θ e Φ forem independentes, então ρ
τ
= 0. A seguir veremos um
resultado importante envolvendo correlação circular.
Teorema 3.1 Sejam (X
1
, Y
1
) e (X
2
, Y
2
) vetores aleatórios independentes com uma
distribuição normal bivariada com variâncias iguais a σ
2
e correlação ρ. Sejam Θ
1
e
Θ
2
variáveis aleatórias angulares definidas por
(X
i
, Y
i
) = R
i
(cos Θ
i
, senΘ
i
), i = 1, 2.
Então a correlação circular entre Θ
1
e Θ
2
é dada por
ρ
τ
=
π
2
16
ρ
2
(1 − ρ
2
)
2
2
F
1
3
2
,
3
2
, 2; ρ
2
2
,
51
onde
2
F
1
é a função hipergeométrica (vide Apêndice B).
A partir desse resultado pode-se definir a função de a utocorrelação do processo
{Θ
t
; t ∈ T } da seguinte maneira
ρ
τ
(k) =
π
2
16
ρ(k)
2
(1 − ρ
2
(k))
2
2
F
1
3
2
,
3
2
, 2; ρ
2
(k)
2
, (3.2)
onde ρ(k) é a função de autocorrelação comum dos processos X
t
e Y
t
.
3.1.2 Processo Arqueado
Seja X
t
uma série temporal univariada de dados lineares. O arqueamento de X
t
em
torno do círculo unitário gera uma série temporal arqueada Θ
t
definida por
Θ
t
= X
t
(mod 2π),
ou sej a, Θ
t
é o resto da divisão de X
t
por 2π. Assim, um processo linear {X
t
; t ∈ T },
que dá origem a um processo arqueado {Θ
t
; t ∈ T }, pode ser decomposto da seguinte
maneira
X
t
= Θ
t
+ 2πk
t
,
onde k
t
é um inteiro não observado. Desta forma, o ajuste desse modelo também
apresenta um problema de falta de dados, que poderá ser solucionado a través do uso
do algoritmo EM.
O arqueamento de um processo {X
t
; t ∈ T } auto-regressivo AR(p) produz um
processo auto-regressivo arqueado (wrapped), denotado por WAR(p).
Segundo [Fisher e Lee, 1983], se (X, Y ) segue uma distribuição normal bivariada
com variâncias σ
2
X
e σ
2
Y
e correlação ρ, então o co eficiente de correlação circular ρ
τ
entre Θ = X(mod 2π) e Φ = Y (mod 2π) é dado por
ρ
τ
=
senh(2ρσ
X
σ
Y
)
senh(2σ
2
X
)senh(2σ
2
Y
)
.
Daí, segue que, se {X
t
; t ∈ T } é um processo AR(p), então a f unção de autocorrelação
circular do processo WAR(p) {Θ
t
; t ∈ T } é dada por
ρ
τ
=
senh[2ρ
k
σ
2
/(1 − φ
1
ρ
1
− φ
p
ρ
p
)]
senh[2σ
2
/(1 − φ
1
ρ
1
− φ
p
ρ
p
)]
,
onde ρ
k
é a autocorrelação de defasagem k do processo {X
t
; t ∈ T }, σ
2
/(1−φ
1
ρ
1
−φ
p
ρ
p
)
é a variância do processo e σ
2
, φ
1
, . . . , φ
p
são os seus parâmetros.
52
3.1.3 Processos Baseados em Funções de Ligação
Uma forma alternativa de se construir um processo angular, a partir de um processo
linear, consiste em utilizar uma função de ligação, isto é, uma função g bijetiva e
monótona que projete a reta real sobre o intervalo (−π, π) e de modo que g(0) = 0.
Duas funções de ligação bastante utilizadas são
g(x) = 2tg
−1
(x)
e a função obtida a partir da probit, que é
g(x) = 2x{Φ(x) − 0, 5},
onde Φ(·) é a função de distribuição da normal padrão.
Se g é uma função de ligação e X uma variável linear então Θ = g(X) é uma
variável angular e, de modo recíproco, para uma variável angular Θ, X = g
−1
(Θ) é
uma variável linear.
A seguir veremos dois modelos circulares construídos a partir de um modelo linear
ARMA(p, q): o modelo com ligação direta - LARMA, o nde os ângulos transformados
seguem um modelo linear ARMA e o modelo com ligação inversa - IAR, utilizando
médias condicionais e a distribuição von Mises [Fisher, 1993].
3.1.3.1 Processo com Ligação Direta - LARMA
Sendo {X
t
; t ∈ T } um processo linear, g uma função de ligaç ão e µ um ponto no círculo
então o processo angular correspondente {Θ
t
; t ∈ T } é definido por
Θ
t
= g(X
t
) + µ.
Dizemos que um processo angular estacionário {Θ
t
; t ∈ T } com média direciona l
µ é um processo auto-regressivo e de médias móveis com ligação (LARMA) quando o
seu processo linear com ligação X
t
= g
−1
(Θ
t
−µ) for um processo ARMA(p,q). Neste
caso, a autocorrelação circular de Θ
t
é dada por
ρ
τ
(k) = ρ
τ
{g(X
t
), g(X
t+k
)},
onde ρ
τ
é o coeficiente de correlação circular definido na seção 3.1.1.
As principais vantagens de um mo delo linear com ligação são as seguintes:
53
- Se o processo {Θ
t
; t ∈ T } for estacionário então o processo resultante {X
t
; t ∈ T }
também é;
- Os parâmetros são facilmente estimados utilizando-se um software.
3.1.3.2 Processo com Ligação Inversa - IAR
Para processos AR, pode-se definir um modelo alternativo através do uso de distribui-
ções condicionais.
Se g é uma função de ligação, µ é um ponto do círculo, κ > 0 e ω
1
, . . . , ω
p
são
coeficientes reais, então o processo com ligação inversa - IAR(p) é definido por
Θ
t
|(θ
t−1
, . . . , θ
t−p
) ∼ VM(µ
t
, κ)
onde
µ
t
= µ + g{ω
1
g
−1
(θ
t−1
− µ) + . . . + ω
p
g
−1
(θ
t−p
− µ)},
ou seja, Θ
t
dado θ
t−1
, . . . , θ
t−p
segue uma distribuição von Mises com média direcional
µ
t
e parâmetro de concentração κ constante. Este processo é bas tante utilizado para
modelar séries dispersas (κ < 2).
3.2 Seleção do Modelo
Em geral, a escolha de uma classe de modelos deve ser feita levando-se em con-
sideração as características observáveis dos dados.
No contexto de séries temporais envolvendo dados ang ulares temos quatro clas-
ses de modelos à nossa disposição (vide seção 3.1), diferentemente do caso de séries
temporais lineares, em que usávamos uma única família de modelos, a saber, os auto-
regressivos integrados e de médias móveis (ARIMA).
A seguir descreveremos algumas características específicas de cada classe de mo-
delos, as quais auxiliarão na escolha do modelo mais adequado.
A principal característica do modelo Gaussiano transformado é a distribuição
uniforme de suas marginais. Caso o processo linear que está sendo transformado tenha
correlação alta, então o processo Gaussiano transformado tenderá a ocupar diferentes
arcos do círculo correspondendo ao grau de variações observado no processo linear
54
altamente correlacionado. Por outro lado, para graus baixos de correlação, o processo
circular tenderá a se espalhar sobre todo o círculo.
A aparência do processo arqueado é afetada tanta pela correlação quanto pela
variância. Se a variância for grande, o processo circular tenderá a se espalhar unifor-
memente sobre o círculo, enquanto que, para valores pequenos da variância, o processo
tenderá a ocupar um arco do círculo.
No caso do processo com ligação direta, valores altos do processo linear tendem
a ser transformados em torno do arco µ + π, gerando um "vazio" na forma como os
dados se distribuem em torno do círculo.
Por fim, para o processo com ligação inversa, quanto menor o valor do parâmetro
de concentração κ, os dados estarão espalhados mais uniformemente sobre o círculo e,
quanto maior for ω, os dados tenderão a permanecer sobre um pequeno arco em torno
da média direcional µ.
Com base nessas características, recomenda-se que, se os dados angulares estive-
rem distribuídos uniformemente sobre o círculo, então o modelo Gaussiano transfor-
mado ou o modelo arqueado são adequados, sendo que para este último, o processo
linear arqueado deverá apresentar variância alta. Se, pelo contrário, os dados tendem
a se agrupar em torno de um arco, então os modelos com ligação direta ou inversa são
mais adequados.
Vale salientar que todos as quatro classes de modelos são bastante flexíveis e
capazes de modelar uma diversidade de comportamento de séries temporais de dados
angulares.
3.3 Identificação do Modelo
Após selecionar uma classe de modelos, devemos identificar um modelo particular
pertencente a essa classe. Uma ferramenta bastante útil nessa etapa é o correlograma
amostral da série. No caso dos modelos Gaussiano transformado, arqueado e com
ligação inversa, o correlograma amostral pode ser calculado utilizando-se a seguinte
expressão para as autocorrelações amostrais
55
ˆρ
τ
(k) =
det
τ−k
t=1
X
t
X
t−τ
det
τ−k
t=1
X
t
X
t
det
τ
t=k+1
X
t
X
t
, (3.3)
onde X
t
= (cos θ
t
, sen θ
t
)
. A idéia é comparar o correlograma amostral com as au-
tocorrelações teóricas de vários mo delos, a fim de identificar um modelo experimental
que se ajuste satisfatoriamente aos dados.
Para o processo Gaussiano transformado, pode-se estimar a função de autocor-
relação ρ(k) do processo linear correspondente, a partir da resolução da equação (3.2)
com ρ
τ
(k) substituído por sua estimativa ˆρ
τ
(k), calculada por (3.3).
Para o processo arqueado, a autocorrelação pode ser estimada através da seguinte
expressão
ρ
τ
(k) = senh(2c
0
ˆρ
k
)/senh(2c
0
)
onde c
0
é a variância do pro cesso linear, a qual poderá ser estimada através da relação
¯
R = exp(−c
0
/2), em que
¯
R é o comprimento médio resultante de θ
t
.
No caso do modelo com ligação direta, o cálculo da média direcional da série
circular é uma estimativa preliminar ˆµ de µ. Utilizando essa estimativa, a série trans-
formada X
t
= g
−1
(Θ
t
− ˆµ) pode ser identificada da maneira usual, isto é, através dos
procedimentos de identificação de séries temporais lineares.
3.4 Ajuste do Modelo
3.4.1 Modelo Gaussiano Arqueado
O processo linear arqueado X
t
pode ser parametrizado através da média µ, da variância
do processo c
0
e das p primeiras cova riâncias, c
1
, . . . , c
p
. O vetor de estatísticas s
= (
X
t
,
X
2
t
, . . . ,
X
t
X
t−p
)
é conjuntamente suficiente para estimar o vetor de
parâmetros (µ, c
0
, c
1
, . . . , c
p
)
[Fisher e Lee, 1994].
Considerando uma série temporal, θ
1
, . . . , θ
τ
, o algoritmo EM deve ser utilizado
em duas etapas:
Passo E → dada uma estimativa γ
N
de (µ, c
0
, c
1
, . . . , c
p
)
calcula-se
s
N
= E[s|γ
N
, θ
1
, . . . , θ
τ
] ;
56
Passo M → calcula-se a estimativa atualizada γ
N+1
, bem como a solução de
s
N
= E[s|γ
N+1
].
A maior dificuldade na aplicação deste algoritmo encontra-se no cálculo da espe-
rança condicional apresentada no passo E. Para s e ter uma idéia dessa dificuldade,
basta observar a quantidade de cálculos computacionais necessários para determi-
nar E[
X
t
X
t−1
|θ
1
, . . . , θ
τ
]. Por exemplo, a esperança condicional de X
t
X
t+1
dados
θ
1
, . . . , θ
τ
é dada por
∞
m
1
=−∞
. . .
∞
m
τ
=−∞
(θ
t
+ 2πm
t
)(θ
t−1
+ 2πm
t−1
) f(θ
1
+ 2πm
1
, . . . , θ
τ
+ 2πm
τ
)
∞
m
1
=−∞
. . .
∞
m
τ
=−∞
f(θ
1
+ 2πm
1
, . . . , θ
τ
+ 2πm
τ
)
,
onde f á densidade do processo AR(p). Como, à medida que j cresce, a influência
de θ
t−j
sobre θ
t
decresce rapidamente, então E[X
t
X
t+1
|θ
1
, . . . , θ
τ
], no caso do processo
AR(1), pode ser aproximada por E[X
t
X
t+1
|θ
1
, θ
t−1
], cuja expressão é, computacional-
mente, mais tratável
∞
m
1
=−∞
∞
m
t−1
=−∞
(θ
t
+ 2πm
t
)(θ
t−1
+ 2πm
t−1
) f
t,t−1
(θ
t
+ 2πm
t
, . . . , θ
t−1
+ 2πm
t−1
)
∞
m
1
=−∞
∞
m
t−1
=−∞
f
t,t−1
(θ
t
+ 2πm
t
, . . . , θ
t−1
+ 2πm
t−1
)
,
onde f
t,t−1
é a densidade conjunta de (X
t
X
t+1
). Essa aproximação é aplicável apenas
para modelos AR de ordem baixa [Fisher e Lee, 1994].
3.4.2 Modelo Gaussiano Transformado
De modo análogo ao processo arqueado, podemos parametrizar os processos lineares X
t
e Y
t
em termos da média µ, da variância do processo c
0
e das p primeiras covariâncias,
c
1
, . . . , c
p
. Novamente, na a plicação do algoritmo EM, a maior dificuldade encontra-se
no cálculo da esperança condicional.
Considerando a representação polar (R
t
cos Θ
t
, R
t
senΘ
t
) do vetor (X
t
, Y
t
), temos
E[X
t
X
t−j
|θ
1
, . . . , θ
τ
] = E[R
t
cos Θ
t
R
t−j
cos Θ
t−j
|θ
1
, . . . , θ
τ
]
= cos θ
t
cos θ
t−j
E[R
t
R
t−j
|θ
1
, . . . , θ
τ
]
57
Se X e Y têm distribuição normal τ-variada com média zero e matriz de covari-
âncias Σ e se (R
i
, θ
i
) á representação polar de (X
i
, Y
i
), então a esperança condicional
de R
t
R
t−j
dados θ
1
, . . . , θ
τ
é
E[R
t
R
t−j
|θ
1
, . . . , θ
τ
] =
∞
0
···
∞
0
r
t
r
t−j
τ
s=1
r
s
exp
−
1
2
r
Ar
dr
1
. . . dr
τ
∞
0
···
∞
0
τ
s=1
r
s
exp
−
1
2
r
Ar
dr
1
. . . dr
τ
,
onde A = (a
ij
), a
ij
= σ
ij
cos(θ
i
−θ
j
) e Σ
−1
= (σ
ij
) [Fisher e Lee, 1994]. Para processos
AR de ordem baixa pode-se aproximar E[R
t
R
t−j
|θ
1
, . . . , θ
τ
] por E[R
t
R
t−j
|θ
t
, θ
t−j
]. Em
particular, para p = 1, E[R
t
R
t−1
|θ
1
, . . . , θ
τ
] pode ser aproximada por
∞
0
···
∞
0
r
2
1
r
2
2
exp
−
1
2σ
2
{r
2
1
+ r
2
2
− 2φ
1
r
1
r
2
cos(θ
t
− θ
t−1
)}
dr
1
dr
2
∞
0
···
∞
0
r
1
r
2
exp
−
1
2σ
2
{r
2
1
+ r
2
2
− 2φ
1
r
1
r
2
cos(θ
t
− θ
t−1
)}
dr
1
dr
2
.
3.4.3 Modelos com Ligação Dir eta e Inversa
Após estimar a média através da média direcional amostral da série, pode-se utilizar
um software para realizar o ajuste do processo AR(p) com ligação direta.
Para ajustar processos AR(p) com ligação inversa, assume-se uma distribuição
marginal para Θ
1
, . . . , Θ
p
. Vamos supor que esta distribuição seja o produto de distri-
buições VM(µ, κ).
A densidade conjunta de θ
1
, . . . , θ
τ
é dada por
τ
t=p+1
f
θ
t
− µ − g[ ω
1
g
−1
{(θ
t−1
− µ)/2} + . . . + ω
p
g
−1
{(θ
t−p
− µ)/2}]
p
t=1
f(θ
t
− µ) ,
onde f é a densidade VM(0, κ).
Conclusão
A análise se séries temporais envolvendo dados angulares é feita de maneira li-
geiramente semelhante à análise de séries temporais envolvendo dados lineares. Vários
modelos utilizados para descrever, analisar e interpretar séries temporais angulares são
obtidos a partir de modelos para séries temporais envolvendo dados lineares, daí a ne-
cessidade de se ter domínio da teoria desenvolvida para a análise de séries temporais
lineares. Os modelos estudados têm ampla aplicação na perfuração de poços petrolífe-
ros direcionais e podem ser utilizados para modelar a posição da broca de perfuração
ao longo da trajetória de um poço petrolífero, considerada como uma série temporal
de dados angulares.
Apêndice A
Demonstração da Desiguadade (3.1)
Para demonstrar essa desigualdade utilizamos dois resultados, enunciados a se-
guir, cujas demonstrações podem ser encontradas em [James, 2002].
A.1 Resultados Utilizados
Desigualdade de Jensen - Seja X uma variável aleatória com E(X) < ∞ e g uma
função convexa. Então
g(E(X)) ≤ E(g(X)).
Desigualdade de Cauchy-Schwarz - Sejam X e Y variáveis aleatórias com segundo
momento finito. Então
E(|XY |) ≤
E(X
2
)E(Y
2
).
A.2 Demostração da Desigualdade (3.1)
Podemos expressar o coeficiente de correlação da seguinte forma
ρ
τ
=
E(XY )
E(X
2
)E(Y
2
)
,
com X = sen(Θ
1
− Θ
2
) e Y = sen(Φ
1
− Φ
2
).
Pela desigualdade de Jensen temos que
|E(XY )| ≤ E(|XY |) , (A.1)
60
já que E(XY ) < ∞ e a função modular é convexa.
Por outro lado, usando a desigualdade de Cauchy-Schwarz, obtemos
E(|XY |) ≤
E(X
2
)E(Y
2
) . (A.2)
De (A.1) e (A.2), segue que
|E(XY )| ≤
E(X
2
)E(Y
2
),
daí,
−
E(X
2
)E(Y
2
) ≤ E(XY ) ≤
E(X
2
)E(Y
2
) ,
logo,
−1 ≤
E(XY )
E(X
2
)E(Y
2
)
≤ 1 ,
ou seja,
−1 ≤ ρ
τ
≤ 1 .
Apêndice B
Demonstração do Teorema (3.1)
Antes de demonstrar esse teorema, vamos apresentar algumas definições e pro-
priedades, que serão úteis durante a demonstração.
B.1 Definições
Definição 1 A função densidade de probabilidade de um vetor aleatório (X,Y) com
distribuição normal padrão bivariada é dada por
f(x, y) =
1
2πσ
2
(1 − ρ
2
)
1/2
exp
−
1
2(1 − ρ
2
)
x
2
+ y
2
− 2ρxy
,
onde ρ é a correlação entre X e Y.
Definição 2 A função gama, é denotada por Γ(·) e definida por
Γ(α) =
∞
0
t
α−1
exp(−t)dt , α > 0 .
Definição 3 A função hipergeométrica, com argumentos α, β, γ e x, é definida por
2
F
1
(α, β, γ; x) =
∞
m=0
α
[m]
β
[m]
γ
[m]
x
m
m!
= 1 +
αβ
γ
x
1!
+
α(α + 1)β(β + 1)
γ(γ + 1)
x
2
2!
+ . . . ,
onde
h
[m]
=
1 se m = 0 ,
m−1
i=0
(h + i) se m > 0 .
62
B.2 Propriedades
1) A função gama satisfaz
(i) Γ(α + 1) = αΓ(α), α > 0 ;
(ii) Γ(n) = (n − 1)! , para qualquer inteiro positivo n ;
(iii) Γ
1
2
=
√
π .
2) Para quaisquer inteiros positivos p e q valem as seguintes fórmulas de integração
(i)
sen
p
x cos
q
x dx = −
sen
p−1
x cos
q+1
x
p + q
+
p−1
p+q
sen
p−2
x cos
q
x dx + c ;
(ii)
cos
q
x dx =
cos
q−1
x senx
q
+
q−1
q
cos
q−2
x dx + c .
B.3 Demonstração do Teorema (3.1)
Seja (R
i
, Θ
i
) a representação polar de (X
i
, Y
i
), ou seja,
(X
i
, Y
i
) = (R
i
cos Θ
i
, R
i
sen Θ
i
), i = 1, 2 .
Como a distribuição de Θ
1
e Θ
2
é uniforme, podemos escrever ρ
τ
da seguinte
maneira
ρ
τ
= 4(E[senΘ
1
senΘ
2
]E[cosΘ
1
cosΘ
2
] − E[senΘ
1
cosΘ
2
]E[senΘ
1
cosΘ
2
]) .
Para calcular cada valor esperado, precisamos encontrar a função densidade con-
junta de (R
1
, Θ
1
, R
2
, Θ
2
). Para tanto, utilizamos o método do jacobiano.
Sendo X = R cos Θ e Y = R sen Θ, f
(X,Y )
e f
(R,Θ)
as densidades de (X, Y ) e
(R, Θ), respectivamente, então
f
(R,Θ)
(r, θ) = f
(X,Y )
(rcosθ, rsenθ) |J|,
onde J é o jacobiano da transformação, calculado da seguinte forma
J =
∂x
∂r
∂x
∂θ
∂y
∂r
∂y
∂θ
=
cos θ −rsenθ
senθ r cos θ
= r cos
2
θ + r sen
2
θ = r ,
63
assim,
f
(R,Θ)
(r, θ) =
r
2π(1 − ρ
2
)
1/2
exp
−1
2(1 − ρ
2
)
[r
2
cos
2
θ + r
2
sen
2
θ − 2ρ r
2
cos θ sen θ]
=
r
2π(1 − ρ
2
)
1/2
exp
−1
2(1 − ρ
2
)
[r
2
− 2ρ r
2
cos θ sen θ]
.
Por hipótese, (X
1
, Y
1
) e (X
2
, Y
2
) são independentes, donde segue que (R
1
, Θ
1
)
e (R
2
, Θ
2
) também são independentes. Daí, a densidade conjunta de (R
1
, Θ
1
, R
2
, Θ
2
)
pode ser encontrada da seguinte maneira
f
(R
1
,Θ
1
,R
2
,Θ
2
)
(r
1
, θ
1
, r
2
, θ
2
) = f
(R
1
,Θ
1
)
(r
1
, θ
1
) f
(R
2
,Θ
2
)
(r
2
, θ
2
)
=
r
1
2π(1 − ρ
2
)
1/2
exp
−1
2(1 − ρ
2
)
[r
2
1
− 2ρ r
2
1
cos θ
1
sen θ
1
]
×
r
2
2π(1 − ρ
2
)
1/2
exp
−1
2(1 − ρ
2
)
[r
2
2
− 2ρ r
2
2
cos θ
2
sen θ
2
]
.
Após algumas manipulações algébricas, obtemos
f
(R
1
,Θ
1
,R
2
,Θ
2
)
(r
1
, θ
1
, r
2
, θ
2
) =
r
1
r
2
(2π)
2
(1 − ρ
2
)
exp
−1
2(1 − ρ
2
)
[r
2
1
+ r
2
2
cos(θ
1
− θ
2
)]
=
r
1
r
2
(2π)
2
(1 − ρ
2
)
exp
−r
2
1
2(1 − ρ
2
)
exp
−r
2
2
2(1 − ρ
2
)
×exp
ρr
1
r
2
cos(θ
1
− θ
2
)
(1 − ρ
2
)
=
r
1
r
2
(2π)
2
(1 − ρ
2
)
exp
−r
2
1
2(1 − ρ
2
)
exp
−r
2
2
2(1 − ρ
2
)
×
∞
n=0
1
n!
r
n
1
r
n
2
ρ
(1 − ρ
2
)
cos(θ
1
− θ
2
)
n
=
1
(2π)
2
(1 − ρ
2
)
∞
n=0
1
n!
r
n+1
1
exp
−r
2
1
2(1 − ρ
2
)
×exp
−r
2
2
2(1 − ρ
2
)
ρ
(1 − ρ
2
)
n
cos
n
(θ
1
− θ
2
)
• Cálculo de E(senΘ
1
cos Θ
2
)
Por definição de valor esperado, temos que
E(senΘ
1
cos Θ
2
) =
π
−π
π
−π
senθ
1
cos θ
2
f
(Θ
1
,Θ
2
)
(θ
1
, θ
2
)dθ
1
dθ
2
=
π
−π
π
−π
senθ
1
cos θ
2
∞
0
∞
0
f
(R
1
,Θ
1
,R
2
,Θ
2
)
(r
1
, θ
1
, r
2
, θ
2
)dr
1
dr
2
dθ
1
dθ
2
64
= C
∞
n=0
1
n!
∞
0
r
n+1
1
exp
−r
2
1
2(1 − ρ
2
)
dr
1
∞
0
r
n+1
2
exp
−r
2
2
2(1 − ρ
2
)
dr
2
×
ρ
(1 − ρ
2
)
n
π
−π
π
−π
senθ
1
cos θ
2
cos
n
(θ
1
− θ
2
)dθ
1
dθ
2
I
1
,
onde C = 1/[(2π)
2
(1 − ρ
2
)].
- Cálculo da integral I
1
Fazendo a seguinte mudança de variáveis
θ
1
= x + y
θ
2
= y
obtemos
I
1
=
S
sen(x + y) cos y cos
n
x dx dy =
S
(senx cos y + seny cos x) cos y cos
n
x dx dy
=
S
cos
n
x senx cos
2
y dx dy
I
11
+
S
cos
n+1
x seny cos y dx dy
I
12
onde S é a nova região de integração devido à mudança de variáveis.
- Cálculo da integral I
11
I
11
=
0
−2π
π
−x−π
cos
n
x senx cos
2
y dy dx +
2π
0
−x+π
−π
cos
n
x senx cos
2
y dy dx
=
0
−2π
cos
n
x senx
y
2
+
sen2y
4
π
−x−π
dx +
2π
0
cos
n
x senx
y
2
+
sen2y
4
−x+π
−π
dx
=
0
−2π
cos
n
x senx
π +
x
2
+
sen2x
4
dx +
2π
0
cos
n
x senx
π −
x
2
−
sen2x
4
dx
= π
0
−2π
cos
n
x senx dx +
1
2
0
−2π
x cos
n
x senx dx +
1
4
0
−2π
cos
n
x senx sen2x dx
+ π
2π
0
cos
n
x senx dx −
1
2
2π
0
x cos
n
x senx dx −
1
4
2π
0
cos
n
x senx sen2x dx
65
como as funções (x cos
n
x senx) e (cos
n
x senx sen2x) são pares, segue que
I
11
= π
0
−2π
cos
n
x senx dx + π
2π
0
cos
n
x senx dx
= π
0
−2π
cos
n
x senx dx +
2π
0
cos
n
x senx dx
= π
2π
−2π
cos
n
x senx dx = π.0 = 0 ,
já que a função (cos
n
x senx) é ímpar.
- Cálculo da integral I
12
I
12
=
0
−2π
π
−x−π
cos
n+1
x seny cos y dy dx +
2π
0
−x+π
−π
cos
n+1
x seny cos y dy dx
=
0
−2π
cos
n+1
x
sen
2
y
2
π
−x−π
dx +
2π
0
cos
n+1
x
sen
2
y
2
−x+π
−π
dx
= −
1
2
0
−2π
cos
n+1
x sen
2
xdx +
1
2
2π
0
cos
n+1
x sen
2
xdx = 0 ,
pois a função (cos
n+1
x sen
2
x) é uma função par. Portanto,
I
1
= 0, ∀n ,
donde segue que
E(senΘ
1
cos Θ
2
) = 0 .
• Cálculo de E(cos Θ
1
cos Θ
2
)
Por definição de valor esperado, temos que
E(cos Θ
1
cos Θ
2
) =
π
−π
π
−π
cos θ
1
cos θ
2
f
(Θ
1
,Θ
2
)
(θ
1
, θ
2
)dθ
1
dθ
2
=
π
−π
π
−π
cos θ
1
cos θ
2
∞
0
∞
0
f
(R
1
,Θ
1
,R
2
,Θ
2
)
(r
1
, θ
1
, r
2
, θ
2
)dr
1
dr
2
dθ
1
dθ
2
= C
∞
n=0
1
n!
∞
0
r
n+1
1
exp
−r
2
1
2(1 − ρ
2
)
dr
1
I
2
∞
0
r
n+1
2
exp
−r
2
2
2(1 − ρ
2
)
dr
2
I
3
×
ρ
(1 − ρ
2
)
n
π
−π
π
−π
cos θ
1
cos θ
2
cos
n
(θ
1
− θ
2
)dθ
1
dθ
2
I
4
.
66
- Cálculo da integral I
k
, k = 2, 3
Fazendo a substituição t = r
2
k
/[2(1 − ρ
2
)], segue que
r
k
= [2(1 − ρ
2
)]
1/2
t
1/2
⇒ dr
k
=
[2(1 − ρ
2
)]
1/2
2t
1/2
dt .
Daí,
I
k
=
∞
0
[2(1 − ρ
2
)]
1/2
t
1/2
n+1
exp(−t)
[2(1 − ρ
2
)]
1/2
2t
1/2
dt
=
[2(1 − ρ
2
)]
2
n+2
2
∞
0
t
(
n
2
+1
)
−1
exp(−t) dt
=
[2(1 − ρ
2
)]
2
n+2
2
Γ
n + 2
2
, k = 2, 3.
Ou seja,
I
2
= I
3
=
[2(1 − ρ
2
)]
2
n+2
2
Γ
n + 2
2
- Cálculo da integral I
4
Fazendo a seguinte mudança de variáveis
θ
1
= x + y
θ
2
= y
obtemos
I
4
=
S
cos(x + y) cos y cos
n
x dx dy =
S
(cos x cos y + senx seny) cos y cos
n
x dx dy
=
S
cos
n+1
x cos
2
y dx dy
I
41
−
S
cos
n+1
x cos
2
y dx dy
I
42
.
onde S é a nova região de integração devido à mudança de variáveis.
- Cálculo da integral I
42
I
42
=
0
−2π
π
−x−π
cos
n
x senx seny cos y dy dx +
2π
0
−x+π
−π
cos
n
x senx seny cos y dy dx
=
0
−2π
cos
n
x senx
sen
2
y
2
π
−x−π
dx +
2π
0
cos
n
x senx
sen
2
y
2
−x+π
−π
dx
67
= −
1
2
0
−2π
cos
n
x sen
3
xdx +
1
2
2π
0
cos
n
x sen
3
xdx
=
1
2
2π
0
cos
n
x sen
3
xdx +
1
2
2π
0
cos
n
x sen
3
xdx =
2π
0
cos
n
x sen
3
xdx
=
−
sen
2
x cos
n+1
x
n + 3
2π
0
+
2
n + 3
2π
0
cos
n
x senxdx =
2
n + 3
−
cos
n+1
x
n + 1
2π
0
= 0.
- Cálculo da integral I
41
I
41
=
0
−2π
π
−x−π
cos
n+1
x cos
2
y dy dx +
2π
0
−x+π
−π
cos
n+1
x cos
2
y dy dx
=
0
−2π
cos
n+1
x
y
2
+
sen2y
4
π
−x−π
dx +
2π
0
cos
n+1
x
y
2
+
sen2y
4
−x+π
−π
dx
=
0
−2π
cos
n+1
x
π +
x
2
+
sen2x
4
dx +
2π
0
cos
n+1
x
π −
x
2
−
sen2x
4
dx
= π
0
−2π
cos
n+1
x dx +
1
2
0
−2π
x cos
n+1
x dx +
1
2
0
−2π
cos
n+2
x senx dx
+ π
2π
0
cos
n+1
x dx −
1
2
2π
0
x cos
n+1
x dx −
1
2
2π
0
cos
n+2
x senx dx
= 2π
0
−2π
cos
n+1
x dx +
0
−2π
x cos
n+1
x dx +
0
−2π
cos
n+2
x senx dx
=
0
−2π
(x + 2π) cos
n+1
x dx +
−
cos
n+3
x
n + 3
0
−2π
=
0
−2π
(x + 2π) cos
n+1
x dx .
Fazendo
u = (x + 2π) ⇒ du = dx ,
dv = cos
n+1
x dx ⇒ v =
cos
n+1
x dx ,
Para calcular v utilizamos, recursivamente, a propriedade 2(ii) enunciada no início
deste apêndice. Observe que a aplicação recursiva dessa propriedade, reduz a potência
(n + 1) do cosseno e, dependendo da pa ridade de n, essa redução resulta em
cos x dx
ou
dx.
Se n for par então
v =
cos
n+1
x dx =
cos
n
x senx
(n + 1)
+
n
(n + 1)
cos
n−1
x dx
=
cos
n
x senx
(n + 1)
+
n
(n + 1)
cos
n−2
x senx
(n − 1)
+
n(n − 2)
(n + 1)(n − 1)
cos
n−3
x dx
68
=
cos
n
x senx
(n + 1)
+
n
(n + 1)
cos
n−2
x senx
(n − 1)
+
n(n − 2)
(n + 1)(n − 1)
cos
n−4
x senx
(n − 3)
+
n(n − 2)(n − 4)
(n + 1)(n − 1)(n − 3)
cos
n−5
x dx
=
cos
n
x senx
(n + 1)
+
n
(n + 1)
cos
n−2
x senx
(n − 1)
+
n(n − 2)
(n + 1)(n − 1)
cos
n−4
x senx
(n − 3)
+ . . . +
n(n − 2)(n − 4) . . . 8.6.4
(n + 1)(n − 1)(n − 3) . . . 9.7.5
cos
2
x senx
3
+
n(n − 2)(n − 4) . . . 8.6.4.2
(n + 1)(n − 1)(n − 3) . . . 9.7.5.3
cos xdx .
Integrando por partes, temos que
I
41
= [uv]
0
−2π
−
0
−2π
vdu = 0 −
0
−2π
vdu =
1
n
cos
n+1
x
(n + 1)
= +
n
(n + 1)(n − 1)
cos
n−1
x
(n − 1)
0
−2π
+ . . . +
n(n − 2)(n − 4) . . . 8.6.4
(n + 1)(n − 1)(n − 3) . . . 9.7.5.3
cos
3
x
3
0
−2π
+
n(n − 2)(n − 4) . . . 8.6.4.2
(n + 1)(n − 1)(n − 3) . . . 9.7.5.3
[cos x]
0
−2π
= 0 .
Portanto, se n for um número par entã o I
41
= 0.
Se n for ímpar então o cálculo de v através da aplicação recursiva da propriedade
2(ii) resulta em
v =
cos
n+1
x dx =
cos
n
x senx
(n + 1)
+
n
(n + 1)
cos
n−2
x senx
(n − 1)
= +
n(n − 2)
(n + 1)(n − 1)
cos
n−4
x senx
(n − 3)
+ . . . +
n(n − 2)(n − 4) . . . 7.5.3
(n + 1)(n − 1)(n − 3) . . . 8.6.4
cos x senx
2
+
n(n − 2)(n − 4) . . . 7.5.3.1
(n + 1)(n − 1)(n − 3) . . . 8.6.4.2
dx .
Daí, para n ímpar, temos que
I
41
= [uv]
0
−2π
−
0
−2π
vdu = 0 −
0
−2π
vdu
1
(n + 1)
cos
n+1
x
(n + 1)
0
−2π
= +
n
(n + 1)(n − 1)
cos
n−1
x
(n − 1)
0
−2π
+ . . . +
n(n − 2)(n − 4) . . . 7.5.3
(n + 1)(n − 1)(n − 3) . . . 8.6.4.2
cos
3
x
3
0
−2π
+
n(n − 2)(n − 4) . . . 7.5.3.1
(n + 1)(n − 1)(n − 3) . . . 8.6.4.2
−
x
2
0
−2π
=
n(n − 2)(n − 4) . . . 7.5.3.1
(n + 1)(n − 1)(n − 3) . . . 8.6.4.2
(2π
2
)
=
n(n − 1)(n − 2)(n − 3)(n − 4) . . . 8.7.6.5.4.3.2.1
(n + 1)[(n − 1)(n − 3) . . . 8.6.4.2]
2
(2π
2
)
=
n!
(n + 1)[(n − 1)(n − 3) . . . 8.6.4.2]
2
(2π
2
).
69
Sendo n ímpar, podemos escrevê-lo na forma n = 2m + 1, assim,
I
41
=
(2m + 1)!
(2m + 2)[(2m)(2m − 2) . . . 8.6.4.2]
2
(2π
2
)
=
(2m + 1)!
2(m + 1)[2
m
m(m − 1) . . . 4.3.2.1]
2
(2π
2
) =
(2m + 1)!
(m + 1)2
2m
[m!]
2
π
2
=
(2m + 1)!
2
2m
(m + 1)m!m!
π
2
=
(2m + 1)!
2
2m
(m + 1)!m!
π
2
=
π
2
2
2m
2m + 1
m
.
Logo,
I
41
=
0 se n for par;
π
2
2
2m
2m+1
m
se n for ímpar, n = 2m + 1.
Portanto,
I
4
= I
41
− I
42
= I
41
− 0 = I
41
=
0 se n for par;
π
2
2
2m
2m+1
m
se n for ímpar, n = 2m + 1.
Substituindo I
2
, I
3
e I
4
na expressão de E(cos Θ
1
cos Θ
2
), obtemos
E(cos Θ
1
cos Θ
2
) =
1
(2π)
2
(1 − ρ
2
)
∞
m=0
1
(2m + 1)!
1
4
2(1 − ρ
2
)
(2m+1+2)
×Γ
2
2m + 1 + 2
2
ρ
2m+1
(1 − ρ
2
)
2m+1
π
2
2
2m
2m + 1
m
=
1
2(1 − ρ
2
)
∞
m=0
2
2m
(1 − ρ
2
)
2m+3
Γ
2
m +
3
2
×
ρ
2m
ρ
(1 − ρ
2
)
2m+1
1
2
2m
1
m!(m + 1)!
=
1
2
ρ (1 − ρ
2
)
∞
m=0
Γ
2
m +
3
2
Γ (m + 2)
(ρ
2
)
m
m!
.
Podemos reescrever Γ(m + 2) da seguinte maneira
Γ(m + 2) = (m + 1)! = (m + 1)m(m − 1) . . . 4.3.2 = 2.3.4 . . . (m − 1)m(m + 1)
= (2 + 0)(2 + 1)(2 + 2) . . . (2 + m − 3)(2 + m − 2)(2 + m − 1)
=
m−1
i=0
(2 + i) = 2
[m]
.
Além disso, pode-se mostrar, por indução, que
Γ
m +
3
2
=
3
2
[m]
√
π
2
.
Assim,
70
E(cos Θ
1
cos Θ
2
) =
π
8
ρ(1 − ρ
2
)
∞
m=0
3
2
[m]
3
2
[m]
2
[m]
(ρ
2
)
m
m!
=
π
8
ρ(1 − ρ
2
)
2
F
1
3
2
,
3
2
, 2; ρ
2
,
onde
2
F
1
3
2
,
3
2
, 2; ρ
2
é a função hipergeométrica definida no início deste apêndice.
• Cálculo de E(senΘ
1
senΘ
2
)
O cálculo de E(senΘ
1
senΘ
2
) é feito de maneira análoga ao anterior, de onde vem
E(senΘ
1
senΘ
2
) =
π
8
ρ(1 − ρ
2
)
2
F
1
3
2
,
3
2
, 2; ρ
2
.
Substituindo, cada valor esperado calculado, na expressão do coeficiente de correlação,
obtemos
ρ
τ
=
π
2
16
ρ
2
(1 − ρ
2
)
2
2
F
1
3
2
,
3
2
, 2; ρ
2
2
,
como queríamos mostrar.
Referências Bibliográficas
[Box e Jenkins, 1976] Box, G. E. P. e Jenkins, G. M (1976). Time Series Analysis -
forecasting and control. San Francisco: Holden-Day.
[Charles et al., 2001] Charles, T., Guéméné, J. M., Vicent, G., Dubrule, O. e TotalFi-
naElf (2001). Experience with the Quantification of Subsurface Uncertainties.
Society of Petroleum Engineers, n
o
68703.
[Fisher, 1993] Fisher, N. I (1993). Statistical Analysis of Circular Data. Cambridge:
University Press.
[Fisher e Lee, 1983] Fisher, N. I. e Lee, A. J (1983). A correlation coefficient for circular
data. Biometrics 70, 327-332.
[Fisher e Lee, 1992] Fisher, N. I. e Lee, A. J (1992). Regression models for an angular
response. Biometrics 48, 665-677.
[Fisher e Lee, 1994] Fisher, N. I. e Lee, A. J (1994). Time series analysis of circular
data. Journal of the Royal Statistical Society, B 56, 327-339.
[James, 2002] James, B. B (2002). Probabilidade: Um curso em Nível Intermediário.
Rio de Janeiro: IMPA.
[Johnson e Kotz, 1969] Johnson, L. N. e Kotz, S (1969). Discrete Distributions. Boston.
[Lima] Lima, H. R. P. L.Fundamentos de Perfuração. Notas de um curso do Centro de
Desenvolvimento de Recursos Humanos da Petrobras.
[Mardia e Jupp, 2000] Mardia, K. V. e Jupp, P. E (2000). Directional Statistics. New
York: John Wiley & Sons.
72
[Morettin e Toloi, 2004] Morettin, P. A. e Toloi, C. M. C (2004). Análise de Séries
Temporais. São Paulo: Edgard Blücher.
[Silva, 2004] Silva, A. M (2004). Um Modelo Estocástico para Previsão de Desvios da
Coluna de Perfuração em Poços Petrolíferos. Monografia de Graduação - Pro-
grama de Recursos Humanos da ANP/MCT - PRH(25). Universidade Federal
de Campina Grande - UFCG. Campina Grande - PB.
[Thomas, 2001] Thomas, J. E (2001). Fundamentos de Engenharia de Petróleo. Editora
Interciência.
[Thomas, 2003] Thomas, G. B (2003). Cálculo. São Paulo: Addison Wesley.
[Wilks, 1995] Wilks, D. S (1995). Statistical Methods in the Atmosferic Sciences. Aca-
demic Press.
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
