Análise de Séries Temporais, Notas de estudo de Estatística
ricardo-petterle-1
ricardo-petterle-1

Análise de Séries Temporais, Notas de estudo de Estatística

115 páginas
50Números de download
1000+Número de visitas
58%de 0 votosNúmero de votos
Descrição
Estatística
100 pontos
Pontos de download necessários para baixar
este documento
Baixar o documento
Pré-visualização3 páginas / 115
Esta é apenas uma pré-visualização
3 mostrados em 115 páginas
Esta é apenas uma pré-visualização
3 mostrados em 115 páginas
Esta é apenas uma pré-visualização
3 mostrados em 115 páginas
Esta é apenas uma pré-visualização
3 mostrados em 115 páginas

ANÁLISE DE SÉRIES TEMPORAIS

RICARDO S. EHLERS

Primeira publicação 2003

Segunda edição publicada em 2004

Terceira edição publicada em 2005

Quarta edição publicada em 2007© RICARDO SANDES EHLERS 2003-2007

Sumário

1 Introdução 1

2 Técnicas Descritivas 6

2.1 Decomposição Clássica . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2 Séries com Tendência . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3 Séries Sazonais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.4 Autocorrelação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.4.1 O Correlograma . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Modelos Probabiĺısticos 18

3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.2 Processos Estacionários . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 A Função de Autocorrelação . . . . . . . . . . . . . . . . . . . . . . . . 20

3.4 Alguns Processos Estocásticos . . . . . . . . . . . . . . . . . . . . . . . 20

3.4.1 Sequência Aleatória . . . . . . . . . . . . . . . . . . . . . . . . 20

3.4.2 Passeio Aleatório . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.4.3 Processos de Média Móveis . . . . . . . . . . . . . . . . . . . . 22

3.4.4 Processos Autoregressivos . . . . . . . . . . . . . . . . . . . . . 24

3.4.5 Modelos Mistos ARMA . . . . . . . . . . . . . . . . . . . . . . 29

3.4.6 Modelos ARMA Integrados . . . . . . . . . . . . . . . . . . . . 30

4 Estimação 33

4.1 Autocovariância e autocorrelação . . . . . . . . . . . . . . . . . . . . . 34

4.2 Ajustando Processos Autoregressivos . . . . . . . . . . . . . . . . . . . 35

4.3 Ajustando Processos Médias Móveis . . . . . . . . . . . . . . . . . . . 39

4.4 Ajustando Processos ARMA . . . . . . . . . . . . . . . . . . . . . . . . 40

4.5 Modelos Sazonais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.6 Adequação do Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.6.1 Análise dos Reśıduos . . . . . . . . . . . . . . . . . . . . . . . . 43

4.6.2 Testes sobre os reśıduos . . . . . . . . . . . . . . . . . . . . . . 44

5 Previsão 52

5.1 Métodos Univariados de Previsão . . . . . . . . . . . . . . . . . . . . . 52

5.1.1 Alisamento Exponencial Simples . . . . . . . . . . . . . . . . . 52

i

ii SUMÁRIO

5.1.2 Método de Holt-Winters . . . . . . . . . . . . . . . . . . . . . . 56

5.2 Previsão em Modelos ARMA . . . . . . . . . . . . . . . . . . . . . . . 58

5.3 Performance Preditiva . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.4 Critérios de Informação . . . . . . . . . . . . . . . . . . . . . . . . . . 64

5.5 Previsões Usando Todos os Modelos . . . . . . . . . . . . . . . . . . . 67

5.6 Previsão Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

6 Modelando a Variância 73

6.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

6.2 Modelos ARCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.3 Modelos GARCH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

6.3.1 Estimação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.3.2 Adequação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

6.4 Volatilidade Estocástica . . . . . . . . . . . . . . . . . . . . . . . . . . 83

7 Modelos Lineares Dinâmicos 86

7.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

7.2 Modelos Polinomiais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

7.2.1 Análise Sequencial e Previsões . . . . . . . . . . . . . . . . . . 89

7.2.2 Variâncias de Evolução e das Observações . . . . . . . . . . . . 91

7.3 Modelo de Crescimento Linear . . . . . . . . . . . . . . . . . . . . . . 94

7.4 Modelos Sazonais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

7.4.1 Modelos sem Crescimento . . . . . . . . . . . . . . . . . . . . . 95

7.4.2 Modelos com Crescimento . . . . . . . . . . . . . . . . . . . . . 96

7.5 Representação de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . 96

7.6 Ilustração . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.7 Modelos de Regressão . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

7.8 Monitoramento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

A Lista de Distribuições 105

A.1 Distribuição Normal . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

A.2 Distribuição Gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

A.3 Distribuição Wishart . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

A.4 Distribuição Gama Inversa . . . . . . . . . . . . . . . . . . . . . . . . . 106

A.5 Distribuição Wishart Invertida . . . . . . . . . . . . . . . . . . . . . . 106

A.6 Distribuição Beta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

A.7 Distribuição de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . 107

A.8 Distribuição t de Student . . . . . . . . . . . . . . . . . . . . . . . . . 107

A.9 Distribuição F de Fisher . . . . . . . . . . . . . . . . . . . . . . . . . . 107

A.10 Distribuição Binomial . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

A.11 Distribuição Multinomial . . . . . . . . . . . . . . . . . . . . . . . . . 108

A.12 Distribuição de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . 108

A.13 Distribuição Binomial Negativa . . . . . . . . . . . . . . . . . . . . . . 108

SUMÁRIO iii

References 110

Caṕıtulo 1

Introdução

Uma série temporal é uma coleção de observações feitas sequencialmente ao longo do

tempo. A caracteŕıstica mais importante deste tipo de dados é que as observações vizi-

nhas são dependentes e estamos interessados em analisar e modelar esta dependência.

Enquanto em modelos de regressão por exemplo a ordem das observações é irrelevante

para a análise, em séries temporais a ordem dos dados é crucial. Vale notar também

que o tempo pode ser substituido por outra variável como espaço, profundidade, etc.

Como a maior parte dos procedimentos estat́ısticos foi desenvolvida para analisar

observações independentes o estudo de séries temporais requer o uso de técnicas espe-

ćıficas. Dados de séries temporais surgem em vários campos do conhecimento como

Economia (preços diários de ações, taxa mensal de desemprego, produção industrial),

Medicina (eletrocardiograma, eletroencefalograma), Epidemiologia (número mensal

de novos casos de meningite), Meteorologia (precipitação pluviométrica, temperatura

diária, velocidade do vento), etc.

Algumas caracteŕısticas são particulares a este tipo de dados, por exemplo,ˆ Observações correlacionadas são mais dif́ıceis de analisar e requerem técnicas espećıficas.ˆ Precisamos levar em conta a ordem temporal das observações.ˆ Fatores complicadores como presença de tendências e variação sazonal ou ćıclica podem ser dif́ıceis de estimar ou remover.ˆ A seleção de modelos pode ser bastante complicada, e as ferramentas podem ser de dif́ıcil interpretação.ˆ É mais dif́ıcil de lidar com observações perdidas e dados discrepantes devido à natureza sequencial.

Terminologia

Uma série temporal é dita ser cont́ınua quando as observações são feitas continuamente

no tempo. Definindo o conjunto T = {t : t1 < t < t2} a série temporal será denotada

1

2 CAPÍTULO 1. INTRODUÇÃO

por {X(t) : t ∈ T}. Uma série temporal é dita ser discreta quando as observações são feitas em tempos espećıficos, geralmente equiespaçados. Definindo o conjunto

T = {t1, . . . , tn} a série temporal será denotada por {Xt : t ∈ T}. Por simplicidade podemos fazer T = {1, 2, . . . , n}.

Note que estes termos não se referem à variável observada X, esta pode assumir

valores discretos ou cont́ınuos. Em muitas situações X pode ser discreta por definição

(e.g. o número de casos notificados de AIDS) porém para efeito de análise estat́ıs-

tica pode ser tratada como continua se os seus valores observados não forem muito

pequenos.

Por outro lado, séries temporais discretas podem surgir de várias formas. Séries

cont́ınuas podem ser discretizadas, i.e. seus valores são registrados a certos intervalos

de tempo. Séries de valores agregados ou acumulados em intervalos de tempo, por

exemplo exportações medidas mensalmente ou quantidade de chuva medida diaria-

mente. Finalmente, algumas séries são inerentemente discretas, por exemplo dividen-

dos pagos por uma empresa aos seus acionistas em anos sucessivos.

Uma série temporal também pode ser multivariada. Se k variáveis são observadas

a cada tempo (por exemplo discreto) denota-se por {X1t, . . . , Xkt, t ∈ T}. Neste caso várias séries correlacionadas devem ser analisadas conjuntamente, ou seja em cada

tempo tem-se um vetor de observações.

Objetivos

Em algumas situações o objetivo pode ser fazer previsões de valores futuros enquanto

em outras a estrutura da série ou sua relação com outras séries pode ser o interesse

principal. De um modo geral, os principais objetivos em se estudar séries temporais

podem ser os seguintes,ˆ Descrição. Descrever propriedades da série, e.g. o padrão de tendência, existên- cia de variação sazonal ou ćıclica, observações discrepantes (outliers), alterações

estruturais (e.g. mudanças no padrão da tendência ou da sazonalidade), etc.ˆ Explicação. Usar a variação em uma série para explicar a variação em outra série.ˆ Predição: predizer valores futuros com base em valores passados. Aqui assume- se que o futuro envolve incerteza, ou seja as previsões não são perfeitas. Porém

devemos tentar reduzir os erros de previsão.ˆ Controle. Os valores da série temporal medem a “qualidade” de um processo de manufatura e o objetivo é o controle do processo. Um exemplo é o controle

estat́ıstico de qualidade aonde as observações são representadas em cartas de

controle. Este tópico não será abordado nestas notas de aula.

3

Abordagensˆ Técnicas Descritivas. Técnicas gráficos, identificação de padrões, etc.ˆ Modelos Probabiĺısticos. Seleção, comparação e adequação de modelos, estima- ção, predição. Ferramenta básica é a função de autocorrelação.ˆ Análise espectral.ˆ Métodos não paramétricos (alisamento ou suavização).ˆ Outras Abordagens. Modelos de espaço de estados, modelos não lineares, séries multivariadas, estudos longitudinais, processos de longa dependência, modelos

para volatilidade, etc.

Sazonalidade

Muitas séries temporais exibem um comportamento que tende a se repetir a cada s

peŕıodos de tempo. Por exemplo, é natural esperar que as vendas mensais de brinque-

dos terão um pico no mês de dezembro e talvez um pico secundário em outubro. Este

padrão possivelmente se repetirá ao longo de vários anos. Vejamos alguns posśıveis

modelos sazonais,

1. Sazonalidade deterministica. Variáveis dummies (binárias). O coeficiente de

cada variável dummy representa o fator sazonal do respectivo mês, trimestre,

etc.

2. Funções trigonométricas.

3. Sazonalidade estocástica:

(a) Variável endógena com defasagem sazonal no modelo (modelos ARMA

periódicos),

(b) modelo ARMA sazonal.

Tipos de Sazonalidadeˆ Aditiva. A série apresenta flutuações sazonais mais ou menos constantes não importando o ńıvel global da série.ˆ Multiplicativa. O tamanho das flutuações sazonais varia dependendo do ńıvel global da série.

No exemplo dos brinquedos, suponha que o aumento esperado nas vendas nos

meses de dezembro é de 1 milhão de reais em relação à média anual. Então as

previsões para os meses de dezembro dos próximos anos deve somar a quantia de 1

milhão de reais à uma média anual para levar em conta esta flutuação sazonal. Isto

é o que se chama de sazonalidade aditiva.

4 CAPÍTULO 1. INTRODUÇÃO

Suponha agora que o aumento esperado nos meses de dezembro seja de 30%.

Então o aumento esperado (em valor absoluto) de vendas em dezembro será pequeno

ou grande dependendo da média anual de vendas ser baixa ou alta. Nas previsões

para os próximos meses de dezembro deve-se multiplicar a média anual pelo fator 1,3.

Isto é o que se chama de sazonalidade multiplicativa.

Tendência

Globalmente, uma série pode exibir tendência de crescimento (ou decrescimento) com

vários posśıveis padrões.ˆ Crescimento linear. Por exemplo, a cada ano o aumento esperado nas vendas de um certo brinquedo é de 1 milhão de reais.ˆ Crescimento exponencial. Por exemplo, a cada ano as vendas de um certo brinquedo aumentam de um fator 1,3.ˆ Crescimento amortecido. Por exemplo, as vendas de um certo brinquedo tem uma aumento esperado de 70% sobre o ano anterior. Se o aumento esperado for

de 1 milhão de reais no primeiro ano, no segundo ano será de 700 mil reais, no

terceiro ano será de 490 mil reais e assim por diante.

Exemplos de Séries Temporais

Como primeira ilustração são apresentadas na Figura 1.1 quatro séries temporais

dispońıveis no pacote R. Nos eixos horizontais aparecem os anos de observação e nos

eixos verticais os nomes das séries (mesmos nomes do R). A Figura 1.1a mostra totais

mensais de passageiros em linhas aéreas internacionais nos EUA entre 1949 e 1960.

Existe uma clara tendência de crescimento bem como um padrão sazonal ao longo

dos anos. A Figura 1.1b mostra a série com o número anual de linces capturados em

armadilhas entre 1821 e 1934 no Canadá. Existe um padrão ćıclico em torno de 10 ou

11 anos. A Figura 1.1c mostra a série com as medições anuais de vazões do Rio Nilo

em Ashwan entre 1871 e 1970. Parece haver alguma alteração estrutural em torno do

ano de 1900. Finalmente a Figura 1.1d mostra a série trimestral do consumo de gás no

Reino Unido entre o primeiro trimestre de 1960 e o quarto trimestre de 1986. Há uma

tendência de crescimento porém a amplitude do padrão sazonal aumenta bastante a

partir de 1971.

Exerćıcios

1. Classifique as seguintes séries temporais quanto ao tempo e quanto a variável

observada.

(a) Registros de maré durante 1 dia.

(b) Medidas de temperatura em uma estação meteorológica.

5

A irP

as se

ng er

s

1950 1954 1958

10 0

30 0

50 0

ly nx

1820 1860 1900

0 20

00 50

00

N ile

1880 1920 1960

60 0

10 00

14 00

U K

ga s

1960 1970 1980

20 0

60 0

10 00

Figura 1.1: (a) Totais mensais de passageiros em linhas aéreas internacionais nos EUA

entre 1949 e 1960, (b) número anual de linces capturados em armadilhas entre 1821

e 1934 no Canadá, (c) medições anuais de vazões do Rio Nilo em Ashwan entre 1871

e 1970, (d) consumo de gás no Reino Unido entre o primeiro trimestre de 1960 e o

quarto trimestre de 1986.

(c) O ı́ndice diário da bolsa de valores de São Paulo.

(d) A inflação mensal medida pelo ı́ndice de preços ao consumidor.

(e) Variação diária de um determinado ı́ndice financeiro, 1 para variação po-

sitiva, -1 para variação negativa ou zero se não ocorreu variação.

(f) Número mensal de novos casos de Dengue em uma determinada região.

2. Dê exemplos de séries temporais continuas que poderiam ser discretizadas (e de

que forma).

Caṕıtulo 2

Técnicas Descritivas

Ao se analisar uma ou mais séries temporais a representação gráfica dos dados se-

quencialmente ao longo do tempo é fundamental e pode revelar padrões de comporta-

mento importantes. Tendências de crescimento (ou decrescimento), padrões ćıclicos,

alterações estruturais, observações aberrantes, etc. são muitas vezes facilmente iden-

tificados. Sendo assim, o gráfico temporal deve ser sempre o primeiro passo e antecede

qualquer análise. Outras ferramentas serão descritas ao longo deste caṕıtulo.

2.1 Decomposição Clássica

Muitas das propriedades observadas em uma série temporal Xt podem ser captadas

assumindo-se a seguinte forma de decomposição

Xt = Tt + Ct +Rt

onde Tt é uma componente de tendência, Ct é uma componente ćıclica ou sazonal e

Rt é uma componente aleatória ou rúıdo (a parte não explicada, que espera-se ser

puramente aleatória). A componente ćıclica se repete a cada intervalo fixo s, i.e.

· · · = Ct−2s = Ct−s = Ct = Ct+s = Ct+2s = . . . .

Assim, variações periódicas podem ser captadas por esta componente.

2.2 Séries com Tendência

Não existe uma definição precisa de tendência e diferentes autores usam este termo de

diferentes formas. Podemos pensar em tendência como uma mudança de longo prazo

no ńıvel médio da série. A dificuldade aqui é definir longo prazo.

A forma mais simples de tendência é

Xt = α+ βt+ ǫt (2.1)

onde α e β são constantes a serem estimadas e ǫt denota um erro aleatório com média

zero. O ńıvel médio da série no tempo t é dado por mt = α+ βt que é algumas vezes

6

2.2. SÉRIES COM TENDÊNCIA 7

chamado de termo de tendência. Porém alguns autores preferem chamar a inclinação

β de tendência, ou seja a mudança no ńıvel da série por unidade de tempo já que

β = mt −mt−1. Note que a tendência na equação (2.1) é uma função determińıstica do tempo e algumas vezes é chamada de tendência global (i.e. vale para toda a série),

em oposição a tendência local.

De um modo geral, uma forma de se lidar com dados não sazonais que contenham

uma tendência consiste em ajustar uma função polinomial,

Xt = β0 + β1t+ · · ·+ βp tp + ǫt.

Uma função linear ou quadrática seria apropriada no caso de uma tendência mono-

tonicamente crescente ou decrescente. Caso contrário polinômios de ordem mais alta

devem ser ajustados.

Outras posśıveis formas de tendência são os crescimentos descritos por uma curva

Gompertz,

log xt = a+ br t

onde a, b e r são parâmetros com 0 < r < 1, ou uma curva Loǵıstica,

xt = a/(1 + be −ct)

onde a, b e c são parâmetros. Estas duas últimas são chamadas curvas S e se apro-

ximam de uma asśıntota quando t → ∞. Neste caso o ajuste pode levar a equações não lineares.

Seja qual for a curva utilizada, a função ajustada fornece uma medida da tendência

da série, enquanto os reśıduos (valores observados – valores ajustados) fornecem uma

estimativa de flutuações locais.

Exemplo 2.1 : A Figura 2.1 mostra as medições anuais de vazões do Rio Nilo em

Ashwan entre 1871 e 1970 juntamente com polinômios de graus 3 e 6 superimpostos.

Os polinômios foram ajustados por mı́nimos quadrados usando os comandos do R a

seguir. A série original com as tendências estimadas aparecem na Figura (2.1).

> mypolytrend = function(y, degree = 1) {

+ n = length(y)

+ x = 1:n

+ X = matrix(NA, n, degree)

+ for (i in 1:degree) X[, i] = x^i

+ a = as.numeric(lm(y ~ X)$coeff)

+ out = cbind(rep(1, n), X) %*% a

+ return(ts(out, start = start(y), freq = frequency(y)))

+ }

> z3 = mypolytrend(Nile, 3)

> z6 = mypolytrend(Nile, 6)

8 CAPÍTULO 2. TÉCNICAS DESCRITIVAS

V az

oe s

1880 1900 1920 1940 1960

60 0

80 0

10 00

12 00

14 00

observado tendencia grau 3 tendencia grau 6

Figura 2.1: Medições anuais de vazões do Rio Nilo em Ashwan entre 1871 e 1970 (pontos), com polinômios de graus 3 e 6 ajustados por minimos quadrados.

Regressão Local

A idéia aqui é estimar para cada t uma equação de regressão polinomial diferente,

por exemplo

x̂t = α̂(t) + β̂(t)t.

Note que as estimativas de α e β dependem do tempo o que da o caráter local das

retas de regressão.

O procedimento conhecido como loess é um procedimento iterativo que a cada

passo aplica a regressão local anterior, calcula os reśıduos xt − x̂t e aplica novamente a regressão local dando peso menor às observações com reśısuos maiores. Este proce-

dimento se repete até atingir convergência.

Exemplo 2.2 : A Figura 2.2 apresenta os mesmos dados da Figura 2.1 sendo que as

curvas superimpostas foram obtidas usando regressão local com os comandos do R a

seguir.

2.2. SÉRIES COM TENDÊNCIA 9

V az

oe s

1880 1900 1920 1940 1960

60 0

80 0

10 00

12 00

14 00

observado tendencia f=1 tendencia f=0.25

Figura 2.2: Medições anuais de vazões do Rio Nilo em Ashwan entre 1871 e 1970 (pontos), tendência estimada via função lowess.

Filtragem

Outro procedimento para analisar séries com tendência é através de filtros lineares.

Um filtro linear converte uma série {xt} em outra {yt} através da seguinte operação linear

yt = s

j=−q

ajxt+j

onde {aj} é um conjunto de pesos. Além disso, como queremos estimar a média local os pesos devem ser tais que

∑s j=−q aj = 1, garantindo assim que min{xt} < yt <

max{xt}. Neste caso a operação é chamada média móvel. Médias móveis são em geral simétricas com s = q e a−r = ar. Por exemplo, se

s = q = 2 temos que

yt = a2xt−2 + a1xt−1 + a0xt + a1xt+1 + a2xt+2.

O caso mais simples é quando todos os pesos aj tem o mesmo valor e devido à restrição

de soma 1 segue que aj = 1/(2q+1), para j = −q, . . . , q. Neste caso, o valor suavizado

10 CAPÍTULO 2. TÉCNICAS DESCRITIVAS

de xt é dado por

yt = 1

2q + 1

q ∑

j=−q

xt+j .

Qualquer que seja o filtro utilizado, yt é uma estimativa da tendência no tempo t

e xt − yt é uma série livre de tendência.

Exemplo 2.3 : A Figura 2.3 apresenta a série com os totais mensais de passageiros

de linhas aéreas internacionais nos EUA, entre 1949 e 1960 (Box, Jenkins and Reinsel,

1976) juntamente com a tendência “estimada” superimposta. Foram aplicados filtros

lineares com médias móveis aproximadamente trimestrais (q = 2) e médias móveis

aproximadamente anuais (q = 5). Os seguintes comandos do R foram usados.

Anos

N um

er o

de p

as sa

ge iro

s (e

m m

ilh ar

es )

1950 1952 1954 1956 1958 1960

0 10

0 20

0 30

0 40

0 50

0 60

0

dados Media Movel q=2 Media Movel q=5

Figura 2.3: Totais mensais de passageiros de linhas aéreas internacionais nos EUA, com a tendência superimposta aplicando médias móveis aproximadamente trimestrais (q = 2) e

médias móveis aproximadamente anuais (q = 5).

Note que, para a aplicação de qualquer filtro simétrico os valores suavizados só

podem ser calculados para t = q + 1, . . . , n− q e assim a série suavizada terá n− 2q valores. Em algumas situações no entanto é importante obter valores suavizados até

o peŕıodo t = n e uma alternativa é utilizar um filtro assimétrico que usa apenas os

valores atual e passados de xt. Por exemplo na técnica conhecida como alisamento

2.3. SÉRIES SAZONAIS 11

exponencial os valores suavizados são dados por

yt = ∞ ∑

j=0

α(1− α)jxt−j

onde 0 < α < 1. Note como, embora todas as observações passadas sejam usadas no

filtro, os pesos α(1−α)j decaem geometricamente com j. Quanto mais próximo de 1 estiver α mais peso será dado às observações mais recentes e quanto mais próximo de

zero mais os pesos estarão distribuidos ao longo da série. Por exemplo se α = 0, 90 a

série filtrada fica yt = 0, 9xt+0, 09xt−1 +0, 009xt−2 + . . . enquanto que para α = 0, 1

temos que yt = 0, 1xt + 0, 09xt−1 + 0, 081xt−2 + . . . .

Este tipo de filtro pode ser utilizado para fazer previsões. Especificamente a

previsão da série original em t+ 1 será o valor filtrado yt (mais detalhes no Caṕıtulo

5).

Diferenciação

Um tipo especial de filtro, muito útil para remover uma componente de tendência

polinomial, consiste em diferenciar a série até que ela se torne estacionária (este con-

ceito será formalizado no Caṕıtulo 3). Para dados não sazonais, a primeira diferença é

em geral suficiente para induzir estacionariedade aproximada. A nova série y2, . . . , yn é formada a partir da série original x1, . . . , xn como

yt = xt − xt−1 = ▽xt.

Note que isto nada mais é do que um filtro (assimétrico) com coeficientes 1 e -1.

Diferenciação de primeira ordem é a mais utilizada sendo que ocasionalmente uma

diferenciação de segunda ordem pode ser requerida, i.e.

yt = ▽ 2xt = ▽(xt − xt−1) = xt − 2xt−1 + xt−2.

Além disso, independente do seu uso para induzir estacionariedade, a diferencia-

ção pode ser muito útil como ferramenta exploratória. Observações discrepantes por

exemplo podem ter um efeito dramático na série diferenciada e uma representação

gráfica é em geral suficiente para identificar tais pontos.

2.3 Séries Sazonais

Uma forma bastante simples de eliminar o efeito sazonal é simplesmente tomar médias

sazonais. Por exemplo, em dados mensais com sazonalidade anual, as médias anuais

estarão livres do efeito sazonal. Embora este procedimento esteja correto muitos dados

serão perdidos e ao invés disto pode-se recorrer mais uma vez às médias móveis.

2.4 Autocorrelação

Uma importante ferramenta para se identificar as propriedades de uma série tem-

poral consiste de uma série de quantidades chamadas coeficientes de autocorrelação

12 CAPÍTULO 2. TÉCNICAS DESCRITIVAS

amostral. A idéia é similar ao coeficiente de correlação usual, i.e. para n pares de

observações das variáveis x e y o coeficiente de correlação amostral é dado por

r =

n ∑

i=1

(xi − x)(yi − y) √

n ∑

i=1

(xi − x)2 n ∑

i=1

(yi − y)2 . (2.2)

Aqui no entanto queremos medir a correlação entre as observações de

uma mesma variável em diferentes horizontes de tempo, i.e. correlações en-

tre observações defasadas 1, 2, . . . peŕıodos de tempo. Assim, dadas n ob-

servações x1, . . . , xn de uma série temporal discreta podemos formar os pares

(x1, x2), . . . , (xn−1, xn). Considerando x1, . . . , xn−1 e x2, . . . , xn como duas variáveis

o coeficiente de correlação entre xt e xt+1 é dado por

r1 =

n−1 ∑

t=1

(xt − x1)(xt+1 − x2) √

n−1 ∑

t=1

(xt − x1)2 n−1 ∑

t=1

(xt+1 − x2)2 (2.3)

onde as médias amostrais são

x1 = n−1 ∑

t=1

xt/(n− 1) e x2 = n ∑

t=2

xt/(n− 1).

Como o coeficiente r1 mede as correlações entre observações sucessivas ele é chamado

de coeficiente de autocorrelação ou coeficiente de correlação serial.

É usual simplificar a equação (2.3) utilizando-se a média de todas as observações,

i.e. x = n ∑

t=1

xt/n já que x1 ≈ x2, e assumindo variância constante. Assim, a versão

simplificada de (2.3) fica

r1 =

n−1 ∑

t=1

(xt − x)(xt+1 − x)

(n− 1) n ∑

t=1

(xt − x)2/n (2.4)

sendo que alguns autores ainda retiram o termo n/(n − 1) que é próximo de 1 para n não muito pequeno. Esta última forma simplificada, sem o termo n/(n − 1) será utilizada neste texto.

A equação (2.4) pode ser generalizada para calcular a correlação entre observações

defasadas de k peŕıodos de tempo, i.e.

rk =

n−k ∑

t=1

(xt − x)(xt+k − x)

n ∑

t=1

(xt − x)2 (2.5)

2.4. AUTOCORRELAÇÃO 13

fornece o coeficiente de correlação de ordem k. Assim como o coeficiente de correlação

usual, as autocorrelações são adimensionais e −1 < rk < 1. Na prática é mais usual calcular primeiro os coeficientes de autocovariância {ck},

definidos por analogia com a fórmula usual de covariância, i.e.

ck = n−k ∑

t=1

(xt − x)(xt+k − x)/n.

Os coeficientes de autocorrelação são então obtidos como rk = ck/c0.

2.4.1 O Correlograma

Um gráfico com os k primeiros coeficientes de autocorrelação como função de k é

chamado de correlograma e pode ser uma ferramenta poderosa para identificar ca-

racteŕısticas da série temporal. Porém isto requer uma interpretação adequada do

correlograma, i.e. devemos associar certos padrões do correlograma como determina-

das caracteŕısticas de uma série temporal. Esta nem sempre é uma tarefa simples e a

seguir são dadas algumas indicações.

Séries aleatórias

A primeira questão que podemos tentar responder através do correlograma é se uma

série temporal é aleatória ou não. Para uma série completamente aleatória os valores

defasados são não correlacionados e portanto espera-se que rk ≈ 0, k = 1, 2, . . . . Suponha que x1, . . . , xn sejam variáveis aleatórias independentes e identicamente

distribuidas com média arbitrárias. Então, pode-se mostrar que o coeficiente de au-

tocorrelação amostral rk é assintoticamente normalmente distribuido, com média e

variância dados por

E(rk) ≈ −1/n e V ar(rk) ≈ 1/n.

(ver Kendall, Stuart, & Ord 1983, Caṕıtulo 48). Portanto, limites de confiança apro-

ximados de 95% são dados por −1/n± 1, 96/√n, que são frequentemente ainda mais aproximados para ±1, 96/√n.

Isto ilustra uma das dificuldades de interpretar o correlograma já que, mesmo para

uma série completamente aleatória, espera-se que 1 em cada 20 coeficientes rk esteja

fora destes limites. Por outro lado, um valor muito grande de rk tem menos chance de

ter ocorrido ao acaso do que um valor que está apenas ligeiramente fora dos limites.

A Figura 2.4 mostra uma série temporal com 100 observações independentes e

identicamente distribuidas geradas no computador juntamente com o seu correlo-

grama. Neste caso os limites de confiança de 95% são aproximadamente ±2/ √ 100 =

±0,2 e podemos notar que 2 dentre as 20 primeiras autocorrelações estão ligeiramente fora destes limites. No entanto isto ocorre em defasagens aparentemente arbitrárias

(12 e 18) e podemos concluir que não há evidência para rejeitar a hipótese de que as

observações são independentes.

14 CAPÍTULO 2. TÉCNICAS DESCRITIVAS

tempo

ob se

rv aç

õe s

0 20 40 60 80 100

− 3

− 1

1

0 5 10 15 20

− 0.

2 0.

4 1.

0

defasagem

au to

co rr

el aç

oe s

Figura 2.4: (a) 100 observações simuladas independentes e identicamente distribuidas. (b) 20 primeiras autocorrelações amostrais.

Correlação de curto-prazo

Uma série temporal na qual uma observação acima da média “tende” a ser seguida

por uma ou mais observações acima da média, similarmente para observações abaixo

da média, é dita ter correlação de curto-prazo. Um correlograma desta série deverá

exibir um valor relativamente grande de r1 seguido por valores que tendem a ficar

sucessivamente menores. A partir de uma certa defasagem k os valores de rk tendem

a ser aproximadamente zero. Na Figura 2.5 temos 50 observações geradas de acordo

com o processo xt = 0, 7xt−1 + ǫt juntamente com o seu correlograma.

Correlação negativa

Se os valores de uma série temporal tendem a se alternar acima e abaixo de um valor

médio, o correlograma desta série também tende a se alternar. O valor de r1 será

negativo enquanto o valor de r2 será positivo já que as observações defasadas de 2

peŕıodos tendem a estar do mesmo lado da média. Esta caracteŕıstica está ilustrada

na Figura 2.6 aonde temos 50 observações simuladas com autocorrelações negativas

juntamente com as 14 primeiras autocorrelações amostrais.

2.4. AUTOCORRELAÇÃO 15

tempo

ob se

rv ac

oe s

0 10 20 30 40 50

− 3

− 1

1

0 5 10 15

− 0.

4 0.

2 0.

8

defasagem

au to

co rr

el ac

oe s

Figura 2.5: (a) 50 observações simuladas com autocorrelações de curto-prazo. (b) 16 primei- ras autocorrelações amostrais.

Séries não estacionárias

Para uma série temporal com tendência os valores de rk não decairão para zero a

não ser em defasagens grandes. Intuitivamente, isto ocorre porque uma observação

de um lado da média tende a ser seguida por um grande número de observações do

mesmo lado (devido à tendência). Neste caso, pouca ou nenhuma informação pode

ser extraida do correlograma já que a tendência dominará outras caracteŕısticas. Na

verdade, como veremos em outros caṕıtulos a função de autocorrelação só tem um

significado para séries estacionárias, sendo assim qualquer tendência deve ser removida

antes do cálculo de {rk}. A Figura 2.7 mostra uma série temporal com 50 observações geradas segundo o

modelo xt = xt−1+ǫt, juntamente com o seu correlograma. Note que a não estaciona-

riedade da série fica evidenciada no correlograma já que as autocorrelações amostrais

decaem muito lentamente.

Variação sazonal

Um padrão sazonal é em geral facilmente identificado no correlograma. De fato, se

uma série temporal contem flutuações sazonais o correlograma irá exibir oscilações na

16 CAPÍTULO 2. TÉCNICAS DESCRITIVAS

tempo

ob se

rv ac

oe s

0 10 20 30 40 50

− 3

− 1

1

0 5 10 15

− 0.

5 0.

5

defasagem

au to

co rr

el ac

oe s

Figura 2.6: (a) 50 observações simuladas com autocorrelações negativas. (b) 15 primeiras autocorrelações amostrais.

mesma frequência. Por exemplo, com observações mensais r6 será grande e negativo

enquanto r12 será grande e positivo.

Na verdade, se o padrão sazonal já for evidente no gráfico da série original o

correlograma trará pouca ou nenhuma informação adicional.

Observações discrepantes

Se uma série temporal contem uma ou mais observações discrepantes (“outliers”) o

correlograma pode ser seriamente afetado. No caso de uma única observação discre-

pante o gráfico de xt contra xt+k terá pontos extremos o que pode viesar os coeficientes

de correlação para zero. Com dois valores discrepantes o efeito pode ser ainda mais

devastador, além de gerar uma correlação espúria quando k é igual à distância entre

os valores.

Exerćıcios

1. Use o R para gerar uma série temporal Yt = b0 + b1t + ǫt, t = 1, . . . , 100, com

b0, b1 6= 0 e ǫt normais e independentes com média µ e variância σ21 se t ≤ 70 mas variância σ22 6= σ21 se t > 70. Usando diferentes valores de α aplique o

2.4. AUTOCORRELAÇÃO 17

tempo

ob se

rv ac

oe s

0 10 20 30 40 50

0. 0

0. 4

0. 8

0 5 10 15 20

− 0.

2 0.

4 1.

0

defasagem

au to

co rr

el ac

oe s

Figura 2.7: (a) 50 observações simuladas segundo um passeio aleatório. (b) 20 primeiras autocorrelações amostrais.

alisamento exponencial e faça um gráfico da série com os valores suavizados.

Comente os resultados.

2. Para cada um dos processos abaixo gere 200 observações. Faça um gráfico da

série e do correlograma.

(a) Serie aleatória, observações iid da distribuição N(0,1).

(b) Serie com tendência estocástica, xt = xt−1 + ǫt, ǫt ∼ N(0, (0, 1)2) (c) Outra serie com tendencia estocastica, xt = xt−1 + ǫt, ǫt ∼ N(1, 52) (d) Serie com correlação de curto-prazo, xt = 0, 7xt−1 + ǫt, ǫt ∼ N(0, 1) (e) Serie com correlações negativas, xt = −0, 8xt−1 + ǫt, ǫt ∼ N(0, 1) (f) Medias moveis, xt = ǫt + 0, 6ǫt−1, ǫt ∼ N(0, 1) (g) passeio aleatorio com desvio Xt = 1 +Xt−1 + ǫt , ǫt ∼ N(0, 1).

Caṕıtulo 3

Modelos Probabiĺısticos

3.1 Introdução

Neste caṕıtulo serão descritos vários modelos adequados para dados de séries tempo-

rais. Tais modelos são chamados de processos estocásticos.

Matematicamente um processo estocástico pode ser definido como uma coleção de

variáveis aleatórias ordenadas no tempo e definidas em um conjunto de pontos T , que

pode ser cont́ınuo ou discreto. Iremos denotar a variável aleatória no tempo t porX(t)

no caso cont́ınuo (usualmente −∞ < t < ∞), e por Xt no caso discreto (usualmente t = 0,±1,±2, . . . ). O conjunto de posśıveis valores do processo é chamado de espaço de estados que pode ser discreto (e.g. o número de chamadas que chegam a uma

central telefônica a cada 2 horas) ou cont́ınuo (e.g. a temperatura do ar em uma

localidade observada em intervalos de 1 hora).

Em análise de séries temporais a situação é bem diferente da maioria dos problemas

estat́ısticos. Embora seja posśıvel variar o tamanho da série observada, usualmente

será imposśıvel fazer mais do que uma observação em cada tempo. Assim, tem-se

apenas uma realização do processo estocástico e uma única observação da variável

aleatória no tempo t denotada por x(t) no caso cont́ınuo e xt, para t = 1, . . . , N no

caso discreto.

Uma maneira de descrever um processo estocástico é através da distribuição

de probabilidade conjunta de X(t1), . . . , X(tk) para qualquer conjunto de tempos

t1, . . . , tk e qualquer valor de k. Esta é uma tarefa extremamente complicada e na

prática costuma-se descrever um processo estocástico através das funções média, va-

riância e autocovariância. Estas funções são definidas a seguir para o caso cont́ınuo

sendo que definições similares se aplicam ao caso discreto.

média µ(t) = E[X(t)]

variância σ2(t) = V ar[X(t)]

autocovariância γ(t1, t2) = E[X(t1)− µ(t1)][X(t2)− µ(t2)]

Note que a função de variância é um caso especial da função de autocovariância

quando t1 = t2. Momentos de ordem mais alta do processo também ser definidos

18

3.2. PROCESSOS ESTACIONÁRIOS 19

mas são raramente utilizados na prática e as funções µ(t) e γ(t1, t2) são em geral

suficientes.

3.2 Processos Estacionários

Uma importante classe de processos estocásticos são os chamados processos estaci-

onários. A idéia intuitiva de estacionariedade foi introduzida no caṕıtulo anterior e

aqui será apresentada a definição formal.

Uma série temporal é dita estritamente estacionária se a distribuição de proba-

bilidade conjunta de X(t1), . . . , X(tk) é a mesma de X(t1 + τ), . . . , X(tk + τ). Ou

seja, o deslocamento da origem dos tempos por uma quantidade τ não tem efeito na

distribuição conjunta que portanto depende apenas dos intervalos entre t1, . . . , tk.

Em particular, para k = 1 a estacionariedade estrita implica que a distribuição

de X(t) é a mesma para todo t de modo que, se os dois primeiros momentos forem

finitos, temos que

µ(t) = µ e σ2(t) = σ2

são constantes que não dependem de t.

Para k = 2 a distribuição conjunta de X(t1) e X(t2) depende apenas da distância

t2 − t1, chamada defasagem. A função de autocovariância γ(t1, t2) também depende apenas de t2 − t1 e pode ser escrita como γ(τ) onde

γ(τ) = E[X(t)− µ][X(t+ τ)− µ] = Cov[X(t), X(t+ τ)]

é chamado de coeficiente de autocovariância na defasagem τ .

Note que o tamanho de γ(τ) depende da escala em que X(t) é medida. Portanto,

para efeito de interpretação, é mais útil padronizar a função de autocovariância dando

origem a uma função de autocorrelação

ρ(τ) = γ(τ)/γ(0)

que mede a correlação entre X(t) e X(t + τ). No caṕıtulo anterior foi apresentado

o seu equivalente emṕırico para séries discretas rk. Note também que o argumento

τ será discreto se a série temporal for discreta e cont́ınuo se a série temporal for

cont́ınua.

Na prática é muito dif́ıcil usar a definição de estacionariedade estrita e costuma-se

definir estacionariedade de uma forma menos restrita.

Definição 3.1. Um processo estocástico {X(t), t ∈ T} é dito ser estacionário de segunda ordem ou fracamente estacionário se a sua função média é constante e sua

função de autocovariância depende apenas da defasagem, i.e.

E[X(t)] = µ e Cov[X(t), X(t+ τ)] = γ(τ).

Nenhuma outra suposição é feita a respeito dos momentos de ordem mais alta. Além

disso, fazendo τ = 0 segue que V ar[X(t)] = γ(0), ou seja a variância do processo

20 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

assim como a média também é constante. Note também que tanto a média quanto a

variância precisam ser finitos.

Esta definição mais fraca de estacionariedade será utilizada daqui em diante já

que muitas propriedades dos processos estacionários dependem apenas da estrutura

especificada pelo primeiro e segundo momentos. Uma classe importante de processos

aonde isto se verifica é a classe de processos normais ou Gaussianos aonde a distribui-

ção conjunta de X(t1), . . . , X(tk) é normal multivariada para todo conjunto t1, . . . , tk.

A distribuição normal multivariada fica completamente caracterizada pelo primeiro e

segundo momentos, i.e. por µ(t) e γ(t1, t2), assim estacionariedade fraca implica em

estacionariedade estrita para processos normais. Por outro lado, µ e γ(τ) podem não

descrever adequadamente processos que sejam muito “não-normais”.

3.3 A Função de Autocorrelação

Foi visto na Seção 2.4 que os coeficientes de autocorrelação amostral de uma série

temporal observada são uma ferramenta importante para descrever a série. Analoga-

mente, a função de autocorrelação teórica (fac) de um processo estocástico estacio-

nário é uma ferramenta importante para assessar suas propriedades. A seguir serão

apresentadas propriedades gerais da função de autocorrelação.

Se um processo estocástico estacionário X(t) tem média µ e variância σ2 então

ρ(τ) = γ(τ)/γ(0) = γ(τ)/σ2

e portanto ρ(0) = 1. As seguintes propriedades são facilmente verificáveis.

1. A correlação entre X(t) e X(t + τ) é a mesma que entre X(t) e X(t − τ), ou seja ρ(τ) = ρ(−τ).

2. −1 < ρ(τ) < 1.

3. Embora um processo estocástico tenha uma estrutura de autocovariância única

o contrário não é verdadeiro em geral. É posśıvel encontrar vários processos com

a mesma função de autocorrelação, o que dificulta ainda mais a interpretação

do correlograma.

3.4 Alguns Processos Estocásticos

Nesta seção serão apresentados alguns processos estocásticos que são utilizados com

frequência na especificação de modelos para séries temporais.

3.4.1 Sequência Aleatória

Um processo em tempo discreto é chamado puramente aleatório se consiste de uma

sequência de variáveis aleatórias {ǫt} independentes e identicamente distribuidas. Isto implica nas seguintes propriedades

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 21

1. E(ǫt) = E(ǫt|ǫt−1, ǫt−2, . . . ) = µ

2. V ar(ǫt) = V ar(ǫt|ǫt−1, ǫt−2, . . . ) = σ2ǫ

3. γ(k) = Cov(ǫt, ǫt+k) = 0, k = ±1,±2, . . . .

Como a média e a função de autocovariância não dependem do tempo o processo é

estacionário em segunda ordem. A função de autocorrelação é simplesmente

ρ(k) =

{

1, k = 0

0, k = ±1,±2, . . . .

Um processo puramente aleatório é as vezes chamado de rúıdo branco e pode ser

útil por exemplo na construção de processos mais complicados. As propriedades

acima podem ser entendidas como ausência de correlação serial e homocedasticidade

condicional (variância condicional constante).

3.4.2 Passeio Aleatório

Seja {ǫt} um processo discreto puramente aleatório com média µ e variância σ2ǫ . Um processo {Xt} é chamada de passeio aleatório se

Xt = Xt−1 + ǫt.

Fazendo-se substituições sucessivas obtém-se que

Xt = Xt−2 + ǫt−1 + ǫt

= Xt−3 + ǫt−2 + ǫt−1 + ǫt ...

= X0 + t

j=1

ǫj

e iniciando o processo em X0 = 0 não é dif́ıcil verificar que

E(Xt) = t

j=1

E(ǫj) = tµ

V ar(Xt) = t

j=1

V ar(ǫj) = tσ 2 ǫ .

Além disso, a função de autocovariância é dada por

Cov(Xt, Xt−k) = Cov(ǫ1 + · · ·+ ǫt−k + · · ·+ ǫt, ǫ1 + · · ·+ ǫt−k) = (t− k)σ2ǫ

e portanto a função de autocorrelação fica

ρt(k) = t− k t

.

22 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

Como a média, a variância e as autocovariâncias dependem de t este processo é não

estacionário. No entanto, é interessante notar que a primeira diferença de um passeio

aleatório é estacionária já que

▽Xt = Xt −Xt−1 = ǫt.

Os exemplos mais conhecidos de séries temporais que se comportam como um

passeio aleatório são os preços de ações em dias sucessivos (ver por exemplo Morettin

e Toloi, 2004).

3.4.3 Processos de Média Móveis

Seja {ǫt} um processo discreto puramente aleatório com média zero e variância σ2ǫ . Um processo {Xt} é chamada de processo de médias móveis de ordem q, ou MA(q), se

Xt = ǫt + β1ǫt−1 + · · ·+ βqǫt−q, (3.1) sendo βi ∈ R, i = 1, . . . , q. Não é dif́ıcil verificar como ficam a média e a variância deste processo,

E(Xt) = E(ǫt) +

q ∑

j=1

βjE(ǫt−j) = 0

V ar(Xt) = V ar(ǫt) +

q ∑

j=1

β2jV ar(ǫt−j) = (1 + β 2 1 + · · ·+ β2q ) σ2ǫ .

Além disso, como Cov(ǫt, ǫs) = σ 2 ǫ para t = s e Cov(ǫt, ǫs) = 0 para t 6= s, a função

de autocovariância é dada por

γ(k) = Cov(Xt, Xt+k)

= Cov(ǫt + β1ǫt−1 + · · ·+ βqǫt−q, ǫt+k + β1ǫt+k−1 + · · ·+ βqǫt+k−q)

=

0 k > q

σ2ǫ

q−k ∑

j=0

βjβj+k k = 0, . . . , q

γ(−k) k < 0

(3.2)

com β0 = 1. Como a média e a variância são constantes e γ(k) não depende de t

o processo é (fracamente) estacionário para todos os posśıveis valores de β1, . . . , βq.

Além disso, se os ǫt’s forem normalmente distribuidos osXt’s também serão e portanto

o processo será estritamente estacionário.

A função de autocorrelação pode ser facilmente obtida de (3.2) como

ρ(k) =

1 k = 0 q−k ∑

j=0

βjβj+k

/ q ∑

j=0

β2j k = 1, . . . , q

0 k > q

ρ(−k) k < 0.

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 23

Note que a função tem um ponto de corte na defasagem q, i.e. ρ(k) = 0 para k >

q. Esta é uma caracteŕıstica espećıfica de processos médias móveis e será útil na

especificação do valor de q na prática (Box & Jenkins 1970, p. 170).

> MAacf <- function(q, beta, lag.max) {

+ sig2x = 1 + sum(beta^2)

+ rho = rep(0, lag.max)

+ for (k in 1:q) {

+ rho[k] = beta[k]

+ if (q - k > 0) {

+ for (j in 1:(q - k)) rho[k] = rho[k] + beta[j] *

+ beta[j + k]

+ }

+ rho[k] = rho[k]/sig2x

+ }

+ return(rho)

+ }

> round(MAacf(q = 2, beta = c(0.5, 0.3), lag.max = 6), 4)

[1] 0.4851 0.2239 0.0000 0.0000 0.0000 0.0000

Vamos analisar agora com mais detalhes o caso particular do processo MA(1). A

função de autocorrelação fica

ρ(k) =

1 k = 0

β1/(1 + β 2 1) k = ±1

0 k > 1.

(3.3)

O processo é estacionário para qualquer valor de β1 mas em geral é desejável impor

restrições para que ele satisfaça uma condição chamada inversibilidade. Considere os

seguintes processos MA(1)

Xt = ǫt + θǫt−1

Xt = ǫt + 1

θ ǫt−1.

Substituindo em (3.3) não é dif́ıcil verificar que estes dois processos diferentes têm

exatamente a mesma função de autocorrelação. Assim, não é posśıvel identificar um

processo MA(1) único a partir da função de autocorrelação. Por outro lado, podemos

fazer substituições sucessivas e reescrever estes dois processos colocando ǫt em função

de Xt, Xt−1, . . . , i.e.

ǫt = Xt − θXt−1 + θ2Xt−2 − θ3Xt−3 + . . .

ǫt = Xt − 1

θ Xt−1 +

1

θ2 Xt−2 −

1

θ3 Xt−3 + . . .

Se |θ| < 1 a primeira série converge e o modelo é dito ser inverśıvel mas a segunda não converge e o modelo é não inverśıvel. Ou seja, a condição de inversibilidade (neste

24 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

caso |θ| < 1) garante que existe um único processo MA(1) para uma dada função de autocorrelação. Outra consequência da inversibilidade é que o processo MA(1)

pode ser reescrito como uma regressão de ordem infinita nos seus próprios valores

defasados.

Para um processo MA(q) esta condição pode ser melhor expressa usando-se o

operador de retardo, denotado por B e definido como

BjXt = Xt−j , para todo j.

A equação (3.1) pode então ser reescrita como

Xt = (1 + β1B + β2B 2 + · · ·+ βqBq)ǫt = θ(B)ǫt

onde θ(B) é um polinômio de ordem q em B. Um processo MA(q) é inverśıvel se as

ráızes da equação

θ(B) = 1 + β1B + β2B 2 + · · ·+ βqBq = 0

estiverem fora do ćırculo unitário. Ou seja, se δ1, . . . , δq são q soluções de θ(B) = 0

então o processo é inverśıvel se |δi| > 1, i = 1, . . . , q. Teremos então 2q modelos com a mesma função de autocorrelação mas somente um deles será inverśıvel.

Finalmente, vale notar que uma constante µ qualquer pode ser adicionada ao lado

direito de (3.1) dando origem a um processo com média µ. O processo continuará

sendo estacionário com E(Xt) = µ e em particular a função de autocorrelação não

será afetada.

3.4.4 Processos Autoregressivos

Suponha que {ǫt} seja um processo puramente aleatório com média zero e variância σ2ǫ . Um processo {Xt} é chamada de processo autoregressivo de ordem p, ou AR(p), se

Xt = α1Xt−1 + · · ·+ αpXt−p + ǫt. (3.4)

Note a similaridade com um modelo de regressão múltipla, onde os valores passados

de Xt fazem o papel das regressoras. Assim, processos AR podem ser usados como

modelos se for razoável assumir que o valor atual de uma série temporal depende do

seu passado imediato mais um erro aleatório.

Por simplicidade vamos começar estudando em detalhes processos de primeira

ordem, AR(1), i.e.

Xt = αXt−1 + ǫt. (3.5)

Note que existe uma estrutura Markoviana no processo AR(1) no sentido de que, dado

Xt−1, Xt não depende de Xt−2, Xt−3, . . . . Fazendo subtituições sucessivas em (3.5)

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 25

obtemos que

Xt = α(αXt−2 + ǫt−1) + ǫt = α 2Xt−2 + αǫt−1 + ǫt

= α2(αXt−3 + ǫt−2) + αǫt−1 + ǫt

= α3Xt−3 + α 2ǫt−2 + αǫt−1 + ǫt

...

= αr+1Xt−r−1 + r

j=0

αjǫt−j .

Se Xt for estacionário com variância finita σ 2 X podemos escrever que

E[Xt − r

j=0

αjǫt−j ] 2 = α2r+2E(X2t−r−1) = α

2r+2σ2X

e se |α| < 1 temos que α2r+2 → 0 quando r → ∞. Portanto, esta condição nos permite escrever Xt como o seguinte processo MA infinito,

Xt = ǫt + αǫt−1 + α 2ǫt−2 + . . .

e assim |α| < 1 é uma condição suficiente para que Xt seja estacionário. Neste caso, reescrevendo o processo k peŕıodos à frente, i.e.

Xt+k = ǫt+k + αǫt+k−1 + · · ·+ αkǫt + . . . (3.6)

note como o efeito de ǫt sobre Xt+k diminui a medida que k aumenta e por isso é

chamado efeito transitório.

Podemos também usar o operador de retardo reescrevendo a equação (3.5) como

(1− αB)Xt = ǫt

ou equivalentemente

Xt = 1

(1− αB)ǫt = (1 + αB + α 2B2 + . . . )ǫt = ǫt + αǫt−1 + α

2ǫt−2 + . . .

Escrevendo o processo AR(1) neste formato de MA infinito fica fácil ver que a sua

média e variância são dados por

E(Xt) = 0 e V ar(Xt) = σ 2 ǫ (1 + α

2 + α4 + . . . ) = σ2ǫ

1− α2 .

A função de autocovariância pode ser obtida usando os resultados acima. Rees-

crevendo a equação (3.6) como

Xt+k = ǫt+k + · · ·+ αk−1ǫt+1 + αkǫt + αk+1ǫt−1 + αk+2ǫt−2 + . . .

pode-se verificar que, para qualquer k = 1, 2, . . . ,

Cov(ǫt + αǫt−1 + α 2ǫt−2 + . . . , ǫt+k + · · ·+ αk−1ǫt+1) = 0.

26 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

Portanto,

E(XtXt+k) = Cov(ǫt + αǫt−1 + α 2ǫt−2 + . . . , α

kǫt + α k+1ǫt−1 + α

k+2ǫt−2 + . . . )

= αkE(ǫ2t ) + α k+2E(ǫ2t−1) + α

k+4E(ǫ2t−2) + . . .

= αkσ2ǫ (1 + α 2 + α4 + . . . ) = αk

σ2ǫ 1− α2 = α

kσ2X = γ(k).

Assim, a função de autocorrelação é ρ(k) = αk para k = 0, 1, 2, . . . . Assim, como a

média e a variância são constantes e γ(k) não depende de t o processo AR(1) com

|α| < 1 é estacionário. Na Figura 3.1 são mostradas graficamente as autocorrelações teóricas de um pro-

cesso AR(1) até a defasagem k = 20 para α igual a 0,8, -0,8, 0,2 e -0,2. Note como

a função de autocorrelação decai rapidamente para zero quando α = 0, 2 e se alterna

entre valores positivos e negativos quando α = −0, 8. Ou seja sempre há um decai- mento exponencial para zero mas este decaimento depende do sinal e magnitude de

α.

5 10 15 20

0. 0

0. 2

0. 4

0. 6

0. 8

5 10 15 20

− 0.

8 −

0. 2

0. 2

0. 6

5 10 15 20

0. 00

0. 10

0. 20

5 10 15 20

− 0.

20 −

0. 10

0. 00

Figura 3.1: As 20 primeiras autocorrelações teóricas de um processo AR(1) com (a) α = 0, 8, (b) α = −0, 8, (c) α = 0, 2 e (d) α = −0, 2.

Generalizando os resultados acima para um processo AR(p) escrevemos novamente

Xt como um processo MA infinito com coeficientes ψ0, ψ1, . . . , i.e.

Xt = ψ0ǫt + ψ1ǫt−1 + ψ2ǫt−2 + · · · = (ψ0 + ψ1B + ψ2B2 + . . . )ǫt = ψ(B)ǫt.

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 27

e em analogia com o caso AR(1) segue que o processo será estacionário se ∑

j ψ 2 j < ∞.

Usando agora o operador de retardo a equação (3.4) fica

(1− α1B − α2B2 − · · · − αpBp)Xt = ǫt ou φ(B)Xt = ǫt

e portanto o processo AR(p) pode ser escrito como

Xt = φ(B) −1ǫt = ψ(B)ǫt.

Assim, os coeficientes ψj podem ser obtidos a partir dos coeficientes αj fazendo-se

(1− α1B − α2B2 − · · · − αpBp)(ψ0 + ψ1B + ψ2B2 + . . . ) = 1

Desenvolvendo-se esta expressão segue que

ψ0 + ψ1B + ψ2B 2 + · · · − α1ψ0B − α1ψ1B2 − α1ψ2B3 − . . .

− α2ψ0B2 − α2ψ1B3 − α2ψ2B4 − . . . ...

− αpψ0Bp − αpψ1Bp+1 − · · · = 1 + 0B + 0B2 + . . .

e agora agrupando em termos de B,B2, . . .

ψ0 + (ψ1 − α1ψ0)B + (ψ2 − α1ψ1 − α2ψ0)B2 + · · · = 1 + 0B + 0B2 + . . .

donde obtém-se os coeficientes MA recursivamente como

ψ0 = 1

ψ1 = ψ0α1

ψ2 = ψ1α1 + ψ0α2

ψ3 = ψ2α1 + ψ1α2 + ψ0α3 ...

ψi = i

j=1

ψi−jαj .

O efeito de ǫt sobre Xt+k é dado por ψk, k = 1, 2, . . . .

Pode-se mostrar que (ver por exemplo Box, Jenkins, & Reinsel 1994) a condição

de estacionariedade do processo Xt é que todas as ráızes de φ(B) = 0 estejam fora

do ćırculo unitário. Em particular, para p = 1 temos que φ(B) = 1−αB = 0 implica que B = 1/α e a condição de estacionariedade fica |α| < 1 conforme já haviamos verificado.

Para reescrever um processo AR(p) em forma vetorial, defina Zt =

(Xt−1, . . . , Xt−p) ′ e portanto

Zt = ΦZt−1 + ut

28 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

sendo a matriz Φ definida como

Φ =

φ1 φ2 . . . φp−1 φp 1 0 . . . 0 0

0 1 . . . 0 0 ...

... ...

... ...

0 0 . . . 1 0

e ut = (ǫt, 0, . . . , 0) ′.

Para obter a função de autocorrelação de um processo AR(p) é algebricamente

mais simples assumir a priori que o processo é estacionário com E(Xt) = 0, V ar(Xt) =

σ2X e Cov(Xt, Xt−k) = γ(k). Neste caso, multiplicando a equação (3.4) por Xt−k, i.e

XtXt−k = α1Xt−1Xt−k + · · ·+ αpXt−pXt−k + ǫtXt−k.

e tomando o valor esperado obtemos que

E(XtXt−k) = γ(k) = α1E(Xt−1Xt−k) + · · ·+ αpE(Xt−pXt−k) = α1γ(k − 1) + · · ·+ αpγ(k − p), k > 0.

Dividindo-se ambos os lados pela variância constante σ2X obtem-se que

ρ(k) = α1ρ(k − 1) + · · ·+ αpρ(k − p), k > 0

chamadas equações de Yule-Walker.

Por exemplo, para um processo AR(1) com coeficiente α segue que ρ(1) = α,

ρ(2) = αρ(1) = α2, . . . , ρ(k) = αk como já haviamos verificado. Para um processo

AR(2) com coeficientes α1 e α2 segue que

ρ(1) = α1ρ(0) + α2ρ(1) ⇒ ρ(1) = α1/(1− α2)

e as outras autocorrelaçõs são obtidas iterativamente como

ρ(k) = α1ρ(k − 1) + α2ρ(k − 2), k ≥ 2

Autocorrelações Parciais

Para um processo AR(p), o último coeficiente αp mede o “excesso de correlação” na

defasagem p que não é levado em conta por um modelo AR(p − 1). Este é chamado de p-ésimo coeficiente de autocorrelação parcial. Assim, variando k = 1, 2, . . . temos

a chamada função de autocorrelação parcial (FACP).

Por outro lado, em um processo AR(p) não existe correlação direta entre Xt e

Xt−p−1, Xt−p−2, . . . e substituindo k = p+ 1, p+ 2, . . . nas equações de Yule-Walker

obtem-se que todos os coeficientes de correlação parcial serão nulos para k > p. Por

exemplo, substituindo-se k = p+ 1 segue que

ρ(p+ 1) = α1ρ(p) + · · ·+ αpρ(1) + αp+1.

O fato de que a facp é igual a zero para k > p é sugerido em Box and Jenkins (1970,

p. 170) como uma ferramenta para determinar a ordem p do processo autoregressivo

para séries temporais observadas.

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 29

3.4.5 Modelos Mistos ARMA

Combinando-se modelos AR e MA pode-se obter uma representação adequada com

um número menor de parâmetros. Processos autoregressivos médias móveis (ARMA)

formam um classe de modelos muito úteis e parcimoniosos para descrever dados de

séries temporais. O modelo ARMA(p, q) é dado por

Xt = α1Xt−1 + · · ·+ αpXt−p + ǫt + β1ǫt−1 + · · ·+ βqǫt−q

onde {ǫt} é um processo puramente aleatório com média zero e variância σ2ǫ . Note que, modelos AR ou MA podem ser obtidos como casos especiais quando p = 0 ou

q = 0. Usando o operador de retardo o modelo pode ser reescrito como

(1− α1B − α2B2 − · · · − αpBp)Xt = (1 + β1B + β2B2 + · · ·+ βqBq)ǫt

ou

φ(B)Xt = θ(B)ǫt.

Os valores de α1, . . . , αp que tornam o processo estacionário são tais que as ráızes

de φ(B) = 0 estão fora do ćırculo unitário. Analogamente, os valores de β1, . . . , βq que tornam o processo inverśıvel são tais que as ráızes de θ(B) = 0 estão fora do

ćırculo unitário.

Vale notar que as funções de autocorrelação e autocorrelação parcial ficam con-

sideravelmente mais complicadas em processos ARMA. De um modo geral, para um

processo ARMA(p, q) estacionário a função de autocorrelação tem um decaimento

exponencial ou oscilatório após a defasagem q enquanto que a facp tem o mesmo

comportamento após a defasagem p (Box & Jenkins 1970, p. 79). Em prinćıpio este

resultado pode ser utilizado para auxiliar na determinação da ordem (p, q) do processo

mas na prática pode ser bastante dif́ıcil distinguir entre decaimentos exponenciais e

oscilatórios através das estimativas destas funções.

A Tabela 3.1 mostra as propriedades teóricas das funções de autocorrelação e au-

tocorrelação parcial para alguns processos estacionários como auxiliar na identificação

do modelo.

Tabela 3.1: Propriedades teóricas da fac e facp.

Processo FAC FACP

série aleatória 0 0

AR(1), α > 0 decaimento exponencial 0, k ≥ 2 AR(1), α < 0 decaimento oscilatório idem

AR(p) decaimento para zero 0, k > p

MA(1) 0, k > 1 decaimento oscilatório

ARMA(p, q) decaimento a partir de q decaimento a partir de p

30 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

3.4.6 Modelos ARMA Integrados

Os modelos discutidos até agora são apropriados para séries temporais estacionárias.

Assim, para ajustar estes modelos a uma série temporal observada é necessário remo-

ver as fontes de variação não estacionárias. Por exemplo, se a série observada for não

estacionária na média pode-se tentar remover a tendência tomando-se uma ou mais

diferenças (esta abordagem é muito utilizada em Econometria).

Um modelo ARMA no qual Xt é substituido pela sua d-ésima diferença ▽ dXt

é capaz de descrever alguns tipos de séries não estacionárias. Denotando a série

diferenciada por

Wt = ▽ dXt = (1−B)dXt

o processo autoregressivo integrado médias móveis denotado ARIMA(p, d, q) é dado

por

Wt = α1Wt−1 + · · ·+ αpWt−p + ǫt + β1ǫt−1 + · · ·+ βqǫt−q ou, equivalentemente

φ(B)(1−B)dXt = θ(B)ǫt. (3.7)

Da equação (3.7) acima pode-se notar que o modelo para Xt é claramente não

estacionário já que o polinômio autoregressivo φ(B)(1−B)d tem exatamente d ráızes sobre o ćırculo unitário, ou d ráızes unitárias. Um processo que se torna estacionário

após d diferenças é dito ser não estacionário homogêneo, ou integrado de ordem d,

I(d).

Na prática valores pequenos são em geral especificados para d, sendo d = 1 o valor

mais frequentemente utilizado e excepcionalmente d = 2. Note também que o passeio

aleatório pode ser considerado um processo ARIMA(0,1,0).

Vale notar que para dados reais um modelo ARIMA (e de fato qualquer modelo)

é no máximo uma aproximação para o verdadeiro processo gerador dos dados. Na

prática pode ser bem dif́ıcil distinguir entre um processo estacionário com memória

longa (e.g. AR(1) com α ≈ 1) e um processo não estacionário homogêneo. Existe uma vasta literatura econométrica sobre testes de ráız unitária (ver por exemplo

Hamilton 1994 e Bauwens, Lubrano, & Richard 1999). Mais recentemente, modelos

da classe ARFIMA (ou ARIMA fracionários) tem sido utilizados para analisar séries

com memória longa. Estes tópicos não serão abordados aqui e o leitor interessado

pode consultar por exemplo Brockwell & Davis (1991) além das referências acima.

Exerćıcios

Nos exerćıcios a seguir {ǫt} é um processo discreto puramente aleatório com média zero e variância σ2ǫ .

1. Encontre a fac do processo Xt = ǫt + 0, 7ǫt−1 − 0, 2ǫt−2.

2. Encontre a fac do processo Xt − µ = 0, 7(Xt−1 − µ) + ǫt.

3. Encontre a fac do processo Xt = 1 3Xt−1 +

2 9Xt−2 + ǫt.

3.4. ALGUNS PROCESSOS ESTOCÁSTICOS 31

4. Se Xt = µ+ ǫt + βǫt−1 mostre que a fac do processo não depende de µ.

5. Reescreva cada um dos modelos abaixo em termos de operador de retardo B e

verifique se o modelo é estacionário e/ou inverśıvel:

(a) Xt = 0, 3Xt−1 + ǫt.

(b) Xt = ǫt − 1, 3 ǫt−1 + 0, 4 ǫt−2. (c) Xt = 0, 5Xt−1 + ǫt − 1, 3 ǫt−1 + 0, 4 ǫt−2. (d) ▽Xt = 0, 3 ▽Xt−1 + ǫt − 0, 6 ǫt−1 (e) Xt = Xt−1 + ǫt − 1, 5ǫt−1

6. Mostre que o processo Xt = Xt−1 + cXt−2 + ǫt é estacionário se −1 < c < 0 e obtenha a fac para c = −3/16.

7. Mostre que o processo Xt = Xt−1 + cXt−2 − cXt−3 + ǫt é não estacionário para qualquer valor de c.

8. Descreva como deve se comportar a função de autocorrelação teórica para os

seguintes processos,

(a) AR(1) estacionário, para α = 0, 1, α = −0, 75 e α = 0, 99. (b) Médias móveis de ordem q.

(c) Como deveriam ficar as funções de autocorrelação e autocorrelação parcial

amostrais que identificam os processos acima?

9. Descreva como deveriam se comportar as funções de autocorrelação e autocor-

relação parcial amostrais para processos AR, MA e ARMA não sazonais.

10. Para o modelo (1−B)(1− 0, 2B)Xt = (1− 0, 5B)ǫt, identifique os valores de p, q, e d e verifique se o processo é estacionário e inverśıvel.

11. Mostre que a função de autocovariância de um processo AR(1) estacionário com

variância σ2X é dada por α kσ2X (Sugestão: use a expressão (3.2) com q → ∞)

12. Verifique se Xt = ∑t

j=1 ǫt é estacionário.

13. Mostre que a fac do processo Xt = aXt−1 + ǫt + bǫt−1 é dada por

ρ(1) = (1 + ab)(a+ b)

1 + b2 + 2ab ρ(k) = aρ(k − 1), k = 2, 3, . . .

14. Obtenha a função de autocovarância do processo

Xt = ǫt + 1

a ǫt−1 +

1

a2 ǫt−2 + · · ·+

1

am ǫt−m

sendo que 0 < a < 1.

32 CAPÍTULO 3. MODELOS PROBABILÍSTICOS

15. Se {Xt} é um processo estacionário obtenha a função de autocovariância de Yt = Xt −Xt−1.

16. Mostre que o processo Xt = (α+1)Xt−1−αXt−2+ ǫt tem exatamente uma raiz unitária e reescreva-o como um processo ARIMA(1,1,0).

17. Obtenha a função de autocorrelação do passeio aleatório Xt = Xt−1 + ǫt com

E(ǫt) = µ, V ar(ǫt) = σ 2 ǫ e Cov(ǫt, ǫs) = 0, ∀t 6= s.

18. Verifique se o processo {Yt} tal que P (Yt = 1) = P (Yt = −1) = 1/2 é estacio- nário. Obtenha sua média, variância e covariância.

19. Sejam os processos Yt = ǫt + θǫt−1, |θ| > 1 e {Xt} tal que Xt = 1 se Yt ≥ 0 e Xt = −1 se Yt < 0. Verifique se {Xt} e {Yt} são estacionários. Calcule a função de autocorrelação de {Xt}.

20. Verifique que o processo Yt = (−1)tǫt é estacionário e que Xt = Yt + ǫt não é estacionário.

21. Se {Xt} e {Yt} são independentes e estacionários verifique se Zt = αXt + βYt, α, β ∈ R também é estacionário.

22. Obtenha a representação MA(∞) de um processo AR(2) estacionário.

23. Obtenha a representação AR(∞) de um processo MA(1) inverśıvel.

Caṕıtulo 4

Estimação

No caṕıtulo anterior foram estudados modelos probabiĺısticos que podem ser utilizados

para descrever dados de séries temporais. Neste caṕıtulo será discutido o problema

de ajustar um modelo aos dados observados. A inferência será baseada na função de

autocorrelação.

Para um processo estacionário {Xt} (t = 1, . . . , n), a função de densidade de probabilidade conjunta de X1, . . . , Xn pode ser sempre fatorada como

p(x1, . . . , xn) = p(x1)p(xn, . . . , x2|x1) = p(x1)p(x2|x1)p(xn, . . . , x3|x2, x1) ...

= p(x1) n ∏

t=2

p(xt|xt−1, . . . , x1).

Em particular para um modelo ARMA(p, q), denotando o vetor de parâmetros por

θ=(α1, . . . , αp, β1, . . . , βq, σ 2 ǫ )

′ e destacando-se a densidade conjunta das p primeiras

realizações segue que

p(x1, . . . , xn|θ) = p(x1, . . . , xp|θ) n ∏

t=p+1

p(xt|xt−1, . . . , x1,θ)

= p(x1, . . . , xp|θ) n ∏

t=p+1

p(xt|xt−1, . . . , xp,θ). (4.1)

A última igualdade vem da estrutura Markoviana da componente autoregressiva. O

segundo termo em (4.1) é a densidade condicional conjunta de xp+1, . . . , xn dados os

valores iniciais x1, . . . , xp e define então uma função de verossimilhança condicional

enquanto p(x1, . . . , xn|θ) define a função de verossimilhança exata. Se for atribuida uma distribuição de probabilidades conjunta também para θ então

pelo Teorema de Bayes é posśıvel obter sua distribuição atualizada após os dados serem

observados (distribuição a posteriori),

p(θ|x) = p(x|θ)p(θ) p(x)

∝ p(x|θ)p(θ).

33

34 CAPÍTULO 4. ESTIMAÇÃO

4.1 Autocovariância e autocorrelação

O coeficiente de autocovariância amostral de ordem k foi definido na Seção 2.4 como

ck = n−k ∑

t=1

(xt − x)(xt+k − x)/n

que é o estimador usual do coeficiente de autocovariância teórico γ(k). As propri-

edades deste estimador não serão detalhadas aqui mas podem ser encontradas por

exemplo em Priestley (1981). Após obter as estimativas de γ(k) os coeficientes de

autocorrelação são então estimados como rk = ck/c0, k = 1, 2, . . . .

Aqui serão consideradas apenas as propriedades de rk quando a amostra vem de

um processo puramente aleatório (propriedades gerais podem ser obtidas em Kendall

et al. 1983, Caṕıtulo 48). Vimos na Seção 2.4.1 que o coeficiente de autocorrelação

amostral rk é assintoticamente normalmente distribuido, com média e variância dados

por

E(rk) ≈ −1/n e V ar(rk) ≈ 1/n.

e os limites de confiança aproximados de 95% frequentemente utilizados são dados

por ±1, 96/√n. No caso geral, limites de 100(1-α)% podem ser construidos como ±qα/2/

√ n sendo qα/2 o percentil α/2 da distribuição normal padrão.

Interpretando o correlograma

No Caṕıtulo 2 foram vistos alguns exemplos de correlogramas associados a caracteŕıs-

ticas de séries temporais observadas. O correlograma é útil também na identificação

do tipo de modelo ARIMA que fornece a melhor representação de uma série observada.

Um correlograma como o da Figura 2.7 por exemplo, aonde os valores de rk decaem

para zero de forma relativamente lenta, indica não estacionariedade e a série precisa

ser diferenciada. Para séries estacionárias o correlograma é comparado com as auto-

correlações teóricas de vários processos ARMA para auxiliar na identificação daquele

mais apropriado. Por exemplo, se r1 é significativamente diferente de zero e todos os

valores subsequentes r2, r3, . . . são próximos de zero então um modelo MA(1) é indi-

cado já que sua função de autocorrelção teórica se comporta assim. Por outro lado, se

r1, r2, r3, . . . parecem estar decaindo exponencialmente então um modelo AR(1) pode

ser apropriado.

Vale notar entretando que a interpretação de correlogramas é um dos aspectos

mais dif́ıceis da análise de séries temporais. A função de autocorrelação parcial é um

importante coadjuvante nesta etapa de identificação se houver termos autoregressivos

no modelo já que seus valores estimados tendem a ficar próximos de zero após a

defasagem p.

Vimos no Caṕıtulo 3 que para um processo ARMA(p, q) estacionário a função

de autocorrelação teórica terá um decaimento exponencial ou oscilatório após a de-

fasagem q enquanto que a função de autocorrelação parcial teórica terá o mesmo

comportamento após a defasagem p. Mas na prática esta distinção entre decaimentos

4.2. AJUSTANDO PROCESSOS AUTOREGRESSIVOS 35

exponenciais e oscilatórios através das estimativas destas funções pode ser bastante

dif́ıcil.

4.2 Ajustando Processos Autoregressivos

Para um processo AR de ordem p com média µ dado por

Xt − µ = α1(Xt−1 − µ) + · · ·+ αp(Xt−p − µ) + ǫt,

e dadas n observações x1, . . . , xn, os parâmetros µ, α1, , . . . , αp podem ser estimados

pelo método de mı́nimos quadrados, i.e. minimizando-se a soma de quadrados

S = n ∑

t=p+1

[(xt − µ)− α1(xt−1 − µ)− · · · − αp(xt−p − µ)]2

com respeito a µ, α1, , . . . , αp. Note que o somatório é de t = p + 1 em diante, mas

esta pequena perda de informação não será importante se a série não for muito curta.

Além disso, se o processo ǫt tiver distribuição normal então as estimativas de mı́nimos

quadrado coincidem com as estimativas de máxima verossimilhança condicionada nas

p primeiras observações.

Alternativamente, dois métodos aproximados podem ser utilizados tomando-se

µ̂ = x. O primeiro ajusta os dados ao modelo

Xt − x = α1(Xt−1 − x) + · · ·+ αp(Xt−p − x) + ǫt,

como se fosse um modelo de regressão linear múltipla.

No segundo método os coeficientes de autocorrelação ρ(k) são substituidos pelas

suas estimativas rk nas p primeiras equações de Yule-Walker. Ou seja, estamos usando

o métodos dos momentos e por isto os estimadores resultantes são assintoticamente

equivalentes aos estimadores de máxima verossimilhança. Assim, temos um sistema

com p equações e p incógnitas α1, . . . , αp, i.e.

r1 = α1 + α2r1 + · · ·+ αprp−1 r2 = α1r1 + α2 + · · ·+ αprp−2

...

rp = α1rp−1 + α2rp−2 + · · ·+ αp

ou equivalentemente, 

r1 r2 ...

rp

=

1 r1 . . . rp−1 r1 1 . . . rp−2 ...

... ...

rp−1 rp−2 . . . 1

α1 α2 ...

αp

Exemplo 4.1 : Usando os comandos do R abaixo vamos simular um processo AR(3)

e usar as equações de Yule-Walker para estimar os coeficientes.

36 CAPÍTULO 4. ESTIMAÇÃO

> x = arima.sim(n = 200, model = list(ar = c(0.6, -0.7, 0.2)))

> r = acf(x, plot = FALSE)$acf[2:4]

> R = diag(3)

> R[1, 2] = R[2, 1] = r[1]

> R[1, 3] = R[3, 1] = r[2]

> R[2, 3] = R[3, 2] = r[1]

> round(solve(R, r), 4)

[1] 0.6659 -0.7513 0.2678

Para estimação por minimos quadrados basta escrever o AR(p) como um modelo

linear usual e resolver um sistema de equações lineares. Definindo-se

y =

xp+1 xp+2 ...

xn

X =

xp . . . x1 xp−1 . . . x2 ...

...

xn−1 . . . xn−p

ǫ =

ǫp+1 ǫp+2 ...

ǫn

α =

α1 α2 ...

αp

podemos reescrever o modelo na forma matricial como

y = Xα+ ǫ, (4.2)

sendo E(ǫ) = 0, V ar(ǫ) = σ2ǫ In−p e In−p a matriz identidade de ordem n − p. A solução de mı́nimos quadrados para os coeficientes α é obtida minimizando-se ǫ′ǫ e é

dada por α̂ = (X ′X)−1X ′y. Usando o valor estimado de α na equação do modelo

calcula-se os reśıduos como y = Xα̂, i.e.

et = xt − p

j=1

α̂jxt−j , t = p+ 1, . . . , n

e a estimativa de mı́nimos quadrados de σ2ǫ é dada por

σ̂2ǫ = 1

n− p

n ∑

t=p+1

e2t .

Note que os reśıduos também foram calculados a partir de t = p+ 1.

Mantendo a representação (4.2) e adicionando a hipótese de normalidade dos erros,

i.e. ǫ ∼ N(0, σ2ǫ In−p) obtém-se uma função de verossimilhança aproximada dada por,

L(α, σ2ǫ ) ∝ (σ2ǫ )−(n−p)/2 exp{−σ−2ǫ (y −Xα)′(y −Xα)/2}.

Neste caso, os EMV de α e σ2ǫ coincidem com os estimadores de mı́nimos quadrados,

∂ log(L(α, σ2ǫ ))

∂α = −σ

−2 ǫ

2

∂(y −Xα)′(y −Xα) ∂α

= −σ −2 ǫ

2

∂(−2α′X ′y +α′X ′Xα) ∂α

= −σ −2 ǫ

2 (−2X ′y + 2X ′Xα).

4.2. AJUSTANDO PROCESSOS AUTOREGRESSIVOS 37

∂ log(L(α, σ2ǫ ))

∂α

α=α̂ = 0 ⇐⇒ α̂ = (X ′X)−1X ′y.

Lembrando que (y −Xα̂)′(y −Xα̂) = ∑n

t=p+1 e 2 t segue que

∂ log(L(α̂, σ2ǫ ))

∂σ2ǫ = −1

2

∂σ2ǫ

(n− p) log(σ2ǫ ) + σ−2ǫ n ∑

t=p+1

e2t

= −1 2

(n− p)σ−2ǫ − σ−4ǫ n ∑

t=p+1

e2t

∂ log(L(α̂, σ2ǫ ))

∂σ2ǫ

σ2ǫ=σ̂ 2 ǫ

= 0 ⇐⇒ σ̂2ǫ = 1

n− p

n ∑

t=p+1

e2t .

Exemplo 4.2 : Para um modelo AR(1) com erros normais a matriz X tem somente

uma coluna e não é dif́ıcil verificar que

X ′X = n ∑

t=2

x2t−1 e X ′y =

n ∑

t=2

xtxt−1.

Portanto, o EMV condicional é dado por

α̂ =

∑n t=2 xtxt−1

∑n t=2 x

2 t−1

e σ̂2ǫ = 1

n− 1

n ∑

t=2

(xt − α̂xt−1)2.

Exemplo 4.3 : Novamente para o modelo AR(1) com erros normais o EMV incon-

dicional é obtido maximizando-se da função de verossimilhança exata. A expressão

(4.1) com p = 1 fica

p(x1, . . . , xn|α, σ2ǫ ) = p(x1|α, σ2ǫ ) n ∏

t=2

p(xt|xt−1, α, σ2ǫ ).

Lembrando que E(Xt) = 0 e V ar(Xt) = σ 2 ǫ /(1 − α2) e razoável assumir que X1 ∼

N(0, σ2ǫ /(1− α2)). Segue então que

L(α, σ2ǫ ) ∝ (

σ2ǫ 1− α2

)−1/2

exp

{

−1− α 2

2σ2ǫ x21

}

×

(σ2ǫ ) −(n−1)/2 exp

{

− 1 2σ2ǫ

n ∑

t=2

(xt − αxt−1)2 }

∝ (1− α2)1/2(σ2ǫ )−n/2 exp {

− 1 2σ2ǫ

(

(1− α2)x21 + n ∑

t=2

(xt − αxt−1)2 )}

.

Maximizar esta expressão (ou seu logaritmo) em relação a α requer algum algoritmo

de otimização numérica (por exemplo métodos de Newton-Raphson). No R podemos

usar a função optim como no Exemplo 4.4.

38 CAPÍTULO 4. ESTIMAÇÃO

Exemplo 4.4 : Foram gerados 200 valores de um processo AR(1) com parâmetros

α = 0, 8 e σ2ǫ = 1. Os comandos abaixo podem ser usados para obter as estimativas de

máxima verossimilhança (incondicional). Note que estamos maximizando o logaritmo

da verossimilhança e verificando a condição de estacionariedade.

> fun = function(theta, x) {

+ s2 = theta[1]

+ alpha = theta[2]

+ if (abs(alpha) >= 1)

+ return(-Inf)

+ n = length(x)

+ e = x[2:n] - alpha * x[1:(n - 1)]

+ Q = (1 - alpha^2) * x[1]^2 + sum(e^2)

+ return(-0.5 * (n * log(s2) - log(1 - alpha^2) + Q/s2))

+ }

> x = arima.sim(n = 200, model = list(ar = 0.8))

> init = c(1, 0.5)

> out = optim(init, fn = fun, method = "BFGS", control = list(fnscale = -1),

+ hessian = T, x = x)

> out$par

[1] 1.0196881 0.8274672

Como o custo computacional de estimar modelos AR não é tão grande uma

abordagem alternativa para determinação de p consiste em estimar modelos de ordem

progressivamente mais alta e calcular a soma de quadrados residual para cada valor

de p. Pode ser posśıvel encontrar o valor de p para o qual a inclusão de termos extras

não melhora sensivelmente o ajuste. Como vimos na Seção 3.4.4 este procedimento

dá origem à função de autocorrelação parcial.

Suponha agora que vamos atribuir uma distribuição de probabilidades para o

vetor de parâmetros θ = (α1, . . . , αp, σ 2 ǫ ). Pelo Teorema de Bayes e usando a verossi-

milhança condicional segue que

p(θ|x) ∝ p(θ) (σ2ǫ )−(n−p)/2 exp{−σ−2ǫ (y −Xα)′(y −Xα)/2}.

Para representar a informação a priori sobre θ pode-se fazer por exemplo, p(θ) =

p(α|σ2ǫ )p(σ2ǫ ) com α|σ2ǫ ∼ N(0, σ2ǫ Ip) ou p(θ) = p(α)p(σ2ǫ ) com α ∼ N(0, Ip). Nos dois casos comumente assume-se que σ2ǫ tem distribuição Gama Inversa, i.e. σ

2 ǫ ∼

GI(a, b) (ver Apêndice A), ou equivalentemente σ−2ǫ ∼ Gama(a, b).

Exemplo 4.5 : No modelo AR(1) com erros normais vamos atribuir as seguintes

distribuições a priori, α ∼ N(0, 1) e σ2ǫ ∼ GI(1, 1). Portanto,

p(α) ∝ exp(−α2/2) e p(σ2ǫ ) ∝ (σ2ǫ )−2 exp(−1/σ2ǫ )

4.3. AJUSTANDO PROCESSOS MÉDIAS MÓVEIS 39

e os comandos abaixo podem ser usados para obter a moda da distribuição a posteriori

conjunta de σ2ǫ e α.

> prior = function(theta) {

+ s2 = theta[1]

+ alpha = theta[2]

+ return(-alpha^2/2 - 1/s2 - 2 * log(s2))

+ }

> post = function(theta, x) fun(theta, x) + prior(theta)

> out = optim(init, fn = post, method = "BFGS", control = list(fnscale = -1),

+ hessian = T, x = x)

> out$par

[1] 1.0095348 0.8262561

Note que as estimativas pontuais nos Exemplos 4.4 e 4.5 são bastante similares.

Nenhuma restrição de estacionariedade foi imposta na distribuição a priori, mas é

posśıvel fazer uma otimização restrita ou mesmo impor esta restrição a priori. No caso

do AR(1) poderiamos atribuir uma distribuição normal truncada ou uma distribuição

uniforme em (-1,1) para o parâmetro α.

4.3 Ajustando Processos Médias Móveis

O problema de estimação dos parâmetros em modelos MA é bem mais complicado do

que em modelos AR. Os erros ǫt são agora funções não lineares complicadas dos parâ-

metros β1, . . . , βq e expressões anaĺıticas para os estimadores não podem ser obtidas.

Assim, métodos computacionais iterativos precisam ser utilizados para minimizar a

soma de quadrados residual.

Dado um modelo MA(q)

Xt = µ+ ǫt + β1ǫt−1 + · · ·+ βqǫt−q

e uma série observada x1, . . . , xn o procedimento iterativo consiste basicamente em

fixar os valores de µ, β1, . . . , βq e calcular os reśıduos

et = xt − µ− β1ǫt−1 − · · · − βqǫt−q

sequencialmente para t = 1, . . . , n assumindo que ǫ0 = ǫ−1 = · · · = ǫ−q+1 = 0 e substituindo ǫt−1, . . . , ǫt−q pelos residuos calculados. Assim,

e1 = x1 − µ e2 = x2 − µ− β1e1 = x2 − µ− β1x1 + β1µ e3 = x3 − µ− β1e2 − β2e1

...

40 CAPÍTULO 4. ESTIMAÇÃO

Dados estes reśıduos pode-se calcular a soma de quadrados residual S(µ,β) = ∑n

t=1 e 2 t . Repetindo este procedimento para µ, β1, . . . , βq variando em uma grade

de pontos pode-se escolher os valores que minimizam a soma de quadrados. Este

procedimento requer o uso de algoritmos eficientes de otimização numérica e nada

garante a sua convergência para um mı́nimo global.

Além das estimativas pontuais, se o processo {ǫt} tem distribuição normal então Box & Jenkins (1970), p. 228 descrevem regiões de confiança para os parâmetros do

modelo. Neste caso, se ǫt ∼ N(0, σ2ǫ ) a função de verossimilhança fica,

L(µ,β, σ2ǫ ) =

n ∏

t=1

(2πσ2ǫ ) −1/2 exp

{

− 1 2σ2ǫ

e2t

}

∝ (σ2ǫ )−n/2 exp {

− 1 2σ2ǫ

n ∑

t=1

e2t

}

.

e os valores de et são calculados como anteriormente. Portanto L(µ,β, σ 2 ǫ ) é uma

função não linear dos parâmetros.

Em termos práticos, se o procedimento de otimização utilizado levar muitas ite-

rações para convergir ou mesmo não convergir deve-se “desconfiar” das estimativas.

Neste caso as estimativas podem ser instáveis no sentido de que adicionando-se ou

removendo-se uma ou duas observações pode-se obter valores muito diferentes. Nesta

situação pode ser computacionalmente mais vantajoso ajustar um modelo AR aos

dados mesmo que o modelo resultante tenha mais parâmetros do que o modelo MA

sugerido pela função de autocorrelação.

4.4 Ajustando Processos ARMA

Os problemas de estimação para modelos ARMA são similares aqueles para modelos

MA no sentido de que um procedimento iterativo precisa ser utilizado. Isto ocorre

porque os erros {ǫt} são funções não lineares complicadas de todos os coeficientes α1, . . . , αp, β1, . . . , βq. Portanto os mesmos comentários da seção anterior são válidos

para procedimentos que levam muitas iterações para convergir, i.e deve-se“desconfiar”

das estimativas. Os residuos são calculados de forma análoga ao modelo MA (ver

Exerćıcio 13).

Outra dificuldade, espećıfica de modelos ARMA, é o problema de cancelamento

de ráızes. Por exemplo considere o modelo ARMA(2,1)

Xt = 2θXt−1 − θ2Xt−2 − φǫt−1 + ǫt

que pode ser reescrito em termos do operador de retardo como

(1− θB)2Xt = (1− φB)ǫt.

Note como θ = φ implica em um modelo AR(1) Xt = θXt−1 + ǫt, ou seja ambos os

modelos implicam exatamento no mesmo comportamento para a série temporal Xt.

4.5. MODELOS SAZONAIS 41

Este é um problema de identificação que fica ainda mais complicado em modelos de

ordem mais alta.

Em termos práticos é dif́ıcil identificar o problema de cancelamento de ráızes a não

ser, como já foi dito, que o procedimento iterativo deverá ter convergência lenta. No

caso particular de um modelo ARMA(1,1) deve-se “desconfiar” quando as estimativas

de α e β são muito similares. Para outros valores de p e q a única sugestão para tentar

minimizar o problema é não incluir muitos parâmetros no modelo.

Exemplo 4.6 : Vamos simular um processo ARMA(1,1) com ráızes similares e veri-

ficar o problema de cancelamento de ráızes.

> x = arima.sim(n = 100, list(ar = 0.7, ma = -0.75))

> arima(x, order = c(1, 0, 1), include.mean = F)

Call:

arima(x = x, order = c(1, 0, 1), include.mean = F)

Coefficients:

ar1 ma1

0.5992 -0.6835

s.e. 0.2946 0.2613

sigma^2 estimated as 1.020: log likelihood = -142.88, aic = 291.75

Note como as estimativas dos coeficientes estão muito diferentes dos valores verdadei-

ros e os erros padrões estão enormes!

4.5 Modelos Sazonais

Muitas séries temporais contém uma componente periódica sazonal que se repete a

cada s observações (s > 1). Por exemplo, com dados mensais e s = 12 tipicamente

espera-se que Xt dependa de Xt−12 e talvez de Xt−24 além de Xt−1, Xt−2, . . . .

Neste caso, tomar a primeira diferença xt − xt−1 não é suficiente para tornar a série (aproximadamente) estacionária. A forma apropriada de diferenciar dados com

padrão sazonal acentuado é tomar diferenças no peŕıodo sazonal. Por exemplo, para

dados mensais a primeira diferença sazonal é

▽12xt = (1−B12)xt = xt − xt−12 e terá variabilidade menor do que a primeira diferença não sazonal ▽xt = xt − xt−1, sendo portanto mais fácil de identificar e estimar.

Em geral, uma diferença sazonal é denotada por ▽s onde s é o peŕıodo sazonal.

A D-ésima diferença sazonal é então denotada por ▽Ds . Combinando-se os dois tipos

de diferenciação obtem-se o operador ▽d▽Ds . Por exemplo, tomando-se 1 diferença

simples e 1 sazonal em uma série mensal tem-se que

▽▽12xt = xt − xt−1 − xt−12 + xt−13

42 CAPÍTULO 4. ESTIMAÇÃO

Box & Jenkins (1970) generalizaram o modelo ARIMA para lidar com sazonalidade

e definiram um modelo ARIMA sazonal multiplicativo, denominado SARIMA, dado

por

φ(B)Φ(Bs)Wt = θ(B)Θ(B s)ǫt (4.3)

onde

φ(B) = (1− α1B − · · · − αpBp) Φ(Bs) = (1− φsBs − · · · − φPBPs)

Wt = ▽ d ▽ D s Xt

θ(B) = (1 + β1B + · · ·+ βqBq) Θ(Bs) = (1 + θsB

s + · · ·+ θQBQs).

Este modelo é chamado SARIMA multiplicativo de ordem (p, d, q)×(P,D,Q)s e parece extremamente complicado à primeira vista mas na prática os valores de d e

D em geral não serão maiores do que 1 e um número pequeno de coeficientes será

suficiente. Por exemplo, com P = 1 temos que

Φ(Bs) = (1− αsBs)

o que significa simplesmente que Wt depende de Wt−s. A série Wt é formada a partir

da série original tomando-se diferenças simples para remover a tendência e diferenças

sazonais para remover a sazonalidade.

Para fixar idéias considere o modelo SARIMA(1,0,0)× (0, 1, 1)12 para dados men- sais. Ou seja temos um termo autoregressivo e um termo média móvel sazonal mode-

lando a primeira diferença sazonal. O modelo pode ser escrito como

(1− αB)(1−B12)Xt = (1− θB12)ǫt

e desenvolvendo os produtos obtemos que

Xt = Xt−12 + α(Xt−1 −Xt−13) + ǫt + θǫt−12.

Assim, Xt depende de Xt−1, Xt−12 e Xt−13 além do erro no tempo t− 12. Para finalizar, ao ajustar um modelo sazonal aos dados a primeira tarefa é es-

pecificar os valores de d e D que tornam a série (aproximadamente) estacionária e

remove a maior parte da sazonalidade. Como já foi dito, estes valores raramente serão

maiores do que 1. Posteriormente os valores de p, P , q e Q devem ser especificados

com base nas funções de autocorrelação e autocorrelação parcial da série diferenciada.

Os valores de P e Q são especificados basicamente a partir de rk, k = s, 2s, . . . . Por

exemplo, para dados mensais se r12 é grande mas r24 é pequeno isto sugere que um

termo média móvel sazonal pode ser adequado.

Após ter identificado, por tentativa, o que parece ser um modelo SARIMA razoável

os parâmetros serão estimados por algum procedimento iterativo similar àqueles pro-

postos para modelos ARMA. Detalhes sobre as rotinas de estimação destes modelos

não serão abordados aqui e podem ser obtidos em Box & Jenkins (1970).

4.6. ADEQUAÇÃO DO MODELO 43

4.6 Adequação do Modelo

Todos os modelos são errados mas alguns são úteis (George Box)

Após identificar a ordem e estimar eficientemente os parâmetros de um modelo é

necessário verificar sua adequação antes de utilizá-lo por exemplo para fazer previsões.

Pode-se fazer testes de sobreajustamento, que consistem em incluir parâmetros extras

no modelo e verificar sua significância estat́ıstica. No caso de modelos ARMA deve-se

incluir um parâmetro de cada vez para evitar o problema de cancelamento de ráızes

mencionado na Seção 4.4.

4.6.1 Análise dos Reśıduos

Após um modelo ter sido ajustado a uma série temporal deve-se verificar se ele fornece

uma descrição adequada dos dados. Assim como em outros modelos estat́ısticos a idéia

é verificar o comportamento dos reśıduos onde

residuo = observação - valor ajustado.

Para os modelos vistos aqui o valor ajustado é a previsão 1 passo a frente de modo

que o reśıduo fica definido como o erro de previsão 1 passo a frente. Por exemplo,

em um modelo AR(1) se α̂ é a estimativa do coeficiente autoregressivo então o valor

ajustado no tempo t é α̂xt−1 e o reśıduo correspondente é et = xt − α̂xt−1. Se o modelo tiver um “bom” ajuste espera-se que os reśıduos se distribuam alea-

toriamente em torno de zero com variância aproximadamente constante e sejam não

correlacionados. Se a variância dos reśıduos for crescente uma transformação logaŕıt-

mica nos dados pode ser apropriada. O fenômeno de “não constância” na variância

é denominado de volatilidade na literatura de séries temporais e pode ser tratado

através de transformações nos dados (e.g. transformações de Box-Cox)1.

Além disso, em modelos de séries temporais os reśıduos estão ordenados no tempo

e é portanto natural tratá-los também como uma série temporal. É particularmente

importante que os reśıduos de um modelo estimado sejam serialmente (i.e. ao longo

do tempo) não correlacionados. Evidência de correlação serial nos reśıduos é uma

indicação de que uma ou mais caracteŕısticas da série não foi adequadamente descrita

pelo modelo.

Consequentemente, duas maneiras óbvias de verificar a adequação do modelo con-

sistem em representar graficamente os reśıduos e o seu correlograma. O gráfico tem-

poral poderá revelar a presença de dados discrepantes, efeitos de autocorrelação ou

padrões ćıclicos enquanto que o correlograma permite uma análise mais detalhada da

estrutura de autocorrelação indicando posśıveis termos faltantes no modelo.

Ou seja, assim como em outros modelos estat́ısticos, a idéia é que os reśıduos

poderão identificar caracteŕısticas que não foram adequadamente modeladas. Por

exemplo, autocorrelações residuais significativas nas defasagens 1 ou 2, ou em defa-

sagens sazonais (e.g. 12 para dados mensais) são uma indicação de que mais termos

1Uma tendência mais recente no entanto consiste em tentar modelar simultaneamente a média e

a variância ao invés de usar transformações.

44 CAPÍTULO 4. ESTIMAÇÃO

médias móveis devem ser incluidos no modelo. Por outro lado, um valor de rk li-

geiramente fora dos limites de confiança em defasagens sem significado óbvio (e.g.

k=5) não é indicação suficiente para se rejeitar o modelo. O mesmo comentário vale

para as autocorrelações parciais dos reśıduos no que diz respeito à inclusão de termos

autoregressivos (sazonais e não sazonais).

4.6.2 Testes sobre os reśıduos

Ao invés de olhar para as autocorrelações residuais individualmente pode-se testar

se um grupo de autocorrelações é significativamente diferente de zero através das

chamadas estat́ısticas Q. Para modelos ARMA Box & Jenkins (1970) sugeriram o

uso do teste de Box-Pierce para as hipóteses

H0 : ρ(1) = · · · = ρ(m) = 0 H1 : ρ(k) 6= 0,para algum k ∈ {1, . . . ,m}.

sendo a estat́ıstica de teste dada por

Q = n m ∑

k=1

r2k.

Na prática o número m de autocorrelações amostrais é tipicamente escolhido entre

15 e 30. Se o modelo ajustado for apropriado então Q terá distribuição aproximada-

mente qui-quadrado com m− p− q graus de liberdade. Assim, valores grandes de Q fornecem indicação contra a hipótese de que as autocorrelações são todas nulas, em

favor da hipótese de que ao menos uma delas é diferente de zero.

O teste de Box-Pierce não tem bom desempenho em amostras pequenas ou mo-

deradas no sentido de que a distribuição se afasta da qui-quadrado. Vários testes

alternativos foram sugeridos na literatura e o mais conhecido é o teste de Ljung-Box,

aonde a estat́ıstica de teste é dada por

Q = n(n+ 2) m ∑

k=1

r2k n− k .

Sua distribuição amostral também é aproximadamente qui-quadrado com

m− p− q graus de liberdade.

Exemplo 4.7 : Considere novamente a série com os totais mensais de passageiros

em linhas aéreas internacionais nos EUA entre 1949 e 1960 que aparece na Figura ??.

Existe uma clara tendência de crescimento bem como um padrão sazonal ao longo

dos anos. Foi feita uma transformação logaritmica nos dados (esta transformação

é sugerida na literatura). Faça os gráficos da FAC amostral da série original, 1a

diferença e 1a diferença sazonal. Os comandos abaixo podem ser utilizados e obtém-

se os gráficos da Figura 4.1.

> y = log(AirPassengers)

> z = cbind(y, diff(y), diff(y, lag = 12))

> yl = c("No de passageiros", "Variacao mensal", "Variacao anual")

4.6. ADEQUAÇÃO DO MODELO 45

Os gráficos anteriores indicam que precisamos tomar 1 diferença simples mais 1

diferença sazonal para tentar induzir estacionariedade aproximada.

> z = diff(diff(y), lag = 12)

> m = acf(z, lag.max = 36, plot = F)

> m$lag = m$lag * 12

Note que há valores grandes nas defasagens 1, 3, 12 e 23 do último gráfico. Isto

pode ser uma indicação de que termos MA sazonais e não sazonais devem ser incluidos

no modelo. Um modelo candidato para o logaritmo da série é SARIMA(0,1,1)x(0,1,1)

e foi estimado usando os comandos abaixo.

> m = arima(y, order = c(0, 1, 1), seasonal = list(order = c(0,

+ 1, 1)))

> m

Call:

arima(x = y, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1)))

Coefficients:

ma1 sma1

-0.4018 -0.5569

s.e. 0.0896 0.0731

sigma^2 estimated as 0.001348: log likelihood = 244.7, aic = -483.4

Como primeiro verificação da adequação do modelo vamos usar a função tsdiag()

que retorna os gráficos dos residuos padronizados, o correlograma e os p-valores do

teste de Ljung-Box para autocorrelações de ordem 1, 2, . . . . O resultado está na Figura

4.3.

Compare estes p-valores com o resultado da função Box.test() que calcula as

estatisticas de Box-Pierce e Ljung-Box para a hipótese nula de independência.

> for (i in 1:10) {

+ b = Box.test(m$residuals, i, type = "Ljung-Box")$p.value

+ print(b)

+ }

[1] 0.8610215

[1] 0.945251

[1] 0.4829255

[1] 0.3663102

[1] 0.4320234

[1] 0.488321

[1] 0.539204

[1] 0.6328112

[1] 0.5096084

[1] 0.5502512

46 CAPÍTULO 4. ESTIMAÇÃO

Testando a Normalidade dos Reśıduos

Para uma variável aleatória X tal que E(X) = µ e V ar(X) = σ2 define-se os coefici-

entes de assimetria e curtose como,

A(X) = E

(

(X − µ)3 σ3

)

e K(X) = E

(

(X − µ)4 σ4

)

respectivamente. A distribuição normal tem assimetria 0 e curtose igual a 3. Substi-

tuindo os momentos teóricos de X pelos seus equivalente amostrais

mj = 1

n

n ∑

t=1

(Xt −X)j

os estimadores da assimetria e curtose são dados por

 = m3 √

m32 e K̂ =

m4 √

m22

respectivamente. Sob a hipótese de normalidade as variáveis aleatórias √

n/6Â e √

n/24(K̂ − 3) são independentes e têm distribuição assintótica N(0, 1) e assim a estat́ıstica

nÂ2

6 +

n(K̂ − 3)2 24

tem distribuição assintótica χ2 com 2 graus de liberdade e pode ser usada para testar

a normalidade de X.

As outras verificações usuais sobre os residuos também devem ser feitas. Por

exemplo, um histograma com curva normal superposta, o gráfico de probabilidades

normais e um teste de normalidade. Os comandos abaixo podem ser utilizados no R.

> z = m$residuals

> d = seq(range(z)[1] - 3 * sd(z), range(z)[2] + 3 * sd(z), 0.001)

> a = shapiro.test(z)

Exerćıcios

1. Calcule as autocorrelações teóricas de um processo MA(Q) puramente sazonal.

2. Faça um esboço do correlograma para uma série com estrutura MA(Q) pura-

mente sazonal, i.e. não existe dependência dentro de um peŕıodo sazonal.

3. Para uma série temporal observada foi identificado o modelo ARIMA(1,1,1).

(a) Escreva o modelo em termos do operador de retardo.

(b) Descreva como deve ter sido o comportamento das funções de autocorrela-

ção e autocorrelação parcial da série original e da série diferenciada.

4. Escreva o modelo SARIMA(0, 0, 1) × (1, 1, 0)12 em termos de operador de re- tardo.

4.6. ADEQUAÇÃO DO MODELO 47

5. Para uma série mensal observada foi identificado e estimado o modelo

SARIMA(1,1,0)×(0,1,0).

(a) Escreva o modelo em termos de operador de retardo.

(b) Descreva como deve ter sido o comportamento das funções de autocorrela-

ção e autocorrelação parcial da série original e da série diferenciada.

(c) Como deve ser o comportamento esperado dos reśıduos em termos de suas

autocorrelações para que o modelo seja adequado?

(d) O que se deve fazer no caso de autocorrelações residuais significativas nas

defasagens 1, 8 e 12 ?

6. Para uma série observada trimestralmente foi identificado e estimado o modelo

SARIMA(1,1,0)×(2,1,1).

(a) Escreva o modelo em termos de operador de retardo.

(b) Descreva como deve ter sido o comportamento das funções de autocorrela-

ção e autocorrelação parcial da série original e da série diferenciada.

(c) O que se deve fazer se a autocorrelação residual na defasagem 4 for signi-

ficativa ?

7. Explique como você estimaria os coeficientes em um modelo ARMA(1,1) utili-

zando as duas primeiras autocorrelações amostrais?

8. Obtenha os estimadores de mı́nimos quadrados para os coeficientes em um mo-

delo AR(2).

9. Escreva as equações de mı́nimos quadrados para o modelo AR(p). Como você

estima a variância dos erros?

10. Em que condições as estimativas de mı́nimos quadrados de um modelo AR(p)

coincidirão com as de máxima verossimilhança?

11. Seja o modelo AR(1) com erros normais.

(a) Obtenha os EMV usando a verossimilhança condicional.

(b) Obtenha os EMV usando a verossimilhança exata com

X1 ∼ N(0, σ2ǫ /(1− α2)).

12. Usando as notas de aula e qualquer outra referência bibliográfica faça um resumo

da análise de reśıduos em séries temporais.

13. Explique como podem ser calculados os reśıduos em um modelo ARMA(p,q).

48 CAPÍTULO 4. ESTIMAÇÃO

> par(mfrow = c(3, 2))

> for (i in 1:3) {

+ plot(z[, i], main = "", xlab = "Anos", ylab = yl[i])

+ m = acf(z[, i], lag.max = 36, plot = F, na.action = na.pass)

+ m$lag = m$lag * 12

+ plot(m, main = "", xlab = "defasagem", ylab = "FAC")

+ }

Anos

N o

de p

as sa

ge iro

s

1950 1954 1958

5. 0

6. 0

0 5 10 15 20 25 30 35

− 0.

2 0.

4 0.

8

defasagem

F A

C

Anos

V ar

ia ca

o m

en sa

l

1950 1954 1958

− 0.

2 0.

0 0.

2

0 5 10 15 20 25 30 35

− 0.

2 0.

4 1.

0

defasagem

F A

C

Anos

V ar

ia ca

o an

ua l

1950 1954 1958

0. 0

0. 2

0 5 10 15 20 25 30 35

− 0.

2 0.

4 1.

0

defasagem

F A

C

Figura 4.1:

4.6. ADEQUAÇÃO DO MODELO 49

> par(mfrow = c(2, 1))

> plot(z, main = "serie com 1 diferenca simples e 1 sazonal", xlab = "Anos",

+ ylab = "")

> plot(m, main = "")

serie com 1 diferenca simples e 1 sazonal

Anos

1950 1952 1954 1956 1958 1960

− 0.

15 0.

00 0.

15

0 5 10 15 20 25 30 35

− 0.

4 0.

2 0.

8

Lag

A C

F

Figura 4.2:

50 CAPÍTULO 4. ESTIMAÇÃO

> tsdiag(m)

Standardized Residuals

Time

1950 1952 1954 1956 1958 1960

− 3

− 1

1 3

0.0 0.5 1.0 1.5

− 0.

2 0.

4 0.

8

Lag

A C

F

ACF of Residuals

2 4 6 8 10

0. 0

0. 4

0. 8

p values for Ljung−Box statistic

lag

p va

lu e

Figura 4.3:

4.6. ADEQUAÇÃO DO MODELO 51

> par(mfrow = c(2, 1))

> hist(z, freq = F)

> lines(d, dnorm(d, 0, sd(z)))

> qqnorm(z)

> qqline(z)

> text(-1.5, 0.05, "Teste de Shapiro-Wilk")

> text(-2, 0.01, paste("p-valor=", round(a$p.value, 4)))

Histogram of z

z

D en

si ty

−0.10 −0.05 0.00 0.05 0.10

0 4

8 12

−2 −1 0 1 2

− 0.

10 0.

05

Normal Q−Q Plot

Theoretical Quantiles

S am

pl e

Q ua

nt ile

s

Teste de Shapiro−Wilk p−valor= 0.1674

Figura 4.4:

Caṕıtulo 5

Previsão

Uma das formas de utilização de um modelo ajustado é para fazer previsões de valores

futuros. Assim, se t é o peŕıodo corrente estamos interessados em prever os valores

de Xt+1, Xt+2, . . . . A previsão de Xt+k, para k = 1, 2, . . . será denotada por x̂t(k) e é

definida como a esperança condicional de Xt+k dados todos os valores passados, i.e.

x̂t(k) = E(Xt+k|xt, xt−1, . . . ). (5.1)

A equação acima é chamada de função de previsão e o inteiro k é chamado de horizonte

de previsão. Pode-se mostrar que esta previsão tem o menor erro quadrático médio

(EQM), E(Xt+k − x̂t(k))2. Na prática temos um número finito de observações e obtemos então que x̂t(k) = E(Xt+k|xt, . . . , x1) que não tem o EQM mı́nimo mas pode ser visto como uma aproximação de (5.1).

Note que se temos uma série temporal observada x1, . . . , xn as previsões podem

ser feitas dentro do peŕıodo amostral e comparadas com os valores observados. Esta

é uma prática bastante comum para checar a performance preditiva do modelo. A

diferença entre os valores previsto e observado, x̂t(k) − xt+k, é chamada de erro de previsão k passos à frente e será denotado por et+k.

5.1 Métodos Univariados de Previsão

Os métodos descritos nesta seção têm um forte apelo intuitivo, decompondo uma série

temporal em componentes de fácil interpretação. Dados os recursos computacionais

dispońıveis atualmente eles também têm a vantagem de serem extremamente simples

de programar e sua utilização ter um custo computacional muito pequeno. Vamos

começar com o caso mais simples, adequado para séries localmente constantes.

5.1.1 Alisamento Exponencial Simples

Dada uma série temporal x1, . . . , xn, não sazonal e sem tendência sistemática, é razoá-

vel tomar a estimativa de xn+1 como uma soma ponderada das observações passadas,

i.e.

x̂n(1) = a0xn + a1xn−1 + . . .

52

5.1. MÉTODOS UNIVARIADOS DE PREVISÃO 53

onde {aj} são os pesos. Parece razoável também dar um peso maior às observações mais recentes do que às observações mais distantes no passado, i.e. a0 > a1 > . . . .

Neste procedimento são adotados pesos que decaem geometricamente a uma taxa

constante dados por

aj = α(1− α)j , j = 0, 1, . . .

onde 0 < α < 1 é chamada de constante de alisamento. Assim, a previsão 1 passo à

frente em t = n fica

x̂n(1) = αxn + α(1− α)xn−1 + α(1− α)2xn−2 + . . . . (5.2)

Naturalmente que na prática haverá um número finito de observações passadas e a

soma acima será também finita. A idéia de que o conteúdo informativo de uma obser-

vação decai com a sua “idade” é bastante intuitivo e o parâmetro α está controlando

o grau de “envelhecimento” deste conteúdo.

A equação (5.2) costuma ser reescrita em forma de equação recursiva. Colocando-

se (1− α) em evidência obtém-se que

x̂n(1) = αxn + (1− α)[αxn−1 + α(1− α)xn−2 + α(1− α)2xn−3 + . . . ] = αxn + (1− α)x̂n−1(1) (5.3)

i.e. uma média ponderada entre a observação mais recente e a previsão anterior.

A equação (5.2) pode ainda ser reescrita na forma de correção de erro. Definindo

en = xn − x̂n−1(1) o erro de previsão 1 passo à frente no tempo n então

x̂n(1) = x̂n−1(1) + αen.

Ou seja, a previsão para t = n + 1 é igual à previsão para t = n que foi feita em

t = n − 1 mais uma proporção do erro cometido. A previsão k-passos a frente é a mesma, i.e x̂n(k) = x̂n(1), k = 2, 3, . . . .

Previsões Dentro da Amostra

Usando x̂0(1) = x1 como previsão inicial em t = 0 e definindo et = xt − x̂t−1(1) os erros de previsão 1 passo à frente, a equação (5.3) pode ser usada recursivamente para

obter as previsões, i.e.

x̂t(1) = αxt + (1− α)x̂t−1(1), t = 1, 2, . . .

Na forma de correção de erro as recursões ficam

x̂t(1) = x̂t−1(1) + αet, t = 1, 2, . . .

Especificação de α

Vale notar que o valor de α não depende da escala em que as observações foram medi-

das, mas sim das propriedades da série temporal. O valor de α deve ser especificado de

modo a refletir a influência das observações passadas nas previsões. Valores pequenos

54 CAPÍTULO 5. PREVISÃO

produzem previsões que dependem de muitas observações passadas. Por outro lado,

valores próximos de 1 levam a previsões que dependem das observações mais recentes

e no caso extremo α = 1 a previsão é simplesmente a última observação. O valor de

α também pode ser estimado a partir dos dados e o critério utilizado é a minimização

da soma de quadrados dos erros de previsão. Ou seja, dado um valor fixo de α e

usando a equação (5.3), calcule

x̂0(1) = x1,

x̂1(1) = αx1 + (1− α)x̂0(1), e2 = x2 − x̂1(1) x̂2(1) = αx2 + (1− α)x̂1(1), e3 = x3 − x̂2(1)

...

x̂n−1(1) = αxn−1 + (1− α)x̂n−2(1), en = xn − x̂n−1(1)

e calcule ∑n

t=2 e 2 t . Repita o procedimento para valores de α variando entre 0 e 1

(digamos com incrementos de 0,1) e selecione o valor que minimiza esta soma de qua-

drados. Na prática, o valor mı́nimo pode ocorrer muito próximo de um dos extremos

do intervalo de variação de α. Isto pode ocorrer quando a soma de quadrados varia

muito pouco na região em torno do mı́nimo. Neste caso faz mais sentido utilizar

valores não tão extremos.

Exemplo 5.1 : No banco de dados do R, a série lh contém as quantidades de um

tipo de hormônio em amostras de sangue coletadas a cada 10 minutos de uma pessoa

do sexo feminino (Diggle 1990). Vamos aplicar o método de alisamento exponencial

simples à esta série fazendo primeiro a seleção do valor de α que minimiza a soma

dos quadrados dos erros de previsão 1 passo a frente. Na Figura 5.1 temos o gráfico

desta soma de quadrados como função de α e o gráfico das previsões 1 passo à frente

juntamente com a série observada.

O valor ótimo obtido foi α = 0, 945 com a soma de erros quadrados igual a 11,86

e os seguintes comandos do R podem ser utilizados para a seleção de α.

> AES = function(x, interval) {

+ e = NULL

+ for (alfa in interval) {

+ e2 = 0

+ prev = x[1]

+ for (i in 2:length(x)) {

+ prev = c(prev, alfa * x[i - 1] + (1 - alfa) * prev[i -

+ 1])

+ e2 = e2 + (x[i] - prev[i])^2

+ }

+ e = c(e, e2)

+ }

+ plot(interval, e, type = "l", xlab = expression(alpha), ylab = "Soma de quadrados

+ e.min = min(e)

5.1. MÉTODOS UNIVARIADOS DE PREVISÃO 55

+ alfa = interval[e == e.min]

+ prev = x[1]

+ for (i in 2:length(x)) prev = c(prev, alfa * x[i - 1] + (1 -

+ alfa) * prev[i - 1])

+ return(list(alfa = alfa, sq2 = e.min, prev = prev))

+ }

> par(mfrow = c(2, 1))

> m = AES(lh, seq(0.1, 0.99, 0.001))

> plot(1:48, m$prev, ylab = "Hormonio", xlab = "Amostras", type = "l")

> points(lh)

0.2 0.4 0.6 0.8 1.0

12 .0

13 .5

α

S om

a de

q ua

dr ad

os d

os e

rr os

0 10 20 30 40

1. 5

2. 5

3. 5

Amostras

H or

m on

io

Figura 5.1: Soma dos quadrados dos erros de previsão 1 passo a frente em função de α. Valores observados (pontos) e previsões 1 passo a frente (linhas) usando o valor ótimo de α.

Exemplo 5.2 : O procedimento do Exemplo 5.1 foi repetido para a série de medidas

anuais de vazões do Rio Nilo entre 1871 e 1970, também do banco de dados do R. Os

resultados estão na Figura 5.2.

56 CAPÍTULO 5. PREVISÃO

> par(mfrow = c(2, 1))

> m = AES(Nile, seq(0.1, 0.99, 0.001))

> plot(1:length(Nile), m$prev, ylab = "", xlab = "", type = "l")

> points(1:length(Nile), Nile)

0.2 0.4 0.6 0.8 1.0

21 00

00 0

27 00

00 0

α

S om

a de

q ua

dr ad

os d

os e

rr os

0 20 40 60 80 100

80 0

10 00

Figura 5.2: Soma dos quadrados dos erros de previsão 1 passo a frente em função de α. Valores observados (pontos) e previsões 1 passo a frente (linhas) usando o valor ótimo de α

5.1.2 Método de Holt-Winters

O procedimento de alisamento exponencial pode ser generalizado para séries que con-

tenham tendência e variação sazonal. Suponha por exemplo que as observações são

mensais e sejam Lt, Tt e It o ńıvel, a tendência e o ı́ndice sazonal no tempo t. Assim,

Tt é o aumento ou redução esperada por mês no ńıvel atual da série.

Suponha que no tempo t os termos (L1, T1, I1), . . . , (Lt−1, Tt−1, It−1) sejam conhe-

cidos. Então, após observar xt os termos Lt, Tt e It são atualizados via alisamento

exponencial. Se a variação sazonal for multiplicativa, i.e. com amplitudes que tendem

a crescer ao longo do tempo, as equações de atualização na forma de recorrência são

5.1. MÉTODOS UNIVARIADOS DE PREVISÃO 57

dadas por

Lt = α(xt/It−12) + (1− α)(Lt−1 + Tt−1) Tt = γ(Lt − Lt−1) + (1− γ)Tt−1 It = δ(xt/Lt) + (1− δ)It−12

e as previsões k peŕıodos à frente são dadas por

x̂t(k) = (Lt + kTt)It−12+k, k = 1, 2, . . . .

No caso de sazonalidade aditiva as equações de atualização para o ńıvel e o ı́ndice

sazonal são modificadas para

Lt = α(xt − It−12) + (1− α)(Lt−1 + Tt−1) It = δ(xt − Lt) + (1− δ)It−12

e as previsões k peŕıodos à frente ficam

x̂t(k) = Lt + kTt + It−12+k, k = 1, 2, . . . .

Aqui temos parâmetros de alisamento, α, γ e δ, para cada componente da série

que são em geral escolhidos no intervalo (0,1) e podem ser estimados minimizando-se

a soma de quadrados dos erros de previsão como na seção anterior. Aqui vale tam-

bém o comentário sobre valores próximos aos extremos devido à soma de quadrados

variar pouco nesta região. Além disso, estes parâmetros não dependem da escala das

observações mas sim das propriedades temporais do ńıvel, tendência e sazonalidade

da série. Valem os mesmos comentários sobre estes valores refletindo a influência das

observações passadas nas previsões de cada componente.

Para o caso particular de séries sem variação sazonal basta utilizar as equações

para Lt e Tt acima (sem o ı́ndice It−12). Ou seja,

Lt = αxt + (1− α)(Lt−1 + Tt−1) Tt = γ(Lt − Lt−1) + (1− γ)Tt−1

e a previsão k passos à frente no tempo t é simplesmente Lt+kTt. Se a série também

não tem uma tendência sistemática retorna-se à equação (5.3), ou seja

Lt = αxt + (1− α)Lt−1

e Lt é a previsão 1 passo à frente (x̂t(1)).

Exemplo 5.3 : A variável UKLungDeaths contém os números mensais de mortes por

doenças do pulmão (bronquite, efisema e asma) no Reino Unido entre janeiro de 1974

e dezembro de 1979. A variável é composta por 3 séries: ambos os sexos (ldeaths),

sexo feminino (fdeaths) e sexo masculino (mdeaths). Aqui será utilizada a função

HoltWinters do R que faz o alisamento exponencial de Holt-Winters com a série lde-

aths. As constantes de alisamento (α , β e γ ) são determinadas minimizando a soma

58 CAPÍTULO 5. PREVISÃO

dos quadrados dos erro de previsão 1 passo à frente. Considere um modelo sazonal

aditivo. O resultado são as constantes de alisamento calculadas e as Estimativas fi-

nais (em t = n) do nivel, tendência e componentes sazonais. Pode-se também obter

as previsões e intervalos de previsão (supondo normalidade) para modelos ajustados

pelo método de Holt-Winters. No gráfico da Figura 5.3 temos a série original com a

série suavizada mais as previsões para os anos de 1980, 1981 e 1982 da série ldeaths.

> data(UKLungDeaths)

> m = HoltWinters(ldeaths, seasonal = "addit")

> p = predict(m, n.ahead = 12, prediction.interval = T)

> plot(m, p)

Holt−Winters filtering

Time

O bs

er ve

d / F

itt ed

1975 1976 1977 1978 1979 1980 1981

50 0

10 00

15 00

20 00

25 00

30 00

35 00

40 00

Figura 5.3: Série original, série suavizada e previsões para o ano de 1980 da série ldeaths via método de Holt-Winters.

5.2 Previsão em Modelos ARMA

Em modelos ARMA as previsões podem ser obtidas usando-se diretamente a equação

do modelo. Assumindo que a equação do modelo seja conhecida a previsão x̂n(k)

é obtida substituido valores futuros dos erros ǫ por zero, valores futuros da série

5.2. PREVISÃO EM MODELOS ARMA 59

Xn+1, Xn+2, . . . pela sua esperança condicional, e valores passados de X e de ǫ pelos

seus valores observados.

Tomemos como exemplo o modelo SARIMA(1, 0, 0) × (0, 1, 1)12. A equação do modelo é dada por

(1− αB)(1−B12)Xt = (1 + θB12)ǫt

ou equivalentemente

Xt = Xt−12 + α(Xt−1 −Xt−13) + ǫt + θǫt−12.

Neste caso, as previsões 1 e 2 passos à frente ficam

x̂n(1) = xn−11 + α(xn − xn−12) + θǫn−11 x̂n(2) = xn−10 + α(x̂n(1)− xn−11) + θǫn−10.

Note como o valor futuro Xn+1 foi substitúıdo na segunda equação pela sua esperança

condicional x̂n(1), i.e. a previsão feita em t = n para t = n + 1. Previsões para

horizontes maiores podem ser obtidas recursivamente.

No caso de modelos autoregressivos AR(p) não é dif́ıcil verificar como fica a função

de previsão.

x̂t(1) = α1xt + · · ·+ αpxt−p+1 x̂t(2) = α1x̂t(1) + · · ·+ αpxt−p+2

...

x̂t(p+ 1) = α1x̂t(p) + · · ·+ αpx̂t(1)

de modo que as previsões para horizontes maiores do que p usam apenas as previsões

anteriores. Para p = 1 por exemplo segue que

x̂t(k) = αx̂t(k − 1) = α2x̂t(k − 2) = · · · = αkxt

Para modelos médias móveis MA(q) também não é dif́ıcil verificar que a equação

de previsão fica

x̂t(1) = β1ǫt + · · ·+ βqǫt−q+1 x̂t(2) = β2ǫt + · · ·+ βqǫt−q+2

...

x̂t(q) = βqǫt

x̂t(q + j) = 0, j = 1, 2, . . .

ou seja,

x̂t(k) =

{

∑q i=k βiǫt+k−i, k = 1, . . . , q

0, k > q

60 CAPÍTULO 5. PREVISÃO

Atualização das Previsões

É interessante notar também como as previsões podem ser atualizadas conforme novas

observações da série forem obtidas. Suponha por exemplo que o valor xn+1 foi obser-

vado. Neste caso a previsão para t = n + 2 ficará condicionada em x1, . . . , xn, xn+1 e pode ser facilmente atualizada para a nova origem n+ 1. Para o modelo SARIMA

visto acima a previsão fica

x̂n+1(1) = E(Xn+2|xn+1, . . . , x1) = xn−10 + α(xn+1 − xn−11) + θǫn−10. (5.4)

Somando e subtraindo αx̂n(1) no lado direito de (5.4) obtemos que

x̂n+1(1) = xn−10 + α(x̂n(1)− xn−11) + α(xn+1 − x̂n(1)) + θǫn−10 = x̂n(2) + α(xn+1 − x̂n(1))

ou seja, a previsão atualizada é a previsão anterior mais uma proporção do erro de

previsão 1 passo à frente em t = n+ 1.

Previsões usando a forma MA

As previsões também podem ser obtidas reescrevendo-se o modelo como um processo

médias móveis de ordem infinita. Neste caso temos que

Xn+k = ǫn+k + ψ1ǫn+k−1 + · · ·+ ψkǫn + ψk+1ǫn−1 + . . .

e fica claro que a previsão k passos à frente é dada por

x̂n(k) = ψkǫn + ψk+1ǫn−1 + . . . . (5.5)

Note que apenas os valores ǫn, ǫn−1, . . . foram utilizados já que a esperança dos valores

futuros é igual a zero. Esta forma é particularmente útil para o cálculo da variância

do erro de previsão. Da equação (5.5) obtemos que o erro de previsão k passos à

frente é dado por

xn+k − x̂n(k) = ǫn+k + ψ1ǫn+k−1 + · · ·+ ψk−1ǫn+1

e portanto a variância do erro de previsão é dada por

V ar(et+k) = (1 + ψ 2 1 + · · ·+ ψ2k−1)σ2ǫ .

O ponto importante a se notar aqui é que, para σ2ǫ fixo, a variância do erro de previsão

aumenta com o horizonte de previsão. Na prática, isto significa ter mais confiança

em previsões de curto prazo do que em previsões de longo prazo.

Até agora não haviamos assumido nenhuma distribuição de probabilidade para os

erros. Assumindo também que a sequência {ǫt} seja normalmente distribuida pode-se

5.2. PREVISÃO EM MODELOS ARMA 61

construir intervalos de confiança para Xt+k simétricos em torno das previsões. Estes

são chamados intervalos de previsão e são dados por

x̂t(k)± zα/2

1 + k−1 ∑

j=1

ψ2j

σ2ǫ .

É claro que neste caso a hipótese de normalidade precisa ser checada.

Finalmente, vale notar que na prática os parâmetros do modelo não são conhecidos

de forma exata e precisam ser estimados. Os valores passados dos erros ǫt também

precisam ser estimados como erros de previsão um passo à frente. Assim, por exemplo

para o modelo SARIMA(1, 0, 0)× (0, 1, 1)12 visto acima teremos que

x̂n(1) = xn−11 + α̂(xn − xn−12) + θ̂ǫ̂n−11

onde o erro de previsão 1 passo à frente em n− 11 é dado por

ǫ̂n−11 = xn−11 − x̂n−12(1).

Além disso, os intervalos de previsão obtidos serão intervalos aproximados devido a

esta substituição.

Exemplo 5.4 : A Figura 5.4 mostra uma série temporal com os totais mensais de

mortes por acidente nos Estados Unidos entre janeiro de 1973 e dezembro de 1978.

Suponha que foi identificado o modelo SARIMA(0,1,1)x(0,1,1). Após a estimação,

análise de reśıduos e verificação da adequação do modelo foram feitas previsões para

o ano de 1979, i.e. previsões 1, 2, . . . , 12 passos à frente. Em julho de 1979 os valores

para os primeiros 6 meses daquele ano foram disponibilizados e aparecem na Figura

5.5 juntamente com as previsões. Note como os valores observados ficaram dentro

dos intervalos de previsão fornecendo assim indicação de que o modelo teve uma boa

performance preditiva. Sendo assim, uma estratégia inicial para o segundo semestre

de 1979 consiste em simplesmente atualizar as previsões. Os comandos do R para este

exemplo são dados a seguir.

Transformações

Em muitas aplicações a série modelada é na verdade uma transformação dos dados

originais, sendo a transformação logaritmica a mais usual. Assim, tanto as previsões

pontuais quanto os intervalos de previsão são obtidos para a série transformada e estes

valores precisam ser transformados novamente para a escala original. A abordagem

mais simples (e geralmente adotada) consiste simplesmente em tomar a transformação

inversa, por exemplo se um modelo foi ajustado para a série Xt = log Yt então ŷn(k) =

exp(x̂n(k)) é a previsão k passos a frente da série original. No entanto deve-se ter em

mente que estas previsões via transformação inversa são em geral viesadas. Felismente

os intervalos de previsão tem boas propriedades e por exemplo quanto se toma o anti-

logaritmo dos limites

x̂n(k)± zα/2 √

var(en+k)

62 CAPÍTULO 5. PREVISÃO

> data(USAccDeaths)

> plot(USAccDeaths, xlab = "Anos", ylab = "Numero de mortes por acidente")

Anos

N um

er o

de m

or te

s po

r ac

id en

te

1973 1974 1975 1976 1977 1978 1979

70 00

80 00

90 00

10 00

0 11

00 0

Figura 5.4: Totais mensais de mortes por acidente nos Estados Unidos entre janeiro de 1973 e dezembro de 1978.

obtém-se um intervalo (geralmente assimétrico) de 100(1−α)% para a previsão pon- tual ŷn(k).

Exemplo 5.5 : Considere novamente a série AirPassengers e faça transformação

logaritmica nos dados (conforme sugerido na literatura). Estime um modelo SA-

RIMA(0,1,1)x(0,1,1) usando os dados até dezembro de 1960 e faça previsões de 1 até

12 meses à frente para o ano de 1961 nas 2 escalas. As previsões e intervalos de previ-

são na escala transformada são dados na Tabela 5.1, enquanto as previsões, intervalos

de previsão e suas semi-amplitudes na escala original são dadas na Tabela 5.2.

5.3 Performance Preditiva

A idéia de verificar a adequação de um modelo em termos dos erros de previsão um

passo à frente foi apresentada na Seção 4.6. Na prática é preciso verificar se os reśıduos

se comportam de maneira aleatória (ou impreviśıvel) em torno de zero e com variância

5.3. PERFORMANCE PREDITIVA 63

previsão li ls

1961 Jan 6.11 6.04 6.18

1961 Feb 6.05 5.97 6.14

1961 Mar 6.17 6.08 6.27

1961 Apr 6.20 6.09 6.31

1961 May 6.23 6.12 6.35

1961 Jun 6.37 6.25 6.49

1961 Jul 6.51 6.38 6.64

1961 Aug 6.50 6.37 6.64

1961 Sep 6.32 6.18 6.47

1961 Oct 6.21 6.06 6.36

1961 Nov 6.06 5.91 6.22

1961 Dec 6.17 6.00 6.33

Tabela 5.1: Previsões e limites inferior (li) e superior (ls) dos intervalos de previsão.

aproximadamente constante, além de serem não correlacionados ao longo do tempo.

Além disso, dois ou mais modelos podem ser comparados segundo a sua perfor-

mance preditiva, ou seja construindo-se medidas baseadas nos erros de previsão. A

maioria dos métodos de previsão baseia-se na idéia de minimizar somas de quadrados

ou de valores absolutos dos erros de previsão e esta é também uma medida usada

para comparar a adequação de modelos alternativos. A idéia então é comparar o erro

quadrático médio ∑

e2t /(n−m) ou erro absoluto médio ∑ |et|/(n−m) para diferentes

modelos, onde m é o número de parâmetros a serem estimados.

Uma estratégia simples de se fazer previsões consiste em tomar a observação mais

recente como a melhor previsão de um valor futuro da série, i.e. x̂t(1) = xt. Note

que esta é a previsão 1 passo à frente de um passeio aleatório. Assim, uma forma

de medir a capacidade preditiva de um modelo consiste em comparar seus erros de

previsão com aqueles do passeio aleatório. Isto pode ser feito através da chamada

estat́ıstica U de Theil definida como

U =

∑n−1 t=1 (xt+1 − x̂t(1))2 ∑n−1

t=1 (xt+1 − xt)2 .

Note que valores maiores do que 1 são uma indicação de que globalmente os erros

de previsão tendem a ser grandes em relação aos erros de um passeio aleatório. Esta

não é uma boa caracteŕıstica e gostariamos que o valor de U fosse sempre menor do

que 1. Vale notar também que neste caso os erros de previsão estão sendo avaliados

independente da escala dos dados.

Finalmente, vale notar que todas as medidas de capacidade preditiva citadas po-

dem ser estendidas para erros de previsão k passos a frente.

Outra prática comum em séries temporais consiste em estimar o modelo excluindo

algumas observações finais e usar o modelo estimado para fazer previsões. Neste caso

as previsões podem ser comparadas com os valores observados. Por exemplo, para uma

64 CAPÍTULO 5. PREVISÃO

prev li ls prev.li ls.prev

1961 Jan 450.42 418.53 484.74 31.89 34.32

1961 Feb 425.72 390.81 463.75 34.91 38.03

1961 Mar 479.01 435.08 527.37 43.93 48.36

1961 Apr 492.40 443.00 547.32 49.41 54.92

1961 May 509.05 453.98 570.81 55.07 61.75

1961 Jun 583.34 516.02 659.45 67.33 76.11

1961 Jul 670.01 588.18 763.23 81.83 93.22

1961 Aug 667.08 581.40 765.38 85.68 98.30

1961 Sep 558.19 483.18 644.85 75.01 86.66

1961 Oct 497.21 427.59 578.17 69.62 80.96

1961 Nov 429.87 367.37 503.01 62.50 73.14

1961 Dec 477.24 405.40 561.81 71.84 84.57

Tabela 5.2: Previsões e limites inferior (li) e superior (ls) e semi-amplitudes dos

intervalos de previsão.

série mensal observada ao longo de 5 anos poderia-se estimar o modelo identificado

usando os primeiros 4 anos e meio (54 observações) e fazer previsões para os últimos

6 meses.

5.4 Critérios de Informação

Em muitas aplicações vários modelos podem ser julgados adequados em termos do

comportamento dos reśıduos. Uma forma de “discriminar” entre estes modelos com-

petidores é utilizar os chamados critérios de informação que levam em conta não

apenas a qualidade do ajuste mas também penalizam a inclusão de parâmetros ex-

tras. Assim, um modelo com mais parâmetros pode ter um ajuste melhor mas não

necessariamente será prefeŕıvel em termos de critério de informação. A regra básica

consiste em selecionar o modelo cujo critério de informação calculado seja mı́nimo.

A regra mais utilizada em séries temporais é o chamado critério de informação de

Akaike, denotado por AIC. A definição mais comumente utilizada é

AIC = −2 log verossimilhança maximizada + 2m1

onde m é o número de parâmetros (em modelos ARMA(p, q) m = p + q + 1). Para

dados normalmente distribuidos e usando-se estimativas de máxima verossimilhança

para os parâmetros pode-se mostrar que

AIC = n log(σ̂2ǫ ) + 2m

onde σ̂2ǫ = (1/n) ∑

ǫ̂2t .

1O fator 2 é somente uma convenção e não irá alterar a seleção do modelo.

5.4. CRITÉRIOS DE INFORMAÇÃO 65

Existem outros critérios de informação que são basicamente modificações do AIC

na forma de penalizar a inclusão de parâmetros extras. O mais famoso deles é o

critério de informação Bayesiano, denotado por BIC e dado por

BIC = −2 log verossimilhança maximizada +m log n.

Note como este critério penaliza bem mais a inclusão de parâmetros do que o AIC e

portanto tende a selecionar modelos mais parcimoniosos.

É sempre bom lembrar que estas medidas não têm nenhum significado quando

olhadas individualmente, i.e. considerando-se um único modelo. Assim, tanto o AIC

quanto o BIC podem assumir valores quaisquer, inclusive valores negativos, já que

eles dependem da forma da função de verossimilhança.

Vale lembrar também que ao usar tais critérios para comparar modelos a esti-

mação precisa ser feita no mesmo peŕıodo amostral de modo que os modelos sejam

comparáveis. Note também que aumentando-se o número de termos autoregressivos

e/ou médias móveis, o valor de m aumenta. Assim se a inclusão de termos adicionais

no modelo não melhorar sensivelmente o ajuste, então o AIC e o BIC (e qualquer

outro critério de informação) serão maiores.

Para uma revisão geral destes e outros critérios de informação no contexto de

séries temporais ver por exemplo Priestley (1981), Caṕıtulo 5.

Identificação Revisitada

Vimos que as duas ferramentas básicas para identificação de modelos da classe ARIMA

são as autocorrelações e autocorrelações parciais amostrais. Esta etapa envolve al-

gum grau de arbitrariedade por parte do pesquisador ao interpretar estas funções,

i.e. comparar subjetivamente seus valores amostrais com os correspondentes valores

teóricos.

Uma abordagem alternativa consiste em usar os critérios de informação de um

forma mais abrangente. Neste caso, um conjunto de posśıveis modelos competidores

é definido a priori e aquele que minimiza o AIC ou BIC é selecionado. Por exemplo,

modelos ARMA(p, q) podem ser estimados sequencialmente variando os valores de p

e q entre 0 e 3 digamos. Note que neste caso teremos 16 posśıveis modelos sendo

comparados e os critérios de informação são agora funções de p e q. Analogamente,

para modelos AR(p) podemos variar o valor de p, digamos entre 1 e 10.

Na prática este procedimento pode ser aplicado de forma semi-automática já que

muitos pacotes estat́ısticos fornecem estes valores. Porém após um modelo ser seleci-

onado a análise residual ainda deve ser feita antes de se passar à etapa das previsões.

Outro problema de ordem prática é que pode haver dois ou mais modelos com AIC

e/ou BIC muito similares de modo que não seja trivial discriminar entre eles. Nestas

situações Burnham & Anderson (1998), Seção 4.2, sugerem o uso de pesos que são

obtidos subtraindo-se o valor associado com o “melhor”modelo. Os pesos relativos ao

AIC são dados por

wk ∝ exp(−∆AIC(k)/2)

66 CAPÍTULO 5. PREVISÃO

sendo ∆AIC(k) = AIC(k) − min(AIC) e k é a ordem do modelo. Estes pesos são então normalizados para somarem 1 de modo que 0 < wk < 1 e a comparação entre

os modelos fica mais fácil. Se M é o número total de modelos a comparação é então

baseada em

w∗i = wi

∑M j=1 wj

, i = 1, . . . ,M.

Por exemplo, para modelos AR(p) os pesos relativos ao AIC são dados por

wp ∝ exp(−∆AIC(p)/2), p = 1, . . . , pmax

sendo ∆AIC(p) = AIC(p)−min(AIC) e pmax deve ser especificado.

Exemplo 5.6 : Na Figura 5.6 é apresentada a série com os totais anuais de linces

canadenses capturados em armadilhas entre 1821 e 1934. Estes dados têm sido mo-

delados na literatura após uma transformação que consiste em tomar o logaritmo na

base 10 e subtrair a média dos dados transformados. Vamos ajustar modelos AR(p)

com p variando de 1 até 5 e calcular os critérios de informação e os respectivos pesos

para cada modelo. Os resultados estão na Tabela 5.3. Note que há falta de concor-

dância entre os critérios de informação quanto ao melhor modelo. Isto pode ser uma

indicação de que na verdade há 2 modelos descrevendo bem os dados. Outro pro-

blema é que o AIC seleciona um modelo com o valor máximo de p e isto pode indicar

a necessidade de considerar mais termos autoregressivos. Repetindo o exercicio com

p variando de 1 a 15 obteve-se a Tabela 5.4.

p AIC pesos AIC BIC pesos BIC

1 1 −242.3913 0.0000 −234.9189 0.0000 2 2 −333.0988 0.1057 −321.8902 0.8137 3 3 −332.7283 0.0878 −317.7835 0.1044 4 4 −335.6596 0.3802 −316.9786 0.0698 5 5 −335.8881 0.4263 −313.4709 0.0121

Tabela 5.3: Critérios de informação AIC e BIC e respectivos pesos para modelos

AR(p) ajustados a série Lynx.

Os comandos do R utilizados no Exemplo 5.6 seguem abaixo.

> y = log10(lynx)

> x = y - mean(y)

> p = 1:15

> n = length(x)

> crit = matrix(0, nrow = length(p), ncol = 5)

> colnames(crit) = c("p", "AIC", "pesos AIC", "BIC", "pesos BIC")

> crit[, 1] = p

> for (k in p) {

+ ar = arima(x, order = c(k, 0, 0), include.mean = F)

5.5. PREVISÕES USANDO TODOS OS MODELOS 67

+ crit[k, 2] = n * log(ar$sigma2) + 2 * (k + 1)

+ crit[k, 4] = n * log(ar$sigma2) + (k + 1) + (k + 1) * log(n)

+ }

> aicp = exp(-(crit[, 2] - min(crit[, 2]))/2)

> bicp = exp(-(crit[, 4] - min(crit[, 4]))/2)

> crit[, 3] = aicp/sum(aicp)

> crit[, 5] = bicp/sum(bicp)

p AIC pesos AIC BIC pesos BIC

1 1 −242.3913 0.0000 −234.9189 0.0000 2 2 −333.0988 0.0000 −321.8902 0.8100 3 3 −332.7283 0.0000 −317.7835 0.1039 4 4 −335.6596 0.0000 −316.9786 0.0695 5 5 −335.8881 0.0000 −313.4709 0.0120 6 6 −334.4484 0.0000 −308.2950 0.0009 7 7 −338.8427 0.0001 −308.9531 0.0013 8 8 −338.8505 0.0001 −305.2247 0.0002 9 9 −338.3849 0.0001 −301.0229 0.0000 10 10 −341.8678 0.0006 −300.7696 0.0000 11 11 −354.5690 0.3581 −309.7346 0.0019 12 12 −354.7117 0.3846 −306.1411 0.0003 13 13 −353.0609 0.1685 −300.7541 0.0000 14 14 −351.0895 0.0629 −295.0465 0.0000 15 15 −349.2335 0.0249 −289.4543 0.0000

Tabela 5.4: Critérios de informação AIC e BIC e respectivos pesos para modelos

AR(p) ajustados a série Lynx.

Finalmente vale notar que se o número de modelos candidatos for muito grande

e a série analisada muito longa o custo computacional deste método pode ser muito

alto. Por exemplo, em modelos SARIMA com pmax = qmax = 5, Pmax = Qmax = 2

e dmax = Dmax = 2 teremos mais de 500 modelos candidatos, sem contar posśıveis

transformações nos dados, diferentes distribuições dos erros, presença de dados dis-

crepantes, alterações estruturais, etc.

5.5 Previsões Usando Todos os Modelos

Suponha que existem k modelos“candidatos”denotados por M1,M2, . . . ,Mk e deseja-

se fazer a previsão de Xn+h. Tratando tanto Xn+h quanto Mi como variáveis aleató-

rias então pelas regras de esperança condicional segue que

x̂n(h) = E(Xn+h|x) = k

i=1

E(Xn+h|x,Mi)P (Mi|x).

68 CAPÍTULO 5. PREVISÃO

Ou seja, podemos escrever a previsão pontual como uma mistura discreta de previsões

pontuais sob cada modelo considerado.

A mesma lógica se aplica a qualquer função de Xn+h, em particular

E(X2n+h|x) = k

i=1

E(X2n+h|x,Mi)P (Mi|x).

que pode ser usado para quantificar a incerteza sobre Xn+h, i.e.

V ar(Xn+h|x) = E(X2n+h|x)− [E(Xn+h|x)]2

= k

i=1

E(X2n+h|x,Mi)P (Mi|x)− [E(Xn+h|x)]2

= k

i=1

[V ar(Xn+h|x,Mi) + E2(Xn+h|x,Mi)]P (Mi|x)− [x̂n(h)]2

Um procedimento para fazer previsões usando todos os modelos estimados consiste

em substituir as probabilidades P (Mi|x) pelos pesos wi padronizados. Por exemplo, para modelos autoregressivos se pmax é o número máximo de defasagens então

E(Xn+h|x) = pmax ∑

i=1

E(Xn+h|x, AR(i))w∗i .

5.6 Previsão Bayesiana

Na prática, os métodos de previsão em modelos ARIMA são aplicados substituindo-

se os parâmetros do modelo pelas suas estimativas pontuais. Porém o fato de não

conhecermos os valores dos parâmetros é mais uma fonte de incerteza em relação as

previsões e que em muitas situações pode ser muito grande para ser ignorada.

No contexto Bayesiano esta incerteza pode ser levada em conta já que a previsão

de valores futuros é feita a partir da distribuição preditiva de Xn+h, que é dada por

p(xn+h|x) = ∫

p(xn+h|x,θ)p(θ|x)dθ.

Neste caso, todos os posśıveis valores de θ estão sendo levados em conta e não apenas

a sua estimativa pontual.

Exerćıcios

1. No alisamento exponencial simples descreva a papel do parâmetro α.

2. No método de Holt-Winters descreva o papel dos parâmetros α, γ e δ.

3. Explique em que situações seriam usados os métodos de Holt-Winters aditivo

ou multiplicativo.

5.6. PREVISÃO BAYESIANA 69

4. Seja o modelo MA(1), Xt = ǫt + θǫt−1.

(a) Obtenha a previsão 1 passo à frente em t = n e mostre que as previsões k

passos à frente para k = 2, 3, . . . são iguais a zero.

(b) Mostre que a variância do erro de previsão k passos à frente é dada por σ2ǫ para k = 1 e (1 + θ2)σ2ǫ para k = 2, 3, . . . .

5. Seja o modelo Xt = 90 + ǫt + 0, 8ǫt−1 + 0, 5ǫt−1.

(a) Obtenha as previsões k passos à frente em t = n.

(b) Obtenha a variância do erro de previsão k passos à frente.

6. Seja o modelo AR(1), Xt = αXt−1 + ǫt.

(a) Mostre que a previsão k passos à frente feita em t = n é dada por αkxn.

(b) Mostre que a variância do erro de previsão k passos à frente é dada por

σ2ǫ (1− α2k)/(1− α2).

7. Para o modelo SARIMA(0, 0, 1)×(1, 1, 0)12 obtenha as previsões no tempo t = n para até 12 peŕıodos à frente em termos das observações e residuos até o tempo

t = n.

8. Seja o modelo (1−B)(1− 0, 2B)Xt = (1− 0, 5B)ǫt.

(a) Obtenha as previsões 1 e 2 passos à frente.

(b) Mostre que as previsões 3 ou mais passos à frente são dadas pela equação

recursiva x̂n(k) = 1, 2x̂n(k − 1)− 0, 2x̂n(k − 2). (c) Obtenha a variância dos erros de previsão 1, 2 e 3 passos à frente.

(d) Obtenha a previsão x̂n(2) e o erro padrão do erro de previsão sabendo que

ǫn = 1, xn = 4, xn−1 = 3 e σ 2 ǫ = 2.

9. Seja o modelo ARIMA(1,0,1) para uma série Xt com média zero.

(a) Reescreva o modelo na forma de choques aleatórios, i.e.

Xt = ǫt + ψ1ǫt−1 + ψ2ǫt−2 + . . .

obtendo uma expressão geral para os coeficientes ψj .

(b) Escreva a expressão da variância do erro de previsão et(k) = xt+k − x̂t(k). (c) Obtenha as previsões x̂t(k) para horizontes k = 1 e k > 1.

10. Sabe-se que se Y ∼ N(µ, σ2) então X = exp(Y ) tem distribuição log-normal com E(X) = exp(µ + σ2/2) e V ar(X) = e2µ+σ

2

(eσ 2 − 1). Se foram obtidas as

previsões k passos à frente de Yt = log(Xt) e assumindo que Yt é normal mostre

que as previsões na escala original são dadas por

X̂t(k) = exp(Ŷt(k) + Vy(k)/2)

com variância

exp(2Ŷt(k) + Vy(k)) [exp(Vy(k))− 1].

70 CAPÍTULO 5. PREVISÃO

11. Deseja-se ajustar um modelo ARMA a uma série temporal estacionária mas

os gráficos das funções de autocorrelação e autocorrelação parcial são pouco

informativos. Descreva um procedimento de identificação alternativo (você tem

um pacote estat́ıstico para fazer as contas).

12. Descreva um procedimento para obter previsões h passos à frente em modelos

autoregressivos com número máximo de defasagens igual a kmax utilizando todos

os modelos estimados. Ilustre situações em que as previsões pontuais médias

devem muito similares (ou muito diferentes) das previsões usando somente o

melhor modelo.

5.6. PREVISÃO BAYESIANA 71

> plot(ts(c(USAccDeaths, pacc$pred), frequency = 12, start = c(1973,

+ 1)), xlab = "Anos", ylab = "Numero de mortes por acidente",

+ ylim = c(6000, 12000))

> abline(v = 1979 - 1/12, lty = 2)

> lines(pacc$pred + 1.96 * pacc$se, lty = 2)

> lines(pacc$pred - 1.96 * pacc$se, lty = 2)

> obs79 = c(7798, 7406, 8363, 8460, 9217, 9316)

> points(1979 + (0:5)/12, obs79, pch = "*")

Anos

N um

er o

de m

or te

s po

r ac

id en

te

1973 1974 1975 1976 1977 1978 1979 1980

60 00

70 00

80 00

90 00

10 00

0 11

00 0

12 00

0

*

*

**

**

Figura 5.5: Previsões para 1979 com observações do primeiro semestre incluidas.

72 CAPÍTULO 5. PREVISÃO

Time

ly nx

1820 1840 1860 1880 1900 1920

0 10

00 20

00 30

00 40

00 50

00 60

00 70

00

Figura 5.6: Totais anuais de linces canadenses capturados em armadilhas entre 1821 e 1934.

Caṕıtulo 6

Modelando a Variância

6.1 Introdução

Nos modelos vistos até aqui a variância dos erros foi assumida constante ao longo

do tempo, i.e. V ar(ǫt) = E(ǫ 2 t ) = σ

2 ǫ . Muitas séries temporais no entanto exibem

peŕıodos de grande volatilidade seguidos de peŕıodos de relativa tranquilidade. Nestes

casos, a suposição de variância constante (homocedasticidade) pode não ser apropri-

ada. Na verdade, embora a variância incondicional dos erros ainda possa ser assumida

constante, sua variância condicional pode estar mudando ao longo do tempo.

Além disso, em muitas situações práticas tem-se interesse em prever a variância

condicional da série além da série propriamente dita. Por exemplo, no mercado de

ações o interesse é não apenas prever a taxa de retorno mas também a sua variância

ao longo de um certo peŕıodo. Esta variância condicional é também chamada de

volatilidade. Algumas referências para este caṕıtulo são Taylor (1986), Franses (1998),

e Tsay (2002).

Exemplo 6.1 : Na Figura 6.1 os gráficos da esquerda apresentam as taxas de câm-

bio diárias da Libra Esterlina, Dolar Canadense, Marco Alemão e Iene Japones, em

relação ao Dolar Americano, enquanto nos gráficos da direita estão os logaritmos das

taxas de variação (retornos diários). O peŕıodo amostral vai de janeiro de 1991 a de-

zembro de 1998. Uma caracteŕıstica comum nestes retornos é que embora as médias

pareçam ser aproximadamente constantes as variâncias mudam ao longo do tempo.

Na Figura 6.2 estão os histogramas com uma curva normal superimposta para os

mesmos dados (retornos). Pode-se notar que muitos valores aparecem nas caudas

das distribuições. Finalmente, na Figura 6.3 temos as autocorrelações amostrais dos

retornos e dos retornos ao quadrado. Note como existe bastante aucorrelação entre os

retornos ao quadrado. Todas estas caracteŕısticas são em geral verificadas em séries

reais de retornos e devem ser levadas em conta pelo modelo.

A idéia aqui é tentar modelar simultaneamente a média e a variância de uma série

temporal. Para fixar idéias, suponha que um modelo AR(1), Xt = αXt−1 + ǫt foi

73

74 CAPÍTULO 6. MODELANDO A VARIÂNCIA

0. 50

0. 60

0. 70

y. Li

br a

E st

er lin

a

1. 2

1. 4

y. D

ol ar

C an

ad en

se

1. 4

1. 6

1. 8

y. M

ar co

A le

m ao

80 10

0 12

0 14

0

0 500 1000 1500 2000

y. Ie

ne J

ap on

es

− 0.

03 0.

00 0.

02

z. Li

br a

E st

er lin

a

− 0.

01 5

0. 00

0 0.

01 5

z. D

ol ar

C an

ad en

se

− 0.

03 0.

00 0.

02

z. M

ar co

A le

m ao

− 0.

04 0.

00

0 500 1000 1500 2000

z. Ie

ne J

ap on

es

Figura 6.1: Taxas de câmbio e retornos diários em relação ao Dolar Americano da

Libra Esterlina, Dolar Canadense, Marco Alemão e Iene Japones, entre janeiro de

1991 a dezembro de 1998.

estimado e deseja-se fazer previsões 1 passo à frente,

x̂t(1) = E(Xt+1|xt) = αxt.

A variância condicional de Xt+1 é dada por

V ar(Xt+1|xt) = V ar(ǫt+1|xt) = E(ǫ2t+1|xt).

Até agora assumimos que E(ǫ2t+1|xt) = σ2ǫ , mas suponha que a variância condicional não seja constante, i.e. E(ǫ2t+1|xt) = σ2t+1. Uma posśıvel causa disto é que os dados se distribuem com caudas muito longas. Para facilitar a notação vamos denotar por

It = {xt, xt−1, . . . , ǫt, ǫt−1, . . . }, ou seja σ2t = V ar(ǫt|It−1).

6.2 Modelos ARCH

Existem várias formas de especificar como a variância condicional (volatilidade) varia

com o tempo. Uma estratégia utilizada para modelar σ2t , proposta em Engle (1982),

6.2. MODELOS ARCH 75

Libra Esterlina

−0.03 −0.01 0.01 0.03

0 20

40 60

80

Dolar Canadense

−0.015 0.000 0.010

0 50

10 0

15 0

Marco Alemao

−0.03 −0.01 0.01 0.03

0 20

40 60

Iene Japones

−0.06 −0.02 0.02

0 20

40 60

Figura 6.2: Histogramas dos retornos diários do Exemplo 6.1.

consiste em assumir que ela depende dos quadrados dos erros passados, ǫt−1, ǫt−2, . . .

através de uma autoregressão. No caso mais simples, faz-se

ǫt = vt

c+ αǫ2t−1 (6.1)

onde {vt} é uma série puramente aleatória com média zero e variância igual a 1 e vt e ǫt são independentes. Segue que a esperança e a variância condicionais são dadas

por,

E(ǫt|It−1) = E(vt) √

c+ αǫ2t−1 = 0

E(ǫ2t |It−1) = σ2t = c+ αǫ2t−1 (6.2)

Neste caso dizemos que a variância segue um processo autoregressivo condicionalmente

heterocedástico de ordem 1, ARCH(1). Note que é necessário impor as restrições c > 0

e α ≥ 0 para que σ2t seja sempre positiva. Quando α = 0 a variância condicional é constante e ǫt é um processo condicionalmente homocedástico.

Além disso queremos garantir a estacionariedade da autoregressão de modo que a

restrição imposta é 0 < α < 1. Note também que (6.2) não inclui um termo de erro

e portanto não é um processo estocástico.

76 CAPÍTULO 6. MODELANDO A VARIÂNCIA

0 5 15 25 −

0. 04

0. 04

Libra Esterlina

0 5 15 25

− 0.

05 0.

10

Libra Esterlina^2

0 5 15 25

− 0.

04 0.

02

Dolar Canadense

0 5 15 25

− 0.

05 0.

10

Dolar Canadense^2

0 5 15 25 −

0. 06

0. 00

Marco Alemao

0 5 15 25

− 0.

05 0.

10

Marco Alemao^2

0 5 15 25

0. 00

0. 10

Iene Japones

0 5 15 25

− 0.

05 0.

10 0.

25

Iene Japones^2

Figura 6.3: Correlogramas dos retornos e retornos ao quadrado no Exemplo 6.1

A esperança e variância incondicionais podem ser obtidas como,

E(ǫt) = E[E(ǫt|It−1)] = 0 V ar(ǫt) = E(ǫ

2 t ) = E[E(ǫ

2 t |It−1)] = c+ αE(ǫ2t−1).

Se o processo é estacionário então E(ǫ2t ) = E(ǫ 2 t−1) = V ar(ǫt) e portanto

V ar(ǫt) = c

1− α.

Além disso,

Cov(ǫt, ǫt+k) = E(ǫtǫt+k) = E[E(ǫtǫt+k)|ǫt+k−1, . . . , ǫt−1] = E[ǫtE(vt+k

c+ αǫt+k−1)] = 0, para k > 0.

Ou seja, ao postular o modelo ARCH(1) estamos assumindo que os {ǫt} são não correlacionados.

Exemplo 6.2 : Para ilustração a Figura 6.4 apresenta dois processos ARCH de ordem

1 simulados a partir de uma sequência {vt} de 200 números aleatórios i.i.d. gerados de uma distribuição N(0, 1). A sequência {ǫt} foi construida usando a equação (6.1)

6.2. MODELOS ARCH 77

com c = 1 e α = 0, 8. Note como a sequência {ǫt} continua tendo média zero mas parece ter tido um aumento de volatilidade em alguns peŕıodos. Em um modelo

AR(1), a forma como esta estrutura nos erros afeta a série original depende do valor

do parâmetro autoregressivo e duas posśıveis situações são mostradas nos gráficos

inferiores da figura. Na Figura 6.5 temos o histograma dos valores {ǫt} gerados, com uma curva normal superimposta, além do gráfico de probabilidades normais (QQ-

plot normal). Note como há um excesso de valores nas caudas ocorrendo com uma

frequência maior do que seria esperado na distribuição normal.

processo aleatório

0 50 100 150 200

− 5

0 5

10

ε(t) = v(t) 1 + 0.8ε(t − 1)2

0 50 100 150 200

− 5

0 5

10

x(t) = 0.5x(t − 1) + ε(t)

0 50 100 150 200

− 5

0 5

10

x(t) = 0.9x(t − 1) + ε(t)

0 50 100 150 200

− 5

0 5

10

Figura 6.4: Processos autoregressivos com erros ARCH(1) simulados.

Basicamente a equação (6.2) nos diz que erros grandes (ou pequenos) em valor

absoluto tendem a ser seguidos por erros grandes (ou pequenos) em valor absoluto.

Portanto o modelo é adequado para descrever séries aonde a volatilidade ocorre em

grupos. Além disso, na equação (6.2) somando ǫ2t e subtraindo σ 2 t de ambos os lados

obtemos que

ǫ2t = c+ αǫ 2 t−1 + νt

com νt = ǫ 2 t − σ2t = σ2t (v2t − 1). Ou seja, o modelo ARCH(1) pode ser reescrito como

um AR(1) estacionário para ǫ2t com erros não normais (v 2 t ∼ χ21 se vt ∼ N(0, 1)) que

têm média zero e variância não constante. Portanto, a função de autocorrelação do

processo {ǫ2t } é dada por ρ(k) = αk e o correlograma amostral deve apresentar um

78 CAPÍTULO 6. MODELANDO A VARIÂNCIA

de ns

id ad

es

−4 −2 0 2 4 6

0. 00

0. 15

0. 30

−3 −2 −1 0 1 2 3

− 4

0 4

Q−Q plot Normal

quantis teoricos

qu an

tis a

m os

tr ai

s

Figura 6.5: Caracteristicas de um processo ARCH(1) simulado.

decaimento exponencial para zero.

Se o processo ARCH(1) for estacionário não é dif́ıcil calcular o seu coeficiente de

curtose que é dado por

κ = E(ǫ4t )

[V ar(ǫt)]2 .

Denotando por E(v4t ) = λ o quarto momento do erro segue que o quarto momento

condicional é dado por

E(ǫ4t |It−1) = E(v4t σ4t |It−1) = λE(σ4t |It−1) = λ(c+ αǫ2t−1)2.

(se assumirmos que vt ∼ N(0, 1) então λ = 3). Portanto, o quarto momento incondi- cional fica,

E(ǫ4t ) = E[E(ǫ 4 t |It−1)] = λE(c2 + α2ǫ4t−1 + 2cαǫ2t−1).

Se o processo é estacionário de quarta ordem então podemos escrever E(ǫ4t ) =

E(ǫ4t−1) = µ4 e portanto,

µ4 = λ(c 2 + α2µ4 + 2cα

c

1− α) = λc 2

(

1 + α

1− α

)

+ λα2µ4

e finalmente,

µ4 = λc2(1 + α)

(1− α)(1− λα2) .

6.2. MODELOS ARCH 79

O coeficiente de curtose então fica,

κ = λc2(1 + α)

(1− α)(1− λα2) (1− α)2

c2 =

λ(1− α2) 1− λα2 , α

2 < 1/λ

que é sempre maior do que λ. Ou seja, qualquer que seja a distribuição de vt o coe-

ficiente de curtose será maior do que a curtose de vt (desde que α > 0 e λ > 1). Em

particular, processos ARCH(1) têm caudas mais pesadas do que a distribuição nor-

mal e são portanto adequados para modelar séries temporais com esta caracteŕıstica.

Séries de retornos, como as do Exemplo 6.1, frequentemente apresentam caudas mais

pesados do que a normal devido ao excesso de curtose.

Previsões da Volatilidade

Suponha que uma série temporal Xt segue um processo ARCH(1), i.e. Xt = vt √ ht,

vt ∼ N(0, 1). As previsões da volatilidade, k passos à frente, no tempo t = n são obtidas como,

ĥn(k) = E(hn+k|In) = c+ αE(X2n+k−1|In). Para k = 1 segue que E(X2n+k−1|In) = X2n+k−1 e para k > 1 temos que

E(X2n+k−1|In) = E(hn+k−1v2n+k−1|In) = E(hn+k−1|In)E(v2n+k−1|In) = E(hn+k−1|In) = ĥn(k − 1)

pois hn+k−1 e vn+k−1 são independentes. As previsões então ficam,

ĥn(k) =

{

c+ αX2n, k = 1

c+ αĥn(k − 1), k = 2, 3, . . .

O Modelo ARCH(p)

Estas idéias podem ser generalizadas para processos mais gerais ARCH(p) em que a

variância condicional depende dos quadrados de p erros passados, i.e.

ǫt = vt

c+ α1ǫ2t−1 + · · ·+ αpǫ2t−p (6.3)

e então a variância condicional é modelada como,

σ2t = E(ǫ 2 t |It−1) = c+ α1ǫ2t−1 + · · ·+ αpǫ2t−p.

Neste caso, para garantir que σ2t seja sempre positiva é necessário impor a seguintes

restrições c > 0 e α1 ≥ 0, . . . , αp ≥ 0 e para garantir estacionariedade é necessário também que as ráızes de 1− α1B − · · · − αpBp = 0 estejam fora do ćırculo unitário. Juntando estas restrições equivale a impor a restrição c > 0 e

∑p i=1 αi < 1.

Analogamente podemos reescrever o modelo ARCH(p) como um modelo AR(p)

para ǫ2t definindo os erros νt como anteriormente, i.e.

ǫ2t = c+ α1ǫ 2 t−1 + · · ·+ αpǫ2t−p + νt.

com νt = σ 2 t (v

2 t − 1).

80 CAPÍTULO 6. MODELANDO A VARIÂNCIA

Identificação

A caracteŕıstica chave dos modelos ARCH é que a variância condicional dos erros ǫt se

comporta como um processo autoregressivo. Portanto deve-se esperar que os reśıduos

de um modelo ARMA ajustado a uma série temporal observada também sigam este

padrão caracteŕıstico. Em particular, se o modelo ajustado for adequado então a FAC

e a FACP dos reśıduos devem indicar um processo puramente aleatório, no entanto

se a FAC dos quadrados dos reśıduos, ǫ̂2t , tiver um decaimento caracteŕıstico de uma

autoregressão isto é uma indicação de que um modelo ARCH pode ser apropriado. A

ordem p do modelo pode ser identificada através da FACP dos quadrados dos reśıduos.

Previsões da Volatilidade

Suponha que uma série temporal Xt segue um processo ARCH (p). As previsões da

volatilidade, k passos à frente, no tempo t = n são obtidas como,

ĥn(k) = E(hn+k|In) = c+ p

i=j

αjE(X 2 n+k−j |In).

Para k ≤ j segue que E(X2n+k−j |In) = X2n+k−j e para k > j temos que

E(X2n+k−j |In) = E(hn+k−jv2n+k−j |In) = E(hn+k−j |In)E(v2n+k−j |In) = E(hn+k−j |In) = ĥn(k − j)

já que hn+k−1 e vn+k−1 são independentes.

6.3 Modelos GARCH

Uma generalização natural dos modelos ARCH consiste em assumir que a variância

condicional se comporta como um processo ARMA, i.e. depende também de seus

valores passados. Fazendo ǫt = vt √ ht onde

ht = c+

p ∑

i=1

αiǫ 2 t−i +

q ∑

j=1

βjht−j

segue que a esperança condicional de ǫt é zero e a variância condicional é

σ2t = ht. Este modelo é chamado ARCH generalizado, ou GARCH, de ordem

(p, q). Aqui as restrições impostas sobre os parâmetros são dadas por c > 0 e ∑p

i=1 αi + ∑q

j=1 βj < 1.

Embora a primeira vista pareça um modelo mais complexo, sua vantagem sobre

os modelos ARCH é basicamente a parcimônia. Assim como um modelo ARMA pode

ser mais parcimonioso no sentido de apresentar menos parâmetros a serem estimados

do que modelos AR ou MA, um modelo GARCH pode ser usado para descrever a

volatilidade com menos parâmetros do que modelos ARCH.

6.3. MODELOS GARCH 81

Em termos de identificação dos valores de p e q as ferramentas básicas são mais

uma vez a FAC e a FACP dos quadrados dos reśıduos. Assim, se o modelo ajustado

for adequado a FAC e a FACP dos reśıduos devem indicar um processo puramente

aleatório, no entanto quando estas funções são aplicadas aos quadrados dos reśıduos

elas devem indicar um processo ARMA(p, q). A identificação pode não ser muito fácil

em algumas aplicações embora na maioria dos casos um modelo GARCH(1,1) seja

suficiente. Na prática recomenda-se também tentar outros modelos de ordem baixa

como GARCH(1,2) e GARCH(2,1).

As previsões da volatilidade em modelos GARCH são obtidas de forma similar

a de um modelo ARMA. Por exemplo, após estimar os parâmetros de um modelo

GARCH(1,1) e assumindo-se que ǫ0 = h0 = 0 pode-se construir as sequências ǫ1, . . . , ǫt e h1, . . . , ht e a previsão 1 passo à frente da volatilidade fica

σ̂2t (1) = c+ αǫ 2 t + βht.

6.3.1 Estimação

Para uma série x1, . . . , xn observada e um modelo GARCH(p, q), denotando-se o vetor

de parâmetros por θ=(c, α1, . . . , αp, β1, . . . , βq) e destacando-se a densidade conjunta

das p primeiras realizações segue que

p(x1, . . . , xn|θ) = p(x1, . . . , xp|θ) n ∏

t=p+1

p(xt|xt−1, . . . , xt−p,θ).

Assumindo normalidade segue que

Xt|xt−1, . . . , xt−p ∼ N(0, ht)

e portanto

p(x1, . . . , xn|θ) = p(x1, . . . , xp|θ) n ∏

t=p+1

(2πht) −1/2 exp(−(1/2)x2t /ht).

Em geral o número de observações será grande o suficiente para que o termo

p(x1, . . . , xp|θ) possa ser desprezado. Por exemplo, para um modelo ARCH(1) a função log-verossimilhança fica

−0.5 n ∑

t=2

[

log(2π) + log(c+ αx2t−1) + x 2 t /(c+ αx

2 t−1)

]

.

Note que algum algoritmo de otimização não linear deverá ser utilizado e nada garante

sua convergência para um ótimo global. No R pode-se usar a função garch do pacote

tseries para fazer a estimação por máxima verossimilhança.

82 CAPÍTULO 6. MODELANDO A VARIÂNCIA

6.3.2 Adequação

Se um modelo ARCH ou GARCH foi ajustado a uma série Xt não correlacionada

então os reśıduos padronizados são dados por

X̃t = Xt√ ht

e formam uma sequência i.i.d. com distribuição normal padrão. Assim, a adequa-

ção do modelo pode ser verificada aplicando os testes usuais de normalidade a estes

residuos padronizados e os testes de aleatoriedade (Box-Pierce e Ljung-Box) aos qua-

drados dos reśıduos.

Exemplo 6.3 : Na parte superior da Figura 6.6 estão os preços diários no fechamento

de um indice de mercado da Alemanha (DAX), entre 1991 e 1998. O interesse é em

analisar os chamados retornos dados por log(xt/xt−1) e estes estão no gráfico inferior

da Figura 6.6. Existe evidência na literatura que modelos GARCH(1,1) conseguem

captar bem os movimentos caracteŕısticos dos retornos.

D A

X

1992 1993 1994 1995 1996 1997 1998

20 00

50 00

re to

rn os

1992 1993 1994 1995 1996 1997 1998

− 0.

10 0.

00

Figura 6.6: Preços diários no fechamento de um indice de mercado da Alemanha (DAX), entre 1991 e 1998 e respectivos retornos.

6.4. VOLATILIDADE ESTOCÁSTICA 83

Usando a função garch no pacote tseries do R o modelo ajustado obtido foi

Yt = vt √

ht, vt ∼ N(0, 1) ht = 0.000005 + 0.068Y

2 t−1 + 0.889ht−1

sendo todos os coeficientes significativos. O teste de Ljung-Box aplicado nos quadra-

dos dos residuos indicou aleatoriedade (p-valor = 0,71), no entanto o teste de norma-

lidade de Jarque-Bera aplicado aos residuos rejeitou a hipótese nula (p-valor<0,001).

Assim a hipótese de normalidade condicional parece estar sendo violada.

Na Figura 6.7 estão os histogramas, gráficos de probabilidades normais dos retor-

nos e reśıduos do modelo GARCH(1,1) estimado, além dos correlogramas dos qua-

drados dos retornos e reśıduos. Use os comandos abaixo para estimar o modelo.

> library(tseries)

> data(EuStockMarkets)

> x = EuStockMarkets

> dax = diff(log(x))[, "DAX"]

> dax.garch = garch(dax, trace = F)

> r = dax.garch$residuals

> round(dax.garch$coef, 6)

a0 a1 b1

0.000005 0.068329 0.889067

Um fato estilizado presente em séries temporais financeiras é que o mercado tem

baixa volatilidade quando está em alta e alta volatilidade quando está em baixa.

Tal assimetria não é levada em conta pelos modelos GARCH e para contornar esta

limitação outros modelos foram propostos na literatura. Por exemplo, no modelo

EGARCH (ou GARCH exponencial) modela-se o logaritmo da volatilidade como,

log(σ2t ) = c+ α

ǫt−1 σt−1

+ γ ǫt−1 σt−1

+ βσ2t−1.

Em termos de estimação uma vantagem deste modelo é que os parâmetros c, α e β são

irrestritos já que estamos modelando o logaritmo da volatilidade. A única restrição é

γ < 0 pois assim a volatilidade aumenta quando ǫt−1 < 0.

6.4 Volatilidade Estocástica

As fórmulas para modelar σ2t vistas até agora foram todas determińısticas, i.e. sem

uma componente de erro aleatório. No entanto, pode ser mais razoável assumir que a

variância condicional varia estocasticamente ao longo do tempo ao invés de determi-

nisticamente, especialmente se existem mudanças abruptas na volatilidade (e.g. como

resultado de greves, guerras, etc.).

84 CAPÍTULO 6. MODELANDO A VARIÂNCIA

DAX

−0.10 −0.05 0.00 0.05 0

10 30

Residuos

−10 −5 0 5

0. 0

0. 2

−3 −2 −1 0 1 2 3

− 0.

10 0.

00 0.

10

DAX

−3 −2 −1 0 1 2 3

− 5

0 5

Residuos

0 5 10 15 20 25 30

0. 0

0. 4

0. 8

DAX^2

0 5 10 15 20 25 30

0. 0

0. 4

0. 8

Residuos^2

Figura 6.7: Histogramas e probabilidades normais dos retornos do indice de mercado da Alemanha (DAX) e reśıduos do modelos GARCH(1,1) e correlogramas dos seus quadrados.

Assim, uma alternativa aos modelos ARCH ou GARCH consiste em assumir que

σ2t segue um processo estocástico. Geralmente modela-se o logaritmo de σ 2 t . Em sua

forma mais simples um modelo de volatilidade estocástica (VE) é dado por

Xt = vt exp(ht/2), vt ∼ N(0, 1) ht = c+ αht−1 + ηt, ηt ∼ N(0, σ2η)

com |α| < 1 e ht = log(σ2t ). Note que não há necessidade de restrições de positividade nos parâmetros pois estamos modelando o logaritmo da volatilidade. O modelo pode

ser estendido para uma estrutura AR(p) em ht, ou seja

Xt = vt exp(ht/2), vt ∼ N(0, 1)

ht = c+

p ∑

i=1

αiht−i + ηt, ηt ∼ N(0, σ2η)

Propriedades

1. E(Xt) = E(vt e ht/2) = E(eht/2)E(vt) = 0, já que ht e vt são independentes.

6.4. VOLATILIDADE ESTOCÁSTICA 85

2. V ar(Xt) = E(X 2 t ) = E(e

ht v2t ) = E(e ht)E(v2t ) = E(e

ht). Mas, como estamos

assumindo que ht é estacionária segue que,

E(ht) = µ e V ar(ht) = σ 2 h = τ

2/(1− φ2)

e a distribuição incondicional do log-volatilidade é ht ∼ N(µ, σ2h). Portanto, eht é log-normal e segue que

V ar(Xt) = E(e ht) = eµ+σ

2

h /2.

3. E(X4t ) = E(v 4 t e

2ht) = 3 e2µ+2σ 2

h e a curtose é dada por

K = 3 e2µ+2σ

2

h

e2µ+σ 2

h

= 3eσ 2

h > 3.

Exerćıcios

1. Um modelo ARIMA foi identificado e estimado para uma série temporal obser-

vada mas há indicação de que a variância condicional deve ser modelada por

um processo GARCH(1,1). Explique como se chegou a esta conclusão.

2. Refaça o exemplo da Figura 6.4 e estime um modelo AR(1) para a sérieXt. Veri-

fique se existe estrutura autoregressiva nos quadrados dos reśıduos e identifique

um modelo ARCH para os erros.

3. Obtenha as previsões 1, 2 e 3 passos a frente para um modelo GARCH(1,2).

4. Descreva duas vantagens de modelos EGARCH sobre modelos GARCH.

Caṕıtulo 7

Modelos Lineares Dinâmicos

A classe de modelos lineares dinâmicos (MLD), também conhecidos como modelos

de espaço de estados tem sido utilizada com sucesso em análise e previsão de séries

temporais. Neste caṕıtulo serão apresentadas as formas mais comumente utilizadas

de MLD, maiores detalhes podem ser obtidos em West & Harrison (1997) e Pole,

West, & Harrison (1994).

7.1 Introdução

Um modelo linear dinâmico pode ser caracterizado pelo seguinte par de equações

yt = F ′ tθt + ǫt

θt = Gtθt−1 + ωt (7.1)

chamadas equações de observação e evolução respectivamente, onde θt denota o vetor

de estados no tempo t, F é um vetor de constantes conhecidadas ou regressores, G

é uma matrix de evolução conhecida. Os erros ǫt e ωt são geralmente assumidos não

correlacionados em todos os peŕıodos de tempo e serialmente não correlacionados com

média zero. Em muitas aplicações práticas pode-se assumir também que ǫt ∼ N(0, σ2ǫ ) e ωt tem distribuição normal multivariada com média zero e matriz de variância-

covariância W t.

A idéia aqui é que a “idade” da informação que se tem sobre θ seja levada em

conta no sentido de que nossa incerteza a respeito de θ deve aumentar com o passar

do tempo. Neste sentido, a forma do modelo é apropriada apenas “localmente” no

tempo e é necessário caracterizar algum tipo de evolução temporal de θ. O que se

tem então é uma sequência de modelos ou um modelo dinâmico parametrizado por θt (o estado do processo no tempo t).

Considere um modelo em que uma variável y está relacionada a uma outra variável

X de acordo com a seguinte forma paramétrica

y = Xθ + ǫ.

Além disso, a incerteza do pesquisador em relação ao parâmetro θ é descrita em

termos de uma distribuição de probabilidades p(θ).

86

7.1. INTRODUÇÃO 87

Em um peŕıodo t qualquer, Dt representa o conjunto de informações dispońıveis

sobre θ. Por simplicidade vamos assumir que Dt = {y1, . . . , yt}. Neste sentido, D0 representa toda a informação inicial (antes de observar os dados) relevante sobre θ

incluindo a própria definição do modelo.

No tempo t − 1, após observar y1, . . . , yt−1, toda a informação sobre o estado do processo está resumida probabilisticamente na distribuição a posteriori p(θt−1|Dt−1). No tempo t, antes de observar yt, toda a informação histórica Dt−1 está resumida

probabilisticamente na distribuição a priori de θt obtida como

p(θt|Dt−1) = ∫

p(θt|θt−1)p(θt−1|Dt−1)dθt−1

que é atualizada após observar yt para a posteriori θt, combinando-se com o modelo

amostral p(yt|θt) via teorema de Bayes

p(θt|Dt) = p(yt|θt)p(θt|Dt−1)

p(yt|Dt−1)

sendo

p(yt|Dt−1) = ∫

p(yt|θt)p(θt|Dt−1)dθt

é a distribuição preditiva de yt. Esquematicamente,

θt−1|Dt−1 −→ θt|Dt−1 −→ θt|Dt Posteriori Priori Posteriori

↓ Yt|Dt−1 Previsão

Estas equações fornecem um sistema de aprendizado sequencial sobre os parâme-

tros do processo (não observáveis) e também uma sequência de distribuições preditivas

(1 passo a frente) para as quantidades observáveis. Porém a sua implementação prá-

tica envolve a resolução de integrais que pode ser um problema de dif́ıcil solução

em casos mais gerais. Um caso particular, onde as equações podem ser escritas em

forma fechada, é o de modelos lineares dinâmicos (MLD) normais onde a distribuição

amostral é definida pela

equação das observações yt = Xtθt + ǫt, ǫt ∼ N(0, Vt)

e os parâmetros se relacionam em peŕıodos sucessivos através da

equação do sistema θt = Gθt−1 + ωt, ωt ∼ N(0,Wt)

onde as sequências ǫt e ωt são independentes, mutuamente independentes e ambas

são independentes da informação inicial θ0|D0 ∼ N(m0, C0). A matriz G descreve a evolução (determińıstica) dos parâmetros. Modelos nesta classe serão analisados nas

próximas seções.

88 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

7.2 Modelos Polinomiais

No MLD mais simples as observações são representadas por

yt = µt + ǫt, ǫt ∼ N(0, Vt)

onde µt é o ńıvel da série no tempo t. A evolução do ńıvel é modelada como um

passeio aleatório simples, i.e.

µt = µt−1 + ωt, ωt ∼ N(0,Wt).

Estas equações podem ser reescritas como

yt|µt ∼ N(µt, Vt) µt|µt−1 ∼ N(µt−1,Wt)

e a informação inicial é µ0|D0 ∼ N(m0, C0). Vamos assumir por enquanto que as vari- âncias Vt e Wt são conhecidas. Este modelo pode ser pensado como uma aproximação

de Taylor de 1a ordem para uma função suave do tempo µ(t) de modo que

µ(t+ δt) = µ(t) + termos de ordem mais alta

e o modelo descreve os termos de ordem mais alta simplesmente como rúıdos de média

zero. Como saber então se este modelo é adequado a uma particular aplicação?

No tempo t, o valor esperado da série k peŕıodos a frente condicional ao ńıvel

atual é

E(Yt+k|µt) = E(µt+k|µt) = E(µt + k

i=1

ωt+i|µt) = µt

e denotando a média da distribuição a posteriori de µt por mt então a função de

previsão é constante

ft(k) = E(Yt+k|Dt) = E[E(Yt+k|µt, Dt)] = E(µt|Dt) = mt, ∀k > 0.

Assim, este modelo é útil para previsões de curto prazo, particularmente quando a

variação das observações (medida por Vt) é muito maior do que a variação do ńıvel

(medida por Wt).

Exemplo 7.1 : Foram gerados 100 valores de um modelo polinomial de primeira

ordem com variâncias constantes (Vt = V e Wt = W ). Na Figura 7.1 estão os valores

gerados com as relações V/W iguais a 20, 2 e 0,2. Seguem os comandos do R para

produção dos gráficos.

> mld.sim = function(n, V, W, mu0) {

+ mu = mu0 + cumsum(rnorm(n, sd = sqrt(W)))

+ obs = mu + rnorm(n, sd = sqrt(V))

+ ts(cbind(obs, mu))

+ }

7.2. MODELOS POLINOMIAIS 89

7.2.1 Análise Sequencial e Previsões

A média inicial m0 é uma estimativa pontual do ńıvel da série e a variância inicial

C0 mede a incerteza associada. Assumindo que µt−1|Dt−1 ∼ N(mt−1, Ct−1), então condicionalmente a Dt−1, µt é a soma de 2 quantidades normais e independentes µt−1 e ωt e portanto é também normal com média e variância dadas por

E(µt|Dt−1) = E(µt−1|Dt−1) + E(ωt|Dt−1) = mt−1 V ar(µt|Dt−1) = V ar(µt−1|Dt−1) + V ar(ωt|Dt−1) = Ct−1 +Wt = Rt

Yt|Dt−1 é também a soma de quantidades normais independentes e portanto tem distribuição normal com

E(Yt|Dt−1) = E(µt|Dt−1) + E(ǫt|Dt−1) = mt−1 V ar(Yt|Dt−1) = V ar(µt|Dt−1) + V ar(ǫt|Dt−1) = Rt + Vt = Qt

Após observar yt, a distribuição atualizada do ńıvel é obtida via teorema de Bayes

combinando-se a verossimilhança

p(yt|µt, Dt−1) = (2πVt)−1/2 exp{−(yt − µt)2/2Vt}

com a priori

p(µt|Dt−1) = (2πRt)−1/2 exp{−(µt −mt−1)2/2Rt}

de modo que

p(µt|Dt) ∝ exp {

−1 2

[

(yt − µt)2 Vt

+ (µt −mt−1)2

Rt

]}

∝ exp {

−1 2

[

µ2t (V −1 t +R

−1 t )− 2µt(V −1t yt +R−1t mt−1)

]

}

∝ exp {

−C −1 t

2 (µ2t − 2µtmt)

}

∝ exp {

−C −1 t

2 (µt −mt)2

}

onde

mt = Ct(V −1 t yt +R

−1 t mt−1)

C−1t = V −1 t +R

−1 t

e todos os termos que não dependem de µt foram colocados na constante de propor-

cionalidade. Portanto, µt|Dt ∼ N(mt, Ct). A média a posteriori pode ser reescrita de 2 formas alternativas definindo-se o

coeficiente adaptativo At = CtV −1 t = Rt/Qt ∈ (0, 1) e o erro de previsão 1 passo a

frente et = yt −mt−1. Assim

mt = (1−At)mt−1 +Atyt = mt−1 +Atet.

Note a similaridade com a equação de previsão do método de alisamento exponencial

simples visto no Caṕıtulo 5. Aqui At faz o papel da constante de alisamento porém

90 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

agora variando no tempo. A variância a posteriori também pode ser reescrita como

função do coeficiente adaptativo como

Ct = Rt −A2tQt < Rt.

Podemos utilizar as equações das observações e de evolução para obter a distri-

buição preditiva k passos a frente. Fazendo substituições sucessivas obtemos que

µt+k = µt + k

j=1

ωt+j

Yt+k = µt +

k ∑

j=1

ωt+j + ǫt+k

e como todos os termos são normais e independentes segue que Yt+k é também normal

com

E(Yt+k|Dt) = E(µt|Dt) = mt

V ar(Yt+k|Dt) = Ct + k

j=1

Wt+j + Vt+k

A função abaixo estima um modelo com tendencia polinomial de 1a ordem fa-

zendo a analise sequencial usando as equações dadas no texto com variâncias fixas e

conhecidas.

> mld = function(Y, V, W, m0, C0) {

+ n = length(Y)

+ m = C = R = Q = f = A = e = ts(rep(NA, length = n), start = start(Y))

+ Y = ts(c(NA, Y), end = end(Y))

+ C[1] = C0

+ m[1] = m0

+ for (t in 2:n) {

+ R[t] = C[t - 1] + W[t]

+ f[t] = m[t - 1]

+ Q[t] = R[t] + V[t]

+ A[t] = R[t]/Q[t]

+ e[t] = Y[t] - f[t]

+ m[t] = m[t - 1] + A[t] * e[t]

+ C[t] = A[t] * V[t]

+ }

+ return(list(m = m, C = C, R = R, f = f, Q = Q))

+ }

Exemplo 7.2 : A função mld pode ser usada para estimar sequencialmente o nivel

da serie de vazões do rio Nilo. Primeiro vamos permitir que o nivel varie bastante ao

7.2. MODELOS POLINOMIAIS 91

longo do tempo especificando um valor grande para W e depois reestimar com pouca

variação temporal (W bem pequeno). Usaremos a variancia amostral da serie como

estimativa de V . Como informação inicial usaremos uma estimativa do nivel igual

a 1000 mas com uma grande incerteza associada. O gráfico da série com os ńıveis

superimpostos aparece na Figura 7.2.

7.2.2 Variâncias de Evolução e das Observações

Tipicamente, Wt é desconhecida. Sua estimação entretanto leva a uma intratabilidade

anaĺıtica que pode ser evitada através de sua especificação subjetiva.

O fator de desconto é o parâmetro básico que controla o grau de “envelhecimento”

da informação de uma observação. Por exemplo, podemos quantificar o envelheci-

mento da informação sobre o parâmetro µt como um aumento de 5% em sua variância

a priori (no tempo t), i.e.

V ar(µt|Dt−1) = (1 + δ)V ar(µt−1|Dt−1) ou Rt = (1 + δ) Ct−1

com δ = 0.05. Por outro lado, informação é em geral medida em termos de precisão

(o inverso da variância) e podemos escrever

Precisão(µt|Dt−1) = (1 + δ)−1Precisão(µt−1|Dt−1) ou R−1t = (1 + δ)−1C−1t−1.

Nesta escala, o fator de desconto λ = (1 + δ)−1 varia entre 0 e 1 e δ = 5% implica

em λ ≈ 0.95. Vale notar que o fator de desconto não depende da escala na qual as observações são medidas.

Se λ = 1 então não existe mudança ao longo do tempo no ńıvel da série e quanto

menor é o valor de λ maiores são as alterações esperadas e maior é a perda de infor-

mação contida em observações mais “antigas”.

Assim, para um valor fixo do fator de desconto λ temos que

Rt = Ct−1/λ = Ct−1 +Wt

ou equivalentemente

Wt = Ct−1

(

1− λ λ

)

= δCt−1.

Como Rt = Ct−1 + Wt podemos interpretar esta especificação intuitivamente como

um aumento de incerteza, ao evoluir de t−1 para t, quantificado como uma proporção de Ct−1.

A sequência de variâncias Vt também é, em geral, desconhecida embora o pesqui-

sador possa ter alguma informação a priori sobre caracteŕısticas desta sequência. Por

exemplo, Vt = V (variância constante e desconhecida), Vt = V kt onde os pesos kt são

conhecidos, Vt = V k(µt) onde k(·) é uma função de variância do ńıvel da série ou em particular Vt = V µ

p t .

Impondo-se uma particular estrutura para a sequência Wt e para a informação

inicial obtem-se um procedimento de atualização sequencial para V além de µt. Para

92 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

isto redefine-se o modelo, agora condicionalmente em V , como

yt = µt + ǫt, ǫt ∼ N(0, V ), µt = µt−1 + ωt, ωt ∼ N(0, V W ∗t ),

µ0|V,D0 ∼ N(m0, V C∗0 )

V −1|D0 ∼ Gama (

n0 2 , n0S0 2

)

ou n0S0V −1 ∼ χ2n0

sendo que m0, C ∗ 0 , n0 e S0 serão especificados. Surgiu assim mais um item na infor-

mação inicial com

E(V −1|D0) = n0/2

n0S0/2 =

1

S0

e S0 é a estimativa pontual a priori da variância V . Com esta definição pode-se

mostrar que a distribuição inicial marginal de µ0 é

µ0|D0 ∼ tn0(m0, C0)

com C0 = S0C ∗ 0 .

Se a distribuição a posteriori (incondicional) do ńıvel em t− 1 é

µt−1|Dt−1 ∼ tnt−1(mt−1, Ct−1)

então pode-se mostrar que as distribuições a priori, preditiva e a posteriori no tempo

t são dadas por

µt|Dt−1 ∼ tnt−1(mt−1, Rt) Yt|Dt−1 ∼ tnt−1(mt−1, Qt) µt|Dt ∼ tnt(mt, Ct)

onde os parâmetros atualizados são

Qt = Rt + St−1

mt = mt−1 +Atet

Ct = (St/St−1)(Rt −A2tQt) nt = nt−1 + 1

ntSt = nt−1St−1 + St−1e 2 t /Qt.

A função mld1 abaixo faz a análise sequencial com a variância das observações

fixa e desconhecida. A especificação de Wt é feita via fator de desconto. Note que

agora tanto o nivel quanto a variância e os graus de liberdade são atualizados sequen-

cialmente.

> mld1 = function(Y, delta, m0, C0, n0, S0) {

+ N = length(Y)

+ m = n = C = R = Q = S = f = A = e = rep(NA, length = N)

+ Y = c(NA, Y)

7.2. MODELOS POLINOMIAIS 93

+ C[1] = C0

+ m[1] = m0

+ S[1] = S0

+ n[1] = n0

+ for (i in 2:N) {

+ n[i] = n[i - 1] + 1

+ R[i] = C[i - 1]/delta

+ f[i] = m[i - 1]

+ Q[i] = R[i] + S[i - 1]

+ A[i] = R[i]/Q[i]

+ e[i] = Y[i] - f[i]

+ S[i] = S[i - 1] + (S[i - 1]/n[i]) * (e[i]^2/Q[i] - 1)

+ m[i] = m[i - 1] + A[i] * e[i]

+ C[i] = A[i] * S[i]

+ }

+ return(list(m = m, C = C, R = R, f = f, Q = Q, n = n, S = S,

+ e = e))

+ }

Exemplo 7.3 : Novamente vamos examinar a série de vazões do rio Nilo, agora

usando diferentes fatores de desconto na função mld1.

> res1 = mld1(y, delta = 0.98, m0 = 1000, C0 = 100, n0 = 1, S0 = 0.01)

> res2 = mld1(y, delta = 0.7, m0 = 1000, C0 = 100, n0 = 1, S0 = 0.01)

Os gráficos na Figura 7.3 mostram a série original, as estimativas do nivel obtidas

com descontos 0,70 e 0,98 e estas mesmas estimativas com um intervalo de ±1, 5 √ Ct.

Os gráficos foram feitos com os seguintes comandos do R,

Os modelos podem ser comparados calculando-se o erro quadrático médio e o

desvio absoluto médio. Usando os comandos abaixo percebe-se que o modelo com

fator de desconto 0,70 é melhor segundo estes critérios.

> eqm = dam = rep(0, 2)

> for (i in 2:length(y)) {

+ eqm[1] = eqm[1] + (y[i] - res1$m[i - 1])^2

+ dam[1] = dam[1] + abs(y[i] - res1$m[i - 1])

+ eqm[2] = eqm[2] + (y[i] - res2$m[i - 1])^2

+ dam[2] = dam[2] + abs(y[i] - res2$m[i - 1])

+ }

> eqm

[1] 2681716 2375484

> dam

[1] 13258.47 11904.16

94 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

7.3 Modelo de Crescimento Linear

Considere agora que a descrição “local” mais apropriada é uma tendência polinomial

de 2a ordem. Um modelo um pouco mais elaborado é então obtido criando-se um

parâmetro extra para descrever o crescimento do ńıvel do processo observado. A

equação das observações fica inalterada, i.e.

yt = µt + ǫt, ǫt ∼ N(0, Vt)

e a evolução do ńıvel e do crescimento é modelada como

µt = µt−1 + βt−1 + ω1t

βt = βt−1 + ω2t.

Usando a representação matricial temos que o vetor de regressão e a matriz de

evolução são dados por

Xt = ( 1 0 ) e Gt =

(

1 1

0 1

)

.

Nesta notação, definindo θt = (µt, βt) ′ obtemos os momentos das distribuições a priori

e preditiva como

E(θt|Dt−1) = at = GE(θt−1|Dt−1) = Gmt−1 = (mt−1 + bt−1, bt−1)′

V ar(θt|Dt−1) = Rt = GCt−1G′ +Wt E(Yt|Dt−1) = ft = Xtat = mt−1 + bt−1

V ar(Yt|Dt−1) = Qt = XtRtX ′t + St−1.

Os momentos da distribuição a posteriori de θt são uma generalização matricial

daqueles obtidos para o modelo anterior,

E(θt|Dt) = mt = at +Atet V ar(θt|Dt) = Ct = (St/St−1)(Rt −AtA′tQt)

Não é dif́ıcil verificar que a função de previsão é dada por

ft(k) = XtG kmt = mt + kbt

sendo quemt e bt são as estimativas pontuais do ńıvel µt e do crescimento βt. Portanto,

assim como no caso anterior, este modelo também é apropriado para previsões de curto

prazo.

As variâncias Wt são mais uma vez especificadas indiretamente através de um

fator de desconto λ. Neste caso,

Rt = GCt−1G ′/λ

implica que

Wt = GCt−1G ′(λ−1 − 1).

7.4. MODELOS SAZONAIS 95

7.4 Modelos Sazonais

Um comportamento periódico ou ćıclico pode ser encontrado em várias séries tem-

porais. É importante que se consiga descrever o padrão sazonal da série através de

quantidades que possam ser estimadas incluindo-se assim este padrão na função pre-

visão. Nos modelos aqui analisados define-se um componente sazonal descrevendo

desvios sazonais em torno de um ńıvel dessazonalizado ou tendência.

7.4.1 Modelos sem Crescimento

A idéia aqui é fazer a superposição de um modelo polinomial de 1a ordem (para o ńıvel

dessazonalizado) com um modelo de efeitos sazonais. As equações das observações e

de evolução são dadas por

yt = µt + φt0 + ǫt, ǫt ∼ N(0, Vt) µt = µt−1 + ωt

φtr = φt−1,r+1 + ωt,r, r = 0, · · · , p− 2 φt,p−1 = φt−1,0 + ωt,p−1

com a restrição ∑p−1

r=0 φtr = 0, ∀t e onde p é o peŕıodo sazonal da série. Por exemplo, p = 12 para uma série com observações mensais e p = 4 para observações trimestrais.

Para fixar idéias, considere uma série trimestral e suponha que t− 1 é o segundo trimestre de um determinado ano. Então o vetor de parâmetros consiste de 4 efeitos

sazonais, um para cada trimestre,

φt−1 =

φt0 φt1 φt2 φt3

=

trim. 2

trim. 3

trim. 4

trim. 1

e ao passar de t − 1 para t ocorre simplesmente uma rotação nos elementos deste vetor,

φt =

φt0 φt1 φt2 φt3

=

trim. 3

trim. 4

trim. 1

trim. 2

.

A função de previsão assume a forma ft(k) = mt+htj onde mt é o valor esperado

do ńıvel dessazonalizado no tempo t + k e htj é o desvio sazonal esperado em torno

deste ńıvel. O desvio utilizado na função de previsão é tal que j é o resto da divisão

k/p. Por exemplo, se p = 12, e φt0 refere-se ao mês de janeiro então a previsão

1 passo a frente (k = 1) feita em dezembro é mt + E(φt0|Dt), com j = 1. Se o horizonte de previsão for k = 2 então j = 2 e o desvio sazonal refere-se a fevereiro,

i.e. ft(2) = mt + E(φt1|Dt).

96 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

7.4.2 Modelos com Crescimento

Novamente a idéia é fazer a superposição de um modelo para os efeitos sazonais

mas agora com um modelo polinomial de 2a ordem onde se tem um parâmetro que

representa o crescimento do ńıvel dessazonalizado.

O modelo pode ser escrito como

yt = µt + φt0 + ǫt, ǫt ∼ N(0, Vt) µt = µt−1 + βt−1 + ωt

βt = βt−1 + ω ∗ t

φtr = φt−1,r+1 + ωt,r, r = 0, · · · , p− 2 φt,p−1 = φt−1,0 + ωt,p−1

com a restrição ∑p−1

r=0 φtr = 0, ∀t. A função de previsão agora assume a seguinte forma

ft(k) = mt + kbt + htj , com

p−1 ∑

j=0

htj = 0

onde htj tem a mesma interpretação anterior.

7.5 Representação de Fourier

Uma forma alternativa de se representar padrões ćıclicos é através de combinações

lineares de funções periódicas. Em particular a utilização de funções trigonométricas

leva a representações de Fourier da sazonalidade.

O modelo (com crescimento) é representado pelas seguintes equações

yt = µt +

p/2 ∑

j=1

αj,t + ǫt, ǫt ∼ N(0, Vt)

µt = µt−1 + βt−1 + ωt,

βt = βt−1 + ω ∗ t ,

αj,t

α∗j,t

 =

[

cos 2πj/p sin 2πj/p

− sin 2πj/p cos 2πj/p

][

αj,t−1 α∗j,t−1

]

+

[

wj,t w∗j,t

]

, j = 1, . . . , p/2− 1

e αj,t = −αj,t−1 + ωj,t para j = p/2. A função de previsão é dada por

ft(k) =

p/2 ∑

j=1

Sjk =

p/2 ∑

j=1

[at,j cos(2πjk/p) + a ∗ t,jsen(2πjk/p)

onde at,j e a ∗ t,j são as estimativas pontuais de coeficientes de Fourier αt,j e α

∗ t,j .

Como no caṕıtulo anterior, as variâncias dos erros de evolução são especificadas

indiretamente através de um fator de desconto. A estratégia recomendada em (Pole,

West, & Harrison 1994) e West & Harrison (1997) consiste em especificar um fator de

7.6. ILUSTRAÇÃO 97

desconto para cada componente do modelo. No modelo com uma tendência polinomial

mais um componente sazonal teremos então 2 fatores de desconto.

Em geral, o fator de desconto do componente sazonal é maior do que o da tendên-

cia. Neste sentido estamos assumindo que o padrão sazonal da série, embora possa

estar sujeito a alterações, é mais estável do que a sua tendência.

7.6 Ilustração

A Figura ?? apresenta o total de vendas trimestrais (em milhares) de perus na Irlanda

entre o primeiro trimestre de 1974 e o terceiro trimestre de 1982. A série exibe

um crescimento sistemático ao longo de todo o peŕıodo juntamente com um padrão

sazonal acentuado. Outra caracteŕıstica interessante é que a forma do padrão sazonal

se alterou a partir de 1978. Vamos fazer a estimação sequencial de um modelo para

os efeitos sazonais superpostos a uma tendência de crescimento linear e verificar o

comportamento das previsões 1 passo a frente.

Suponha que a informação a priori foi acessada examinando-se as vendas dos anos

anteriores a 1974. Esta informação está resumida na Tabela 7.1. Note a restrição de

soma zero na especificação a priori dos efeitos sazonais e também que a especificação

equivalente em termos de fatores sazonais seria 11, 19, 19 e 11 para os fatores e

(11+19+19+11)/4 = 15 para o ńıvel.

Tabela 7.1: Informação a priori.

Componente Média (Desvio padrão)

Nı́vel 15 (0.75)

Crescimento 0 (0.3)

Efeito sazonal 1 -4 (0.5)

Efeito sazonal 2 4 (0.5)

Efeito sazonal 3 4 (0.5)

Efeito sazonal 4 -4 (0.5)

D.P. das observações 1 com 1 g.l.

A performance preditiva do modelo foi investigada para fatores de desconto va-

riando nos intervalos (0.9,1.0) para a tendência e (0.6,1.0) para os fatores sazonais.

Estes intervalos estão coerentes com a idéia de que espera-se um padrão sazonal mais

estável do que a tendência. Entretanto os valores encontrados após esta “busca” fo-

ram 0.90 para a tendência e 0.80 para os fatores sazonais. Uma idéia intuitiva é a

alteração no padrão sazonal ocorrida em 1978 deve ter contribuido para este resultado

at́ıpico.

Os 2 gráficos a seguir apresentam as previsões pontuais (1 passo a frente) junta-

mente com intervalos de 90% de probabilidade e os valores observados da série. O

primeiro gráfico refere-se ao modelo estático (ambos os fatores de desconto iguais a 1).

98 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

Note que a mudança no padrão sazonal ocorre muito lentamente no modelo estático e

no final da série o padrão estimado é apenas ligeiramente diferente do padrão inicial.

Já no modelo dinâmico o padrão sazonal evolui para uma forma completamente di-

Tabela 7.2:

Descontos EQM DAM LLIK

0.90 e 0.80 3.11 1.34 -71.1

1.00 e 1.00 4.23 1.64 -77.6

ferente melhorando a performance preditiva. Este fato pode ser notado por inspeção

visual e é confirmado pelos indicadores na Tabela 7.2.

A explicação intuitiva para este fato, lembrando da definição de fator de desconto,

é que no modelo dinâmico um peso maior é dado para as observações mais recentes

ao fazer previsões. Com isto a alteração no padrão sazonal é “incorporada” mais rapi-

damente do que no modelo estático. As previsões de vendas para o quarto trimestre

de 1982 e para 1983 também levarão em conta os diferentes padrões sazonais do final

da série.

7.7 Modelos de Regressão

Para completar o nosso modelo dinâmico podemos pensar em incluir na equação

das observações efeitos de variáveis regressoras. Considere por exemplo a regressão

linear da varável yt em uma coleção de p variáveis independentes X1t, . . . , Xpt. Se

um termo constante for incluido no modelo então X1t = 1, ∀t. Denotando o vetor de regressão e o vetor de coeficientes de regressão no tempo t por Xt = (X1t, . . . , Xpt) e

θt = (θ1t, . . . , θpt) ′ respectivamente então as equações do modelo são dadas por

yt = Xtθt + ǫt, ǫt ∼ N(0, Vt) θt = θt−1 + ωt, ωt ∼ N(0,Wt).

Assim, os coeficientes da regressão evoluem segundo um passeio aleatório, como

no modelo polinomial de 1a ordem, i.e., a matriz de evolução G = Ip. O vetor de

regressão é formado pelas próprias variáveis regressoras e note que a equação das

observações pode ser reescrita como

yt =

p ∑

i=1

θitXit + ǫt

de modo que o modelo pode ser visto como uma superposição de p regressões simples

pela origem.

Todas as distribuições envolvidas são análogas aos casos anteriores e as equações

dadas na Seção 2.3 podem ser utilizadas para obter os momentos das distribuições a

7.8. MONITORAMENTO 99

priori, preditiva e a posteriori fazendo-se G = Ip. Assim,

at = mt−1

Rt = Ct−1 +Wt

ft = Xtmt−1

e as outras equações permanecem inalteradas.

É interessante notar como fica a função de previsão ft(k) neste caso. Primeiro

reescreva a equação de evolução para θt+k fazendo k substituições sucessivas obtendo

θt+k = θt + k

j=1

ωt+j

de modo que

at+k = mt

Rt+k = Ct +

k ∑

j=1

Wt+j .

Então, usando a equação das observações obtemos que

ft(k) = Xt+kmt

Qt+k = Xt+kRt+kX ′ t+k + St.

Assim, a previsão pontual k passos a frente é a própria função de regressão avaliada

na estimativa dos coeficientes no tempo t e nos valores futuros dos regressores (que

nem sempre estão dispońıveis).

A sequência de variâncias Wt é mais uma vez estruturada usando um fator de

desconto.

7.8 Monitoramento

Ao comparar sequencialmente as previsões com os valores observados pode-se julgar

a adequação relativa de modelos alternativos com base em sua performance preditiva.

Observações ocorrendo nas caudas da distribuição preditiva são sempre posśıveis

por definição porém improváveis. Quanto mais afastada em uma das caudas mais

improvável é a observação. É preciso então estabelecer um critério para julgar que

tipo de inconsistência entre observação e previsão deve ser sinalizada pelo sistema.

No entanto, sinalizar uma observação como improvável apenas indica uma posśıvel

deficiência geral do modelo. É preciso saber em que sentido o modelo é deficiente,

i.e. verificar que modelos alternativos, com diferentes distribuições preditivas, teriam

uma performance melhor. O fator de Bayes, definido a seguir, é a ferramenta utilizada

para fazer esta comparação de modelos.

100 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

Se pA(yt|Dt−1) é a densidade preditiva 1 passo a frente de um modelo alternativo então o fator de Bayes é definido como

Ht = p(yt|Dt−1) pA(yt|Dt−1)

,

i.e. a razão das densidades preditivas avaliadas no valor observado yt.

Outra forma de comparar a performance preditiva de dois modelos é considerer um

grupo de observações ao invés de uma única e se basear no fator de Bayes acumulado

Ht(k) = p(yt, . . . , yt−k+1|Dt−k) pA(yt, . . . , yt−k+1|Dt−k)

= p(yt|Dt−1)p(yt−1, . . . , yt−k+1|Dt−k)

pA(yt|Dt−1)pA(yt−1, . . . , yt−k+1|Dt−k)

= HtHt−1(k − 1) = k−1 ∏

j=0

Ht−j .

Pode-se assim sinalizar evidências de alteração lenta na estrutura da série. A idéia

é que, individualmente, estas evidências não são suficientes para se “questionar” as

previsões do modelo em uso mas quando consideradas conjuntamente a evidência

acumulada pode ser grande e deve ser sinalizada. A questão agora é como construir

um sistema de monitoramento automático da série a partir destas idéias intuitivas.

Quando as observações estão cada vez mais afastadas das previsões então um fator

de Bayes individual Ht pode não ser suficientemente pequeno e precisa ser acumu-

lado para indicar alguma evidência contra o modelo padrão. Neste caso, o monitor

identifica o grupo mais discrepante de observações consecutivas calculando Vt and lt da seguinte forma,

Vt = min 1≤k≤t

Ht(k) = Ht(lt)

sendo calculado sequencialmente com as seguintes recursões,

Vt = Htmin{1, Lt−1} e lt = {

lt−1 + 1, se Lt−1 < 1

1, se Lt−1 ≥ 1

conforme mostrado em West (1986).

O modelo padrão é aceito como sendo satisfatório até a ocorrência de um valor

Lt menor do que um valor pré-especificado τ < 1 (o limite inferior para aceitação

de Lt) quando a ocorrência de uma descontinuidade na série é sinalizada. Se lt = 1

então uma única observação discrepante é identificada como a causa mais provável de

falha, embora o ińıcio de uma mudança também seja uma possibilidade. Por outro

lado, lt > 1 indica que uma mudança começou a ocorrer lt periods atrás em t− lt+1. Além disso, se uma mudança estrutural lenta está ocorrendo na série as observações

mais recentes indicarão evidência contra o modelo padrão que não será suficiente para

fazer Lt < τ . Assim, para aumentar a sensibilidade do monitor a estas mudanças uma

descontinuidade deve ser sinalizada se lt > 3 ou 4.

Para especificar o modelo alternativo assume-se que as densidades preditivas são

normais com média comum ft e variâncias Qt e Qt/ρ onde 0 < ρ < 1, de modo que o

fator de Bayes fica

Ht =

1

ρ exp

{

−(yt − ft) 2

2Qt (1− ρ)

}

=

1

ρ exp

{

−1 2 (1− ρ)e2t

}

7.8. MONITORAMENTO 101

onde et é o erro de previsão um passo a frente padronizado.

A escolha de ρ pode ser facilitada reescrevendo-se o fator de Bayes como

Ht = exp(−0.5 log ρ+ (1− ρ)e2t ).

Claramente Ht = 1 ou equivalentemente e 2 t = −(log ρ)/(1 − ρ) indica nenhuma evi-

dência para discriminar entre os modelos. O valor de ρ, pode ser escolhido de modo

a fornecer o valor máximo de |et| que não indica evidence contra o modelo padrão. Por exemplo, ρ ∈ (0.1, 0.3) implica que a evidência contra o modelo padrão deve ser acumulada para 1.3 < |et| < 1.6 que são aproximadamente os percentil 0.90 e 0.95 distribuição normal padrão.

É claro que para ρ fixo, a evidência contra o modelo padrão aumenta com |et|. West & Harrison (1997) ilustraram como a escolha de ρ tem pouca influência quando

o erro se torna muito grande em relação ao modelo alternativo. Este pode ser visto

como um modelo geral no sentido de levar em conta vários tipos de mudanças além

de observações discrepantes. Essencialmente, este procedimento pode ser visto como

um método exploratório gerando informação sobre o tipo e o peŕıodo mais provável

de mudança estrutural.

102 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

> w = c(0.05, 0.5, 5)

> g = list(col = 1:2, xlab = "tempo", ylab = "y")

> par(mfrow = c(2, 2))

> for (i in w) {

+ ts.plot(mld.sim(100, 1, i, 25), gpars = g, main = paste("V/W=",

+ 1/i))

+ }

V/W= 20

tempo

y

0 20 40 60 80 100

21 23

25 27

V/W= 2

tempo

y

0 20 40 60 80 100

18 22

26 30

V/W= 0.2

tempo

y

0 20 40 60 80 100

10 20

30 40

Figura 7.1: 100 valores simulados do modelo polinomial de 1a ordem com (a) V/W = 20, (b) V/W = 2, (c) V/W = 0, 2.

7.8. MONITORAMENTO 103

> y = Nile

> n = length(y)

> res = mld(y, V = rep(var(y), n), W = rep(50, n), m0 = 1000, C0 = 1000)

> plot(y, xlab = "Anos", ylab = "Mediç~oes", type = "p")

> lines(res$m, col = 2)

> lines(res$m - 2 * sqrt(res$C), col = 2, lty = 1)

> lines(res$m + 2 * sqrt(res$C), col = 2, lty = 1)

> res = mld(y, V = rep(var(y), n), W = rep(0.05, n), m0 = 1000,

+ C0 = 1000)

> lines(res$m, col = 4)

> legend(1940, 1350, c("obs", "W=50", "W=.05"), col = c(1, 2, 4),

+ bty = "n")

Anos

M ed

iç õe

s

1880 1900 1920 1940 1960

60 0

80 0

10 00

12 00

14 00

obs W=50 W=.05

Figura 7.2:

104 CAPÍTULO 7. MODELOS LINEARES DINÂMICOS

Desconto 0.98

Time

1880 1900 1920 1940 1960

60 0

80 0

10 00

12 00

14 00

70

60 0

80 0

10 00

12 00

14 00

Serie original

1880 1920 1960

60 0

10 00

14 00

1880 1920 1960

60 0

10 00

14 00

obs desconto=.98 desconto=.70

Figura 7.3:

Apêndice A

Lista de Distribuições

Neste apêndice são listadas as distribuições de probabilidade utilizadas no texto para

facilidade de referência. São apresentadas suas funções de (densidade) de probabili-

dade além da média e variância. Uma revisão exaustiva de distribuições de probabili-

dades pode ser encontrada em Johnson et al. (1994), Johnson et al. (1995) e Johnson

et al. (1992).

A.1 Distribuição Normal

X tem distribuição normal com parâmetros µ e σ2, denotando-se X ∼ N(µ, σ2), se sua função de densidade é dada por

p(x|µ, σ2) = (2πσ2)−1/2 exp[−(x− µ)2/2σ2], −∞ < x < ∞,

para −∞ < µ < ∞ e σ2 > 0. Quando µ = 0 e σ2 = 1 a distribuição é chamada normal padrão. A distribuição log-normal é definida como a distribuição de eX .

No caso vetorial, X = (X1, . . . , Xp) tem distribuição normal multivariada com

vetor de médias µ e matriz de variância-covariância Σ, denotando-se X ∼ N(µ,Σ) se sua função de densidade é dada por

p(x|µ,Σ) = (2π)−p/2|Σ|−1/2 exp[−(x− µ)′Σ−1(x− µ)/2]

para µ ∈ Rp e Σ positiva-definida.

A.2 Distribuição Gama

X tem distribuição Gama com parâmetros α e β, denotando-se X ∼ Ga(α, β), se sua função de densidade é dada por

p(x|α, β) = β α

Γ(α) xα−1e−βx, x > 0,

para α, β > 0.

E(X) = α/β e V (X) = α/β2.

105

106 APÊNDICE A. LISTA DE DISTRIBUIÇÕES

Casos particulares da distribuição Gama são a distribuição de Erlang, Ga(α, 1), a

distribuição exponencial, Ga(1, β), e a distribuição qui-quadrado com ν graus de

liberdade, Ga(ν/2, 1/2).

A.3 Distribuição Wishart

Diz-se que uma matriz aleatória Ω (n × n) segue uma distribuição Wishart com parâmetro Σ e ν graus de liberdade, denotando-se Ω ∼ W (Σ, ν), se sua função de densidade é dada por,

p(Ω|Σ, ν) ∝ |Ω|(ν−n−1)/2 exp(−(1/2)tr(ΣΩ))

sendo ν ≥ n, Σ positiva-definida e tr(A) indica o traço de uma matriz A. Uma propriedade útil é que AΩA′ ∼ W (AΣA′, ν).

A.4 Distribuição Gama Inversa

X tem distribuição Gama Inversa com parâmetros α e β, denotando-se

X ∼ GI(α, β), se sua função de densidade é dada por

p(x|α, β) = β α

Γ(α) x−(α+1)e−β/x, x > 0,

para α, β > 0.

E(X) = β

α− 1 e V (X) = β2

(α− 1)2(α− 2) .

Não é dif́ıcil verificar que esta é a distribuição de 1/X quando X ∼ Ga(α, β).

A.5 Distribuição Wishart Invertida

Diz-se que uma matriz aleatória Ω (n× n) segue uma distribuição Wishart-Invertida com parâmetro Σ e ν graus de liberdade, denotando-se Ω ∼ WI(Σ, ν) se sua função de densidade é dada por,

p(Ω|Σ, ν) ∝ |Ω|−(ν+n+1)/2 exp(−(1/2)tr(ΣΩ))

sendo ν ≥ n, Σ positiva-definida e tr(A) indica o traço de uma matriz A. Não é dif́ıcil verificar que Ω−1 ∼ W (Σ, ν). Outra propriedade é que AΩA′ ∼ WI(AΣA′, ν).

A.6 Distribuição Beta

X tem distribuição Beta com parâmetros α e β, denotando-se X ∼ Be(α, β), se sua função de densidade é dada por

p(x|α, β) = Γ(α+ β) Γ(α)Γ(β)

xα−1(1− x)β−1, 0 < x < 1,

A.7. DISTRIBUIÇÃO DE DIRICHLET 107

para α, β > 0.

E(X) = α

α+ β e V (X) =

αβ

(α+ β)2(α+ β + 1) .

A.7 Distribuição de Dirichlet

O vetor aleatório X = (X1, . . . , Xk) tem distribuição de Dirichlet com parâmetros

α1, . . . , αk, denotada por Dk(α1, . . . , αk) se sua função de densidade conjunta é dada

por

p(x|α1, . . . , αk) = Γ(α0)

Γ(α1), . . . ,Γ(αk) xα1−11 . . . x

αk−1 k ,

k ∑

i=1

xi = 1,

para α1, . . . , αk > 0 e α0 = ∑k

i=1 αi.

E(Xi) = αi α0

, V (Xi) = (α0 − αi)αi α20(α0 + 1)

, e Cov(Xi, Xj) = − αiαj

α20(α0 + 1)

Note que a distribuição Beta é obtida como caso particular para k = 2.

A.8 Distribuição t de Student

X tem distribuição t de Student (ou simplesmente t) com média µ, parâmetro de escala

σ e ν graus de liberdade, denotando-se X ∼ tν(µ, σ2), se sua função de densidade é dada por

p(x|ν, µ, σ2) = Γ((ν + 1)/2)ν ν/2

Γ(ν/2) √ πσ

[

ν + (x− µ)2

σ2

]−(ν+1)/2

, x ∈ R,

para ν > 0, µ ∈ R e σ2 > 0.

E(X) = µ, para ν > 1 e V (X) = ν

ν − 2 , para ν > 2.

Um caso particular da distribuição t é a distribuição de Cauchy, denotada por

C(µ, σ2), que corresponde a ν = 1.

A.9 Distribuição F de Fisher

X tem distribuição F com ν1 e ν2 graus de liberdade, denotando-se X ∼ F (ν1, ν2), se sua função de densidade é dada por

p(x|ν1, ν2) = Γ((ν1 + ν2)/2)

Γ(ν1/2)Γ(ν2/2) ν ν1/2 1 ν

ν2/2 2 x

ν1/2−1(ν2 + ν1x) −(ν1+ν2)/2

x > 0, e para ν1, ν2 > 0.

E(X) = ν2

ν2 − 2 , para ν2 > 2 e V (X) =

2ν22(ν1 + ν2 − 2) ν1(ν2 − 4)(ν2 − 2)2

, para ν2 > 4.

108 APÊNDICE A. LISTA DE DISTRIBUIÇÕES

A.10 Distribuição Binomial

X tem distribuição binomial com parâmetros n e p, denotando-se X ∼ bin(n, p), se sua função de probabilidade é dada por

p(x|n, p) = (

n

x

)

px(1− p)n−x, x = 0, . . . , n

para n ≥ 1 e 0 < p < 1.

E(X) = np e V (X) = np(1− p)

e um caso particular é a distribuição de Bernoulli com n = 1.

A.11 Distribuição Multinomial

O vetor aleatório X = (X1, . . . , Xk) tem distribuição multinomial com parâmetros n

e probabilidades θ1, . . . , θk, denotada por Mk(n, θ1, . . . , θk) se sua função de probabi-

lidade conjunta é dada por

p(x|θ1, . . . , θk) = n!

x1!, . . . , xk! θx11 , . . . , θ

xk k , xi = 0, . . . , n,

k ∑

i=1

xi = n,

para 0 < θi < 1 e ∑k

i=1 θi = 1. Note que a distribuição binomial é um caso especial

da multinomial quando k = 2. Além disso, a distribuição marginal de cada Xi é

binomial com parâmetros n e θi e

E(Xi) = nθi, V (Xi) = nθi(1− θi), e Cov(Xi, Xj) = −nθiθj .

A.12 Distribuição de Poisson

X tem distribuição de Poisson com parâmetro θ, denotando-se X ∼ Poisson(θ), se sua função de probabilidade é dada por

p(x|θ) = θ xe−θ

x! , x = 0, 1, . . .

para θ > 0.

E(X) = V (X) = θ.

A.13 Distribuição Binomial Negativa

X tem distribuição de binomial negativa com parâmetros r e p, denotando-se X ∼ BN(r, p), se sua função de probabilidade é dada por

p(x|r, p) = (

r + x− 1 x

)

pr(1− p)x, x = 0, 1, . . .

A.13. DISTRIBUIÇÃO BINOMIAL NEGATIVA 109

para r ≥ 1 e 0 < p < 1.

E(X) = r(1− p)/p e V (X) = r(1− p)/p2.

Um caso particular é quando r = 1 e neste caso diz-se que X tem distribuição geo-

métrica com parâmetro p.

Referências

Bauwens, L., Lubrano, M. & Richard, J. (1999). Bayesian Inference in Dynamic

Econometric Models. Oxford University Press.

Box, G. E. P. & Jenkins, G. M. (1970). Time Series Analysis, Forecasting and

Control. Holden-Day, San Francisco, California.

Box, G. E. P., Jenkins, G. M. & Reinsel, G. C. (1994). Time Series Analysis:

Forecasting and Control (Third ed.). Englewood Cliffs NJ: Prentice-Hall.

Brockwell, P. & Davis, R. (1991). Time Series: Theory and Methods (2nd ed.).

New York: Springer-Verlag.

Burnham, K. P. & Anderson, D. R. (1998). Model Selection and Inference: A

Practical Information-Theoretic Approach. Springer: New York.

Diggle, P. (1990). Time Series: A Biostatistical Introduction. Oxford University

Press: New York.

Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates

of the variance of United Kingdom inflation. Econometrica 50, 987–1007.

Franses, P. H. (1998). Time Series Models for Business and Economic Forecasting.

Cambridge University Press.

Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press.

Johnson, N. L., Kotz, S. & Balakrishnan, N. (1994). Continuous Univariate Dis-

tributions (2nd ed.), Volume 1. John Wiley, New York.

Johnson, N. L., Kotz, S. & Balakrishnan, N. (1995). Continuous Univariate Dis-

tributions (2nd ed.), Volume 2. John Wiley, New York.

Johnson, N. L., Kotz, S. & Kemp, A. W. (1992). Univariate Discrete Distributions

(2nd ed.). John Wiley, New York.

Kendall, M. G., Stuart, A. & Ord, J. K. (1983). Advanced theory of statistics (4th

ed.), Volume 3. Griffin: London.

Pole, A., West, M. & Harrison, J. (1994). Applied Bayesian Forecasting and Time

Series Analysis. Texts in Statistical Sciences. Chapman & Hall.

Priestley, M. B. (1981). Spectral Analysis and Time Series. London: Academic

Press.

Taylor, S. (1986). Modelling Financial Time Series. Wiley.

110

References. 111

Tsay, R. S. (2002). Analysis of Financial Time Series. Wiley.

West, M. & Harrison, P. J. (1997). Bayesian Forecasting and Dynamic Models.

Springer Verlag, New York.

Até o momento nenhum comentário
Esta é apenas uma pré-visualização
3 mostrados em 115 páginas